format: text

text:
# Asymptotic Stability of Traveling Wave Solutions for Perturbations with Algebraic Decay

## 1. Introduction

The topic of this note is the class of regularized scalar conservation laws in one spatial dimension

$$u_t+\beta u_{xxx}\alpha u_{xx}+g(u)_x=0$$

where subscripts denote partial derivatives. The regularization is due to the presence of viscous terms ($`\alpha >0`$) and dispersive terms ($`\beta 0`$). The case $`g(u)=u^2/2`$ is typical and has received much attention. If $`\alpha >0=\beta `$, this is known as Burgers equation. If $`\alpha =0<\beta `$, this is essentially the Korteweg - deVries equation. The case $`\alpha ,\beta >0`$ thus is referred to with a canonical all-Dutch name; it also been studied extensively, as has been the case of general $`g`$. I want to study the stability of traveling wave solutions of the form $`u(x,t)=\varphi (xct)`$ with respect to perturbations of the initial data, in the cases $`\alpha >0=\beta `$ and $`\alpha ,\beta >0`$. Here $`c`$ is the speed at which the wave profile $`\varphi `$ travels to the right (if $`c>0`$). Only monotone wave profiles $`\varphi `$ will be considered. The question is whether the solution approaches a traveling wave in some sense. A natural setting for this is a spatial coordinate system that moves along with the expected wave profile at speed $`c`$.

Since all shifted wave profiles also give rise to traveling wave solutions, one can only expect that the solution will converge to some shifted profile $`\varphi (xcth)`$. Since the quantity $`_{}\left(u(x,t)v(x,t)\right)𝑑x`$ is independent of $`t`$ for any pair of solutions $`u`$ and $`v`$ for which it is finite, the shift $`h`$ must be such that $`_{}\left(u(x,0)\varphi (xh)\right)𝑑x=0`$. It is easy to see that the quantity on the left hand side is an affine function of $`h`$, and thus $`h`$ can be determined explicitly and can be considered a known quantity.

Stability in this sense was first studied in , for the viscous case $`\beta =0`$. These authors noted that one cannot expect a rate of convergence that holds for all classes of perturbations. However, in 1976, Sattinger showed in  that an exponential rate of decay holds in a moving coordinate system if the perturbation of the initial value decays exponentially in space. Since this class of perturbations is somewhat restrictive, one may ask the question what the consequences of algebraic decay of the initial perturbation are. It was shown in ,  and  that in this case algebraic decay for the initial perturbation can be ”traded in” for some temporal algebraic decay. Heuristically, the equation for the perturbation behaves like $`w_tw_x=0`$ near $`x=\pm \mathrm{}`$ in this situation, that is, the solution near $`\mathrm{}`$ behaves like $`w(x,t)w(x+t,0)`$. Therefore, if $`w(x,0)e^x`$, then $`w(x,t)e^te^x`$ near $`x=\mathrm{}`$ as $`t\mathrm{}`$ (exponential decay in exponentially weighted norms), and if $`w(x,0)x^k`$, then $`w(x,t)x^{mk}t^m`$ near $`x=\mathrm{}`$ as $`t\mathrm{}`$ (algebraic decay with a weaker algebraic weight). For finite $`x`$, diffusion dominates and leads to exponential decay. Similar results were shown for case $`\beta >0`$ in , assuming the wave profile $`\varphi `$ is monotone. The proofs rely on a detailed study of the spectrum of the linearized problem or on ad hoc energy estimates. A more refined analysis relying on Green’s function estimates is carried out in .

The goal of this note is to show that this ”trade-off” follows whenever exponential decay holds in a setting with spatial exponential weights together with (simple) stability in a setting without weights. Showing these two conditions is often easier than analyzing the full problem in a setting with polynomials weights. On the other hand, for viscous conservation laws, this interpolation argument can only be applied in the ”totally compressive case” where all characteristics run into the shock for the inviscid problem. This is always true for the scalar case, but systems of regularized conservation laws from physical situations usually do not fall in this category; see  and  for the much more complicated theory for this case. On the other hand, the argument given in this note is not restricted to the case of one space dimension.

The paper is organized as follows. In the next section, I show a general interpolation result for linearized problems. In section 3, this is applied to study the generalized Burgers equation. The main result is a sharp decay result in spaces with algebraic spatial decay. In section 4, a similar result is shown for the generalized Korteweg-deVries - Burgers equation. Appendices A and B contain results for the corresponding linearized problems. Appendix C contains a simple integral inequality that is used for the passage from linear to nonlinear stability.

Here is some notation that is used throughout the paper. Let $`\mathrm{\Omega }^n`$ be an unbounded measurable set. For $`1p<\mathrm{},k>0,\rho >0`$ let us define the function spaces

$`L^{p,k}`$ $`=\{uL^p(\mathrm{\Omega })|{\displaystyle _\mathrm{\Omega }}|u(x)|^p(1+|x|)^{kp}𝑑x=u_{p,k}^p<\mathrm{}\}`$
$`L^{p,\rho }`$ $`=\{uL^p(\mathrm{\Omega })|{\displaystyle _\mathrm{\Omega }}|u(x)|^pe^{\rho p|x|}𝑑x=u_{p,\rho }^p<\mathrm{}\}`$
with their natural norms. The usual modifications are made to define $`L^{\mathrm{},k}`$ and $`L^{\mathrm{},\rho }`$. The norms on the unweighted $`L^p`$-spaces are denoted by $`_p`$. The set $`\mathrm{\Omega }`$ is specified to be the real line in sections 3 and 4 and in appendices A and B, thus it does not appear further in the notation. Constants are denoted by the same letter $`C`$ whose value may change from line to line, depending only on values that can be expressed in terms of quantities mentioned in the assumptions of a result. If constants have indices, their values remain constant throughout a proof.

## 2. A Linear Interpolation Result

Suppose we are given an operator $`S`$, not necessarily linear, which maps the space $`L^p(\mathrm{\Omega })`$ into itself and also $`L^{p,\rho }(\mathrm{\Omega })`$ into itself, with estimates

$`S(u)S(v)_p`$ $`C_0uv_p\text{for all}u,vL^p(\mathrm{\Omega })`$ $`2.1a`$$`2.1b`$
$`S(u)_{p,\rho }`$ $`C_0e^tu_{p,\rho }\text{for all}uL^{p,\rho }(\mathrm{\Omega })`$
for some constants $`\rho ,C_0,t>0`$. The main result of this section says that $`S`$ also maps $`L^{p,k}`$ into $`L^{p,l}`$ for $`0<l<k`$ and gives an estimate for this mapping.

I shall prove this in detail for the case $`p<\mathrm{}`$. A different proof will be given for the case where $`S`$ is linear and $`p=\mathrm{}`$.

###### Theorem 2.1

Let $`1p<\mathrm{}`$. Under the above conditions, $`S`$ maps $`L^{p,k}`$ into $`L^{p,l}`$ for all $`0<l<k`$, and there exists $`C_1>0`$, depending only on $`p`$, such that for all $`vL^{p,k}`$

$$Sv_{p,l}C_0C_1\rho ^{kl}(1+t)^{lk}v_{p,k}.$$
$`2.2`$
###### Demonstration Proof

Let us assume $`\rho =1`$ and continue to write $`_{p,\rho }`$ for the corresponding norm with exponential weight. Fix $`p[1,\mathrm{})`$ and define for $`r0`$

$$m_p(r)=\frac{r}{\left(1+r^{\frac{p}{p1}}\right)^{\frac{p1}{p}}}$$

if $`p>1`$ and $`m_1(r)=\mathrm{min}(1,r)`$. Clearly $`m_1(r)2m_p(r)2m_1(r)`$ for all $`r0`$ and all $`p`$. For $`s`$ and $`uL^p(\mathrm{\Omega })`$ define the functional

$$K(s,u)=\left(_\mathrm{\Omega }(|u(x)|m_p(e^{s+|x|}))^p𝑑x\right)^{1/p}.$$

This is clearly an equivalent norm on $`L^p`$ and in fact a modified K-functional (), namely

$$K(s,u)=\underset{vL^{p,\rho }}{inf}\left(uv_p^p+e^{sp}v_{p,\rho }^p\right)^{1/p}.$$
$`2.3`$
Indeed,

$`\underset{vL^{p,\rho }}{inf}\left(uv_p^p+e^{sp}v_{p,\rho }^p\right)`$ $`=\underset{vL^{p,\rho }}{inf}{\displaystyle _\mathrm{\Omega }}\left(|u(x)v(x)|^p+e^{sp+|x|p}|v(x)|^p\right)𝑑x`$
$`{\displaystyle _\mathrm{\Omega }}\underset{\zeta }{inf}(|u(x)\zeta )|^p+e^{sp+|x|p}|\zeta |^p)dx`$
$`={\displaystyle _\mathrm{\Omega }}\left(|u(x)|m_p(e^s+|x|)\right)^p𝑑x`$
$`=\left(uv_0_p^p+e^{sp}v_0_{p,\rho }^p\right).`$
Here

$$v_0(x)=\frac{u(x)}{1+e^{(s+|x|)\frac{p}{p1}}}$$

for $`p>1`$. For $`p=1`$ one sets $`v_0(x)=u(x)`$ for $`s+|x|0`$ and $`v_0(x)=0`$ otherwise, and the last equation is again valid. Then $`v_0L^{p,\rho }`$, and (2.3) follows. The definition shows immediately that $`K(,u)`$ is differentiable, non-decreasing, and bounded above by $`u_p`$ for all $`p`$. Also,

$$2^p_\mathrm{\Omega }|u(x)|^pm_1(e^{s+|x|})^p𝑑xK(s,u)^p_\mathrm{\Omega }|u(x)|^pm_1(e^{s+|x|})^p𝑑x.$$

Next fix also $`k>0`$ and set

$$h_k(s)=\{\begin{array}{cc}& e^s(s0)\hfill \\ & (1s)^{kp1}(s<0)\hfill \end{array}.$$
$`2.4`$
and

$$u_{}^p=_{\mathrm{}}^{\mathrm{}}K(s,u)^ph_k(s)𝑑s$$

whenever this quantity is finite. The next claim is that

$$L^{p,k}=\left\{u\right|u_{}^p<\mathrm{}\}$$

and that $`_{}`$ is an equivalent norm on this space. Indeed,

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}K(s,u)^ph_k(s)𝑑s`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle _\mathrm{\Omega }}|u(x)|^p\mathrm{min}(1,e^{sp+|x|p})h_k(s)𝑑x𝑑s`$
$`={\displaystyle _\mathrm{\Omega }}|u(x)|^p\left({\displaystyle _{\mathrm{}}^{|x|}}e^{sp+|x|p}(1s)^{kp1}𝑑s+{\displaystyle _{|x|}^0}(1s)^{kp1}𝑑s+{\displaystyle _0^{\mathrm{}}}e^{sp}𝑑s\right)𝑑x`$
$`{\displaystyle _\mathrm{\Omega }}|u(x)|^pC\left(1+|x|^{kp}\right)𝑑xCu_{p,k}^p.`$
Reversely,

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}K(s,u)^ph_k(s)𝑑s`$ $`2^p{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle _\mathrm{\Omega }}|u(x)|^p\mathrm{min}(1,e^{sp+|x|p})h_k(s)𝑑x𝑑s`$
$`2^p{\displaystyle _\mathrm{\Omega }}|u(x)|^p{\displaystyle _{|x|}^0}(1s)^{kp1}𝑑s𝑑x`$
$`Cu_{p,k}^p.`$
Let now $`S`$ be an operator satisfying (2.1a,b). Then

$`K(s,Su)^p`$ $`=\underset{vL^{p,\rho }}{inf}\left(Suv_p^p+e^sv_{p,\rho }^p\right)`$
$`\underset{vL^{p,\rho }}{inf}\left(SuSv_p^p+e^sSv_{p,\rho }^p\right)`$
$`C_0\underset{vL^{p,\rho }}{inf}\left(uv_p^p+e^{st}v_{p,\rho }^p\right)`$
$`=C_0K(st,u)^p.`$
Let $`0<l<k`$, let $`uL^{p,\rho }`$, and set $`H_l(r)=_r^{\mathrm{}}h_l(\tau )𝑑\tau `$ and $`k(s)=\frac{d}{ds}K(s,u)^p0`$. Then

$`Su_{p,l}^p`$ $`C{\displaystyle _{\mathrm{}}^{\mathrm{}}}K(s,Su)^ph_l(s)𝑑s`$
$`CC_0{\displaystyle _{\mathrm{}}^{\mathrm{}}}K(st,u)^ph_l(s)𝑑s`$
$`=CC_0{\displaystyle _{\mathrm{}}^{\mathrm{}}}k(s)H_l(s+t)𝑑s.`$
An elementary calculation shows that $`H_l(s+t)CH_k(s)(1+t)^{lk}`$ for all $`s`$ and $`t`$. One can therefore estimate further

$$Su_{p,l}^pCC_0(1+t)^{lk}_{\mathrm{}}^{\mathrm{}}k(s)H_k(s)𝑑s=C_1C_0(1+t)^{lk}u_{p,k}.$$

This proves the theorem in the case $`\rho =1`$. The general case follows by a scaling argument.

The proof can be modified to extend to the case $`p=\mathrm{}`$. However, I prefer to give an alternative proof in this case. It extends to subspaces of $`L^{\mathrm{}}`$ that are closed under multiplication with smooth functions that are bounded together with their derivatives. Such spaces include $`X=BC^m(\mathrm{\Omega })`$ and $`X=L^{\mathrm{}}(\mathrm{\Omega })UC(\mathrm{\Omega })`$, where $`BC^m`$ is the space of $`m`$-times differentiable functions with bounded derivatives and $`UC`$ is the set of uniformly continuous functions on $`\mathrm{\Omega }`$. The result is formulated for the cases $`X=L^{\mathrm{}}(\mathrm{\Omega })`$, $`X=BC^0(\mathrm{\Omega })`$, and $`X=L^{\mathrm{}}(\mathrm{\Omega })UC(\mathrm{\Omega })`$, equipped with the supremum norm. Accordingly let $`X_k=XL^{\mathrm{},k}`$ and $`X_\rho =XL^{\mathrm{},\rho }`$, equipped with their natural norms $`_{\mathrm{},k}`$ and $`_{\mathrm{},\rho }`$.

###### Theorem 2.2

Let $`S:XX`$, $`S:X_\rho X_\rho `$ be a linear operator for which (2.1a,b) holds. Then $`S`$ maps $`X_k`$ into $`X_l`$ for all $`0<l<k`$, and there exists a constant $`C_1>0`$ depending on $`k`$ and $`\rho `$ such that for all $`vX_k`$

$$Sv_{\mathrm{},l}\frac{C_0C_1}{(1+t)^{kl}}v_{\mathrm{},k}.$$
$`2.5`$
###### Demonstration Proof

Let $`vL^{\mathrm{},k}`$ be given. Thus $`(1+|x|)^k|v(x)|A`$ almost everywhere for some smallest constant $`A`$. Let $`w=Sv`$. The goal is to show that $`(1+t)^{kl}(1+|x|)^k|w(x)|C_0C_1A`$ almost everywhere for some constant $`C_1`$. Let $`R0`$, to be chosen later. Choose $`\phi C^{\mathrm{}}()`$ with $`\phi (r)=0`$ for $`r0`$, $`\phi (r)=1`$ for $`r1`$, and $`\phi ^{}0`$. Set $`v_1(x)=v(x)\phi (|x|R)`$ and $`v_2=vv_1`$. Then $`v_1_{\mathrm{}}A(1+R)^k`$ and $`v_2_{\mathrm{},\rho }C_0C(1+R)^ke^{\rho R}`$ with $`C=C(k,\rho )0`$. Thus

$$Sv_1_{\mathrm{}}C_0(1+R)^kA\text{and}Sv_2_{\mathrm{},\rho }C_0C(1+R)^ke^{\rho Rt}A.$$

Let $`x\mathrm{\Omega }`$. Then

$`(1+t)^{kl}(1+|x|)^l|w(x)|`$ $`(1+t)^{kl}(1+|x|)^l\left(|Sv_1(x)|+|Sv_2(x)|\right)`$
$`C_0(1+t)^{kl}(1+|x|)^l(1+R)^kA(1+e^{\rho R\rho |x|t}).`$
Set $`\sigma =\mathrm{min}({\displaystyle \frac{l}{k}},{\displaystyle \frac{kl}{\rho k}})`$ and choose $`R=\sigma t^{\frac{kl}{k}}|x|^{\frac{l}{k}}\frac{t}{\rho }+|x|`$. One can estimate further

$$(1+t)^{kl}(1+|x|)^l|w(x)|2C_0(1+t)^{kl}(1+|x|)^l(1+R)^kAC_0C_1A.$$

This proves the theorem.

At first glance, it is surprising that the algebraic decay estimates are independent of the exponential spatial weight that appears in the assumptions. The scaling argument used at the end of the proof explains this phenomenon and shows where the constant $`\rho `$ reappears in the result. The special case $`\mathrm{\Omega }=(\mathrm{},0]`$ with the right shift operator $`S_tv(x)=v(xt)`$ for $`x0`$ shows that all estimates in the two theorems are sharp, up to the values of the constants. Indeed, the use of K-functionals in the proof of theorem 2.1 makes the argument resemble the direct proof for this special case.

## 3. Scalar Viscous Conservation Laws

In this section scalar viscous conservation laws of the following form are considered:

$$u_tu_{xx}+f(u)_x=0(x,t>0).$$
$`3.1`$
Here $`f:`$ is a $`C^2`$-function with uniformly bounded derivatives. A traveling wave solution $`u(x,t)=\varphi (xct)`$ with speed $`c`$ and limiting behavior $`\varphi (r)\varphi _\pm `$ as $`r\pm \mathrm{}`$ is easily seen to exist if and only if $`c`$ equals the slope of the line segment connecting the points $`(\varphi _+,f(\varphi _+))`$ and $`(\varphi _{},f(\varphi _{}))`$ and the graph of $`f`$ lies entirely above or below this line segment. Accordingly, the wave profile $`\varphi `$ is decreasing or increasing. After rescaling the dependent variable $`u`$, adding a linear function to $`f`$, and changing to a moving coordinate system, one can assume that $`c=0,\varphi _{}=1,\varphi _+=0`$. The prototypical example is $`f(r)=r^2r`$, with the wave profile $`\varphi (r)=\left(1+e^r\right)^1`$.

The focus is the convergence of solutions of (3.1) to some translate $`\varphi (h)`$ as $`t\mathrm{}`$, where $`h`$ is known. Let us therefore assume that $`h=0`$. It was shown in  that this convergence in the uniform sense follows if $`f`$ is uniformly convex and $`_0^{\mathrm{}}|u(x,0)|𝑑x+_{\mathrm{}}^0|u(x,0)1|𝑑x`$ is finite. If $`f`$ is merely $`C^2`$-smooth but not necessarily convex, the same conclusion holds in the $`L^1`$-sense even for the case of quasilinear diffusion (). In the seminal paper , it was shown that convergence at an exponential rate holds in spaces with exponentially weighted norms, namely

$$u(,t)\varphi _{\mathrm{},ϵ}=O(e^{\delta t})$$
$`3.2`$
for some $`ϵ,\delta >0`$, if the quantity on the left hand side is sufficiently small for $`t=0`$. This holds for arbitrary non-convex $`f`$, assuming only that

$$f^{}(0)0f^{}(1).$$
$`3.3`$
In several recent papers, the stability of wave profiles in spaces with polynomial weights was discussed. Assuming only (3.2) it was shown in  that if

$$u(,0)\varphi _{\mathrm{},k+m}=\delta $$
$`3.4a`$
is sufficiently small, then

$$u(,t)\varphi ()_{\mathrm{},k}C(1+t)^{m/2}\delta $$
$`3.4b`$
for integers $`k,m`$ satisfying $`k1,\mathrm{\hspace{0.17em}2}mk+1`$, or $`k=1,m2`$. A comparable result for this situation from  and  assumes that

$$\mathrm{\Psi }_{2,\alpha }<\mathrm{}$$
$`3.5a`$
where $`\mathrm{\Psi }(x)=_{\mathrm{}}^x\left(u(z,0)\varphi (z)\right)𝑑z`$, and that a suitable unweighted $`L^2`$ \- norm of $`u(,0)\varphi `$ is small. The conclusion then is essentially that

$$u(,t)\varphi _{\mathrm{}}C(1+t)^\alpha .$$
$`3.5`$
Here $`\alpha >0`$ is arbitrary. Both results show a trade-off between the spatial decay of the initial data and the temporal decay of the solution. A much more detailed and general result in  implies that spatial decay of the antiderivative $`\mathrm{\Psi }`$ dominates the temporal decay of $`u(,t)\varphi `$ for finite $`x`$, while the spatial decay of $`u(,0)\varphi `$ dominates the temporal and spatial decay as both $`x`$ and $`t`$ go to $`\mathrm{}`$, with the canonical trade-off.

The main result for this situation assumes also (3.3).

###### Theorem 3.1

Let $`k>1,\mathrm{\hspace{0.17em}0}<m<k`$ be real numbers, and set $`\mathrm{\Psi }(x)=_{\mathrm{}}^x\left(u(z,0)\varphi (z)\right)𝑑z`$. There exists a constant $`C_0`$ such that if

$$\mathrm{\Psi }_{\mathrm{},k}=ϵ(u_0)$$
$`3.6a`$
is sufficiently small, then the solution $`u`$ exists for all $`t>0`$, and for all $`0<m<k`$

$$u(,t)\varphi _{\mathrm{},m}C_0(1+t)^{mk}ϵ(u_0).$$
$`3.6b`$
###### Demonstration Proof

The proof follows a pattern which will be repeated in the next section. A formal linearization of the problem is introduced (step 0), a function space setting is defined, and a solution $`u`$ is produced, using results about the linearized equation from Appendix A (step 1). After specifying the short time behavior of the solution (step 2), suitable a priori estimates are shown with the aid of Lemma C.1, which finishes the proof.

Step 0. Let $`u`$ be a solution of (3.1), and set $`v(x,t)=_{\mathrm{}}^x\left(u(z,t)\varphi (z)\right)𝑑z`$. Thus $`v(x,0)=\mathrm{\Psi }(x)`$, and $`v`$ solves the equation

$$v_t(x,t)v_{xx}(x,t)=f(\varphi (x))f(\varphi (x)+v_x(x,t))$$
$`3.7`$
or

$$v_t(,t)+Lv(,t)=F(v)(,t)$$

where formally

$`Ly(x)`$ $`=y_{xx}(x)+f^{}(\varphi (x))y_x(x)`$ $`3.8`$$`3.9`$
$`F(y)(x)`$ $`=f(\varphi (x))+f^{}(\varphi (x))y_x(x)f(\varphi (x)+y_x(x)=f^{\prime \prime }(\zeta (x))y_x^2(x).`$
Here $`\zeta (x)`$ is a number between $`\varphi (x)`$ and $`\varphi (x)+y_x(x)`$.

Step 1. Recall the notation and the results from Appendix A: $`X=L^{\mathrm{}}()C^0()`$, $`X_0=\{yX|y\text{is uniformly continuous}\}`$ are Banach spaces, both equipped with the supremum norm. The operator $`L`$ acting in $`X`$ is defined in (3.8), and the restriction of $`L`$ to $`X_0`$ is also denoted by $`L`$. Then $`L`$ generates an analytic semigroup $`(S(t))_{t0}`$ in $`X_0`$ and in various weighted spaces, especially in the spaces $`X_\rho `$ and $`X_k`$. From the properties of $`f`$ and the definition of $`F`$ in (3.8), one immediately deduces the estimates

$`F(y_1)F(y_2)_{\mathrm{},k}`$ $`Cy_{1,x}y_{2,x}_{\mathrm{},k}`$ $`3.10a`$$`3.10b`$
$`F(y)_{\mathrm{},k}`$ $`Cy_x_{\mathrm{},k/2}^2`$
for any $`k0`$ and suitable $`y,y_1,y_2`$. Suppose that $`\mathrm{\Psi }X_kX_0`$ is given and satisfies the assumptions of Theorem 3.1. Let us seek a continuous $`X_0`$-valued solution $`v(,t)`$ of the integral equation

$$v(,t)=S(t)\mathrm{\Psi }+_0^tS(ts)F(v)(,s)𝑑s.$$
$`3.11`$
A contraction argument, using (3.10a), produces such a solution on some finite time interval, and the solution satisfies (3.7) in the classical sense if it is sufficiently smooth. Due to the global Lipschitz property in (3.10a), this solution exists for all times. The existence and uniqueness arguments hold in all spaces $`X_k`$, $`0<km`$, and therefore $`v(,t)`$ belongs to all these spaces for all $`t`$. Set $`w(x,t)=u(x,t)\varphi (x)=v_x(x,t)`$, and let us write $`N(t)=\mathrm{max}(1,t^{1/2})`$.

Step 2. Let us now characterize the behavior of $`w`$ for $`0<t1`$ in more detail. The goal is to show the estimate

$$w(,t)_{\mathrm{},m}CN(t)ϵ(u_0)(0<t1)$$
$`3.12`$
for all $`0<mk`$. For this purpose consider the integral equation for $`w`$ which is obtained by differentiating (3.1), namely

$$w(,t)=[S(t)\mathrm{\Psi }]_x+_0^t[S(ts)F(v)(,s)]_x𝑑s.$$
$`3.13`$
For the “free term” $`[S(t)\mathrm{\Psi }]_x`$ in (3.13), (3.12) is just estimate (A.6). Using (3.12), (A.6), and (3.10b) one then derives the inequality

$$\frac{w(,t)_{\mathrm{},m}}{N(t)}Cϵ(u_0)+_0^tC\frac{N(ts)N(s)}{N(t)}\frac{w(,s)_{\mathrm{},m}}{N(s)}𝑑s.$$
$`3.14`$
A standard argument for linear integral inequalities implies that $`\frac{w(,t)_{\mathrm{},m}}{N(t)}C`$ on $`(0,1]`$, i.e. (3.12).

Step 3. Finally, let us prove estimate (3.6b). The first thing to notice is again that the estimate (3.6b) in Theorem 3.1 holds for the “free” term $`[S(t)\mathrm{\Psi }]_x`$, since $`[S(t)\mathrm{\Psi }]_x_{\mathrm{},m}CS(t1/2)\mathrm{\Psi }_{\mathrm{},m}C(1+t)^{mk}ϵ(u_0)`$. Consider first the special case $`m=k/2`$ and thus $`mk=m`$. Define the quantity

$$\gamma (t)=\underset{1st}{sup}s^mw(,s)_{\mathrm{},m}$$
$`3.15`$
for $`t1`$. Then $`\gamma `$ is continuous, $`\gamma (1)=w(,1)_mC_0ϵ(u_0)`$ for some constant $`C_0`$, and $`F(v)(,s)_{\mathrm{},m}Cs^k\gamma ^2(s)`$ for all $`s`$. Then (3.13) implies for $`t1`$

$`w(,t)_{\mathrm{},m}`$ $`Ct^mϵ(u_0)+{\displaystyle _0^1}C(1+ts)^mN(ts)w(,s)_{\mathrm{},m}𝑑s`$
$`+{\displaystyle _1^t}C(1+ts)^mN(ts)s^k\gamma ^2(s)𝑑s`$
$`C_1t^mϵ(u_0)+C_2t^m\gamma ^2(t)`$
by Lemma C.1 and estimate (A.6). Therefore $`\gamma (t)C_1ϵ(u_0)+C_2\gamma ^2(t)`$ for all $`t1`$. If $`4C_1C_2ϵ(u_0)<1`$ and $`2C_0C_2ϵ(u_0)<1`$ (thus $`\gamma (1)<\left(2C_2\right)^1`$), then $`\gamma (t)C_0ϵ(u_0)`$ for all $`t`$ by an elementary algebra argument. This is the desired estimate for $`m=k/2`$, and it holds if $`ϵ(u_0)`$ is sufficiently small. If $`0<m<k`$ is arbitrary, then one can use the estimates for $`w`$ on $`(0,1]`$ and for $`[S(t)\mathrm{\Psi }]_x`$ to obtain for $`t1`$

$`w(,t)_{\mathrm{},m}`$ $`Ct^{mk}ϵ(u_0)+{\displaystyle _0^1}C(1+ts)^{mk}N(ts)w(,s)_{\mathrm{},m}𝑑s`$
$`+{\displaystyle _1^t}C(1+ts)^{mk}N(ts)F(v)(,s)_{\mathrm{},k}𝑑s`$
$`Ct^{mk}ϵ(u_0)+{\displaystyle _1^t}C(1+ts)^{mk}N(ts)w,s)_{\mathrm{},k/2}^2ds`$
$`Ct^{mk}ϵ(u_0)+Cϵ(u_0)^2{\displaystyle _1^t}(1+ts)^{mk}N(ts)(1+s)^k𝑑s`$
$`Ct^{mk}ϵ(u_0)`$
by the estimate for $`w,s)_{\mathrm{},k/2}`$ that was just established and by Lemma C.1. The theorem is now completely proved.

## 4. Generalized Korteweg-DeVries - Burgers Equations

Let us now look at the partial differential equation

$$u_t\alpha u_{xx}+u_{xxx}+g(u)_x=0(x,t>0).$$
$`4.1`$
The parameter $`\alpha `$ is positive, and $`g`$ is $`C^2`$-smooth. The case $`g(u)=(p+1)^1u^{p+1}`$ with integer $`p>0`$ is a model for long wave propagation in media with dissipation and dispersion. The special case $`p=1`$ is known as Korteweg-DeVries - Burgers  equation. It reduces to the Korteweg-DeVries equation if $`\alpha =0`$. Under certain conditions, the equation admits monotone traveling wave solutions $`u(x,t)=\varphi (xct)`$ with speed $`c`$ that connect the end states $`\varphi _\pm =lim_{r\pm \mathrm{}}\varphi (r)`$. Such a wave profile must satisfy the third order ordinary differential equation

$$c\varphi ^{}+g(\varphi )^{}+\varphi ^{\prime \prime \prime }\alpha \varphi ^{\prime \prime }=0.$$
$`4.2`$
An example is $`g(r)=2r(r1)(br)`$ with $`b2`$, which has the wave profile $`\varphi (r)=\left(1+e^r\right)^1`$ for the parameter $`\alpha =2b1`$ and the speed $`c=0`$. General profiles (not necessarily monotone) have been constructed in  and . It is known that monotone profiles exist for $`g(u)=(p+1)^1u^{p+1}`$ and $`\alpha 2\sqrt{pc}`$. A slightly more situation is the setting for the next result.

###### Proposition 4.1

Let $`gC^2`$ be strictly convex. A monotone wave profile $`\varphi `$ for (4.1) exists if and only if

$`c`$ $`={\displaystyle \frac{g(\varphi _+)g(\varphi _{})}{\varphi _+\varphi _{}}}`$ $`4.3a`$$`4.3b`$$`4.3c`$
$`\alpha `$ $`2\sqrt{g^{}(\varphi _{})c}`$
$`\varphi _+`$ $`<\varphi _{}.`$
The profile $`\varphi `$ must therefore be monotonically decreasing.

###### Demonstration Proof

Suppose $`\varphi `$ is a monotone wave profile with limits $`\varphi _\pm `$ at $`r=\pm \mathrm{}`$. Clearly $`c\varphi +g(\varphi )+\varphi ^{\prime \prime }\alpha \varphi ^{}=const.`$ and thus $`c\varphi _{}+g(\varphi _{})=c\varphi _++g(\varphi _+)`$, implying (4.3a). Set $`\psi (z)=\varphi _{}\varphi (z)`$ and

$$f(r)=g(\varphi _{})crg(\varphi _{}r).$$
$`4.4`$
Then $`f`$ is concave, and $`\alpha \psi ^{}\psi ^{\prime \prime }=f(\psi )`$. This is the equation for a wave profile $`\psi `$ of the Fisher - Kolmogorov-Petrovskii-Piskunov (F-KPP) equation $`v_tv_{xx}=f(v)`$ that travels to the right with speed $`\alpha `$ and has limits $`\psi _{}=\varphi _{}\varphi _+,\psi _+=0`$. It is known () that such a monotone wave profile for concave $`f`$ exists if and only if $`\alpha 2\sqrt{f^{}(0)}=2\sqrt{g^{}(\varphi _{})c}`$. In this case, $`\psi _{}>\psi _+`$ and therefore $`\varphi _{}>\varphi _+`$, since $`f`$ is positive between $`\psi _{}`$ and $`\psi _+`$. Thus (4.3b) and (4.3c) are true. Conversely, let (4.3a-c) hold. Define $`f`$ as in (4.4). By well-known results about the F-KPP equation, there exists a unique decreasing wave profile $`\psi `$ with $`\psi (0)=(\varphi _{}\varphi _+)/2`$ that moves to the right with speed $`\alpha `$. Then $`\varphi (z)=\varphi _{}\psi (z)`$ is a monotone wave profile for (4.1) with $`\varphi (\pm \mathrm{})=\varphi _\pm `$.

It is easy to see that in fact $`\varphi ^{}<0`$ on $``$. After rescaling the independent variable, adding a linear function to $`g`$, and changing to a moving coordinate system, one can assume that $`c=0,\varphi _{}=1`$, and $`\varphi _+=0`$. As in section 3, the focus is on the convergence of solutions of (4.1) to some translate $`\varphi (h)`$ as $`t\mathrm{}`$. As before, we can assume that $`h=0`$ and define $`\mathrm{\Psi }(x)=_{\mathrm{}}^x\left(u(z,0)\varphi (z)\right)𝑑z`$. In , (4.1) was discussed in the case $`g(x)=x^2`$, and it was shown that the difference $`u(,t)\varphi `$ converges to 0 in $`L^2`$ together with its derivatives if this difference is small in $`L^{2,k}`$ for some $`k>1`$ at $`t=0`$ and if sufficiently many derivatives of its derivatives are small in $`L^2`$. One of the results in  says that this convergence in fact holds if the initial difference is small just in $`L^2`$. The main result in  also covers the case $`g(u)=u^2`$ and states essentially that if $`\mathrm{\Psi }()_{2,k}`$ is sufficiently small, then $`u(,t)\varphi _{2,m}=O(t^{mk+ϵ})`$. Here $`ϵ=0`$ if $`2m2k`$ is an integer, and it is positive but arbitrarily small otherwise. Derivatives of $`u(,t)\varphi `$ are shown to decay at higher rates.

The main result of this section assumes that

$$g^{}(0)<0<g^{}(1)\text{and}g^{\prime \prime }(r)>0\text{for all}x.$$
$`4.5`$
###### Theorem 4.2

Let $`k>1`$ be a real number. There exists a constant $`C_0`$ such that if

$$\mathrm{\Psi }_{2,k}=ϵ(u_0)$$
$`4.6a`$
is sufficiently small, then the solution $`u`$ exists for all $`t>0`$, and for all $`0<m<k`$ and $`t1`$

$`u(,t)\varphi _{4,m}`$ $`C(1+t)^{mk}ϵ(u_0)`$ $`4.7a`$$`4.7b`$
$`u(,t)\varphi _{2,m}`$ $`C(1+t)^{mk}ϵ(u_0)`$
###### Demonstration Proof

The proof follows the scheme used in section 3. Two different function space settings ($`L^4`$ and $`L^2`$) are used to handle the quadratic nonlinearity. Estimates (B.7a, b) connect these settings.

Step 0. Let $`u`$ be a solution of (4.1), and set $`v(x,t)=_{\mathrm{}}^x\left(u(z,t)\varphi (z)\right)𝑑z`$. Thus $`v(x,0)=\mathrm{\Psi }(x)`$, and $`v`$ satisfies

$$v_t(,t)+Lv(,t)=G(v)(,t)$$

where now formally

$`Ly(x)`$ $`=\alpha y_{xx}(x)+y_{xxx}(x,t)+g^{}(\varphi (x))y_x(x)`$ $`4.8`$$`4.9`$
$`G(y)(x)`$ $`=g(\varphi (x))+g^{}(\varphi (x))y_x(x)g(\varphi (x)+y_x(x)=g^{\prime \prime }(\zeta (x))y_x^2(x).`$
Step 1. The operator $`L`$ acting in $`L^2()`$ is defined in (4.8), and $`L`$ generates a $`C_0`$ semigroup $`S(t)_{t0}`$ in this space by Appendix B. Let us note that $`\frac{d}{dx}g^{}(\varphi (x))=g^{\prime \prime }(\varphi (x))\varphi ^{}(x)<0`$, and thus $`S`$ is a contraction semigroup by (B.10). The semigroup can be restricted to the weighted spaces $`L^{2,k}`$ and $`L^{2,\rho }`$ ($`\rho <\alpha /3`$). From the properties of $`g`$ and the definition of $`G`$ in (4.8), one immediately deduces the estimates

$`G(y_1)G(y_2)_{2,k}`$ $`Cy_{1,x}y_{2,x}_{2,k}`$ $`4.10a`$$`4.10b`$
$`G(y)_{2,k}`$ $`Cy_x_{4,k/2}^2`$
for any $`k0`$ and suitable $`y,y_1,y_2`$, with some universal constant $`C`$. Suppose now that $`\mathrm{\Psi }L^{2,k}`$ is given. A continuous $`L^2`$-valued solution $`v(,t)`$ of the integral equation

$$v(,t)=S(t)\mathrm{\Psi }+_0^tS(ts)G(v)(,s)𝑑s$$
$`4.11`$
is again found by a contraction argument, using (4.10a), and due to the global Lipschitz property in (4.10a), this solution exists for all times. The existence and uniqueness arguments hold in all spaces $`L^{2,m}`$, $`0<mk`$, and therefore $`v(,t)`$ belongs to all these spaces for all $`t`$. The solution satisfies (4.7) in the classical sense if it is sufficiently smooth. Although the linear part of (4.7) does not enjoy maximal regularity properties, more smoothness for the solution follows easily from smoothness of the data, if the equation is differentiated and the results in Appendix B are used. Set $`w(x,t)=u(x,t)\varphi (x)=v_x(x,t)`$ and write $`N_0(t)=\mathrm{max}(1,t^{1/2})`$ and $`N_1(t)=\mathrm{max}(1,t^{5/8})`$.

Step 2. As before, let us next characterize the behavior of $`w`$ for $`0<t1`$ in more detail. The goal is to show the estimates

$`w(,t)_{2,m}CN_0(t)ϵ(u_0)`$ $`4.12a`$$`4.12b`$
$`w(,t)_{4,m}CN_1(t)ϵ(u_0)`$
for $`0<t1`$. For this purpose consider the integral equation for $`w`$ which is obtained by differentiating (4.1), i.e. (3.13) with $`F`$ replaced by $`G`$. For the “free term” $`[S(t)\mathrm{\Psi }]_x`$ in (3.13), (4.12a, b) follows directly from (B.7a, b). To show (4.12a), one uses (B.7a) and (3.10b) to derive inequality (3.14) with $`w(,t)_{\mathrm{},m}/N(t)`$ replaced everywhere by $`w(,t)_{2,m}/N_0(t)`$. For (4.12b), one uses (B.7b) and arrive at an inequality like (3.14) for $`w(,t)_{4,m}/N_1(t)`$. Standard arguments for linear integral inequalities imply that $`w(,t)_{2,m}/N_0(t)+w(,t)_{4,m}/N_1(t)C`$ on $`(0,1]`$, i.e. (4.12a, b).

Step 3. Finally, let us prove estimate (4.6b). Start by observing that (B.14a,b) are just estimates (4.6a,b) for the “free” term $`[S(t)\mathrm{\Psi }]_x`$. Consider again first the special case $`m=k/2`$ and thus $`mk=m`$. Define the quantity

$$\gamma (t)=\underset{1st}{sup}s^mw(,s)_{4,m}$$
$`4.13`$
for $`t1`$. Then $`\gamma `$ is continuous, $`\gamma (1)=w(,1)_{4,m}C_0ϵ(u_0)`$ for some $`C_0`$, and $`G(v)(,s)_{2,m}Cs^k\gamma ^2(s)`$ for all $`s`$. Equation (4.11), estimate (4.12) and Lemma C.1 imply that for $`t1`$

$`w(,t)_{4,m}`$ $`Ct^mϵ(u_0)+{\displaystyle _0^1}C(1+ts)^mN_1(ts)w(,s)_{2,m}𝑑s`$
$`+{\displaystyle _1^t}C(1+ts)^mN_1(ts)s^k\gamma ^2(s)𝑑s`$
$`C_1t^mϵ(u_0)+C_2t^m\gamma ^2(t).`$
Therefore $`\gamma (t)C_1ϵ(u_0)+C_2\gamma ^2(t)`$ for all $`t1`$. As in section 3, $`\gamma (t)C_0ϵ(u_0)`$ follows if $`ϵ(u_0)`$ is sufficiently small, which is the desired estimate. If $`0<m<k`$ is arbitrary, then one obtains as in section 3 for $`t1`$

$`w(,t)_{4,m}`$ $`Ct^{mk}ϵ(u_0)+{\displaystyle _0^1}C(1+ts)^{mk}N_1(ts)w(,s)_{2,m}𝑑s`$
$`+{\displaystyle _1^t}C(1+ts)^{mk}N_1(ts)G(v)(,s)_{2,k}𝑑s`$
$`Ct^{mk}ϵ(u_0)+{\displaystyle _1^t}C(1+ts)^{mk}N_1(ts)w,s)_{4,k/2}^2ds`$
$`Ct^{mk}ϵ(u_0).`$
Finally for $`t1`$,

$`w(,t)_{2,m}`$ $`Ct^{mk}ϵ(u_0)+{\displaystyle _0^t}C(1+ts)^{mk}N_0(ts)w(,s)_{4,k/2}^2𝑑s`$
$`Ct^{mk}ϵ(u_0)+{\displaystyle _0^t}C(1+ts)^{mk}N_0(ts)(1+s)^kϵ^2(u_0)𝑑s`$
$`Ct^{mk}ϵ(u_0).`$
The theorem is now completely proved.

Let us note in concluding that estimates (B.7a,b) can easily be modified to

$$S(t)\varphi _{\mathrm{},k}Ct^{1/4}\varphi _{2,k},[S(t)\varphi ]_x_{\mathrm{},k}Ct^{3/4}\varphi _{2,k}$$

for $`0<t1`$. The last argument in the proof of Theorem 4.1 then implies that also

$$u(,t)\varphi _{\mathrm{},m}Ct^{mk}ϵ(u_0).$$

## Appendix A: Linearization of Scalar Viscous Conservation Laws

In this section, properties of solutions of the equation

$$u_tu_{xx}+cu_x+du=0(x,t>0)$$
$`A.1`$
are collected that are used in the main part of this paper. Here $`c`$ and $`d`$ are suitable coefficient functions which are bounded together with their first derivatives. The results are mostly well known.

Define the Banach spaces $`X=L^{\mathrm{}}()C^0()`$ and $`X_0=\{yX|y\text{is uniformly continuous}\}`$, both equipped with the supremum norm $`_{\mathrm{}}`$. Define the operator $`L`$ acting in $`X`$ by $`L\phi =\phi _{xx}+c\phi _x+d\phi `$ for $`\phi D(L)=\{yXy^{\prime \prime }X\}`$. The restriction of $`L`$ to $`X_0`$ will also be denoted by $`L`$. Then $`L`$ generates a $`C_0`$ semigroup $`(S(t))_{t0}`$ in $`X_0`$ which can be constructed as a perturbation of the heat semigroup. In addition, the estimates hold

$`S(t)\phi _{\mathrm{}}`$ $`\phi _{\mathrm{}}(0t<\mathrm{})\text{if}d=0`$ $`A.2`$$`A.3`$
$`\left[S(t)\phi \right]_x_{\mathrm{}}`$ $`{\displaystyle \frac{C}{\sqrt{t}}}\phi _{\mathrm{}}(0<t1)`$
The first estimate is the maximum principle, the second follows from estimates for fundamental solutions in . The semigroup can be extended to act on $`X`$, and the extension will also be denoted by $`S(t)`$. The extension still satisfies the estimates above. It is only strongly continuous for $`t>0`$, but this is irrelevant for the purposes of this paper.

Let us next examine the behavior of this semigroup in subspaces of $`X_0`$ that are defined by means of weight functions. Let $`w:[1,\mathrm{})`$ be smooth, with $`w(\pm \mathrm{})=\mathrm{}`$, and with the first three derivatives of $`x\mathrm{log}w(x)`$ bounded. Set $`X_w=\{yX|ywX\}`$ with the norm $`y_{\mathrm{},w}=wy_{\mathrm{}}`$. There is a natural bijection $`R:X_wX,Ry=wy`$, inducing the (formal) conjugate $`\stackrel{~}{S}(t)=RS(t)R^1`$ of $`S(t)`$ which again acts on $`X`$. A straight forward computation shows that the restriction of $`\stackrel{~}{S}(t)`$ to $`X_0`$ has an infinitesimal generator $`\stackrel{~}{L}`$ of the same form as $`L`$. The estimates in  apply to this more general case and imply that $`\stackrel{~}{S}(t)`$ acts on $`X_0`$ as a $`C_0`$ \- semigroup. Thus $`S(t)`$ acts on each $`X_w`$ with estimates $`S(t)y_{\mathrm{},w}Ce^{Mt}y_{\mathrm{},w}`$ for $`0t<\mathrm{}`$ and $`\left[S(t)y\right]_x_{\mathrm{},w}{\displaystyle \frac{C}{\sqrt{t}}}y_{\mathrm{},w}`$ for $`0<t1`$. The constants $`C`$ and $`M`$ now depend also on $`w`$. The same notation for $`S(t)`$ will be used, whether it acts on $`X`$ or on $`X_w`$.

Now let $`\rho ,k>0`$ and consider specifically the spaces $`X_\rho `$ and $`X_k`$ that correspond to the weight functions $`w(x)=\mathrm{cosh}(\rho x)`$ and $`w(x)=\left(1+x^2\right)^{k/2}`$ (as in section 1). The work in  and  implies that for sufficiently small $`\rho `$ there exists $`ϵ=ϵ(\rho )>0`$ and $`C>0`$ such that for all $`\phi X_\rho `$ and all $`t0`$

$$S(t)\phi _{\mathrm{},\rho }Ce^{ϵt}\phi _{\mathrm{},\rho }.$$
$`A.4`$
Theorem 2.2 and (A.2) then imply that for all $`\phi X_k`$ and all $`0<m<k`$

$$S(t)\phi _{\mathrm{},m}C_1(1+t)^{mk}\phi _{\mathrm{},k}.$$
$`A.5`$
Using (A.3), it also follows that for all $`t>0`$

$$\left[S(t)\phi \right]_x_{\mathrm{},m}N(t)C_1(1+t)^{mk}\phi _{\mathrm{},k}$$
$`A.6`$
with $`N(t)=\mathrm{max}\{1,t^{1/2}\}`$. Finally, since $`[S(t)\phi ]_x`$ satisfies a parabolic equation of the same form as (A.1), the estimates hold for $`0<t1`$

$$[S(t)\phi ]_x_{\mathrm{},k}CN(t)[\phi ]_x_{\mathrm{},k}\text{and}[S(t)\phi ]_x_{\mathrm{}}C[\phi ]_x_{\mathrm{}}.$$
$`A.7`$
## Appendix B: Linearization of Generalized Korteweg-De Vries - Burgers Equations

This appendix collects properties of solutions of the equation

$$u_t+u_{xxx}\alpha u_{xx}+cu_x+du=0(x,t>0)$$
$`B.1`$
which are used in the main part of the paper. Here $`\alpha >0`$, and $`c`$ and $`d`$ are suitable smooth coefficient functions. Proofs will only be indicated.

If $`c=d=0`$, the spatial Fourier transform $`\widehat{u}(,t)`$ of the solution $`u`$ of (B.1) is given by

$$\widehat{u}(\xi ,t)=e^{(i\xi ^3\alpha \xi ^2)t}\widehat{u}(\xi ,0)$$
$`B.2`$
This defines a contraction semigroup $`S_0(t)_{t0}`$ in $`L^2`$. Moreover, for $`\rho <\alpha /3`$, $`S_0`$ maps $`L^{2,\rho }`$ into itself, since Fourier transforms of functions in $`L^{2,\rho }`$ have analytic extensions into the strip $`\{z||\mathrm{}(z)|<\rho \}`$ that are square integrable on its boundary and since the multiplier in (B.2) is bounded by $`Ce^{(\alpha 3\rho )|\mathrm{}(\xi )|^2t}`$ on any such strip. Thus $`S_0`$, restricted to any such $`L^{2,\rho }`$, also generates a $`C_0`$ semigroup there. In addition, since the multiplier $`\xi ^2te^{(i\xi ^3\alpha \xi ^2)t}`$ is similarly bounded, there are the estimates for all $`t`$

$$[S_0(t)\phi ]_{xx}_2\frac{C}{t}\phi _2\text{and}[S_0(t)\phi ]_{xx}_{2,\rho }\frac{C}{t}\phi _{2,\rho }$$
$`B.3`$
and from standard calculus estimates

$$[S_0(t)\phi ]_x_2\frac{C}{\sqrt{t}}\phi _2\text{and}[S_0(t)\phi ]_x_{2,\rho }\frac{C}{\sqrt{t}}\phi _{2,\rho }.$$
$`B.4`$
By Theorem 2.1, $`S_0`$ therefore maps each $`L^{2,k}`$ into itself and satisfies similar estimates for $`0<k<\mathrm{}`$.

By a standard perturbation argument, one obtains existence and uniqueness of solutions of (B.1)for general coefficients $`c,d`$ that are bounded together with their first and second derivatives. The resulting semigroup $`S(t)_{t0}`$ maps $`L^2,L^{2,\rho }(\rho <1/3),L^{2,k}(0<k<\mathrm{})`$ into itself, and the estimates hold

$`S(t)\phi _{2,k}`$ $`C\phi _{2,k}`$ $`B.5a`$$`B.5b`$$`B.5c`$
$`[S(t)\phi ]_x_{2,k}`$ $`{\displaystyle \frac{C}{\sqrt{t}}}\phi _{2,k}`$
$`[S(t)\phi ]_{xx}_{2,k}`$ $`{\displaystyle \frac{C}{t}}\phi _{2,k}.`$
Next note that for differentiable initial data $`u(,0)`$, the derivative $`u_x`$ also satisfies an equation of the form (B.1). Thus there is also the estimate

$$[S(t)\phi ]_x_2C\phi _x_2(0<t1)$$
$`B.6`$
Finally, the estimates

$`S(t)\phi _{4,k}Ct^{1/8}\phi _{2,k}(0<t1)`$ $`B.7a`$$`B.7b`$
$`[S(t)\phi ]_x_{4,k}Ct^{5/8}\phi _{2,k}(0<t1)`$
for all $`\phi L^{2,k}`$ follows from (B.5a,b,c) and the calculus inequality $`v_{4,k}^4Cv_{2,k}^3\left(v_{2,k}+v_x_{2,k}\right)`$.

Let us next turn to estimates for $`S(t)`$ for large $`t`$, in the special case where $`d=0`$. Consider a general weight function $`w>0`$ with the properties

$$|w^{}|\kappa w,|w^{\prime \prime }|+|w^{\prime \prime \prime }|C,\text{with}\kappa <\alpha /3,C>0.$$

A calculation shows that for $`u(,t)=S(t)\phi `$ with $`\phi L^{2,\rho },\rho <\alpha /3`$, and $`t>0`$

$$0=\frac{d}{dt}w^2u+u_x^2\left(3|w^{}|^2+2\alpha w^2\right)u^2\left(\alpha (w^2)^{\prime \prime }+(w^2)^{\prime \prime \prime }+(cw^2)^{}\right).$$
$`B.8`$
Consider first the case $`w=1`$. Then this identity implies

$$\frac{d}{dt}u^2c^{}u^2.$$
$`B.9`$
Therefore if $`c^{}0`$, then

$$S(t)\phi _2\phi _2$$
$`B.10`$
as was observed in . Next fix the assumptions

$$c^{}(x)<0(x),\underset{x\mathrm{}}{lim}c(x)=c_L>0>c_R=\underset{x+\mathrm{}}{lim}c(x).$$
$`B.11`$
Under these assumptions, for all sufficiently small $`\rho >0`$ there are $`C,\gamma >0`$ such that

$$S(t)\phi _{2,\rho }Ce^{\gamma t}\phi _{2,\rho }.$$
$`B.12`$
Indeed, consider the weights $`w(x)=\sqrt{\mathrm{cosh}(\rho (xx_0))}`$ where $`x_0`$ is such that $`c(x_0)=0`$. From (B.8) one then obtains

$$\frac{d}{dt}w^2u^2F_\rho w^2u^2$$
$`B.13`$
where

$$F_\rho (x)=\alpha \rho ^2+c^{}(x)+\left(\rho c(x)+\rho ^3\right)\mathrm{tanh}\rho (xx_0).$$

Let us show that for all sufficiently small $`\rho >0`$, $`F_\rho <\gamma <0`$ on $``$. Set $`2A=\mathrm{min}(c_L,c_R)`$. If $`M`$ is sufficiently large (depending on $`A`$) and $`\rho ^2<A`$, then for $`|xx_0|>M`$

$$F_\rho (x)\alpha \rho ^2+(\rho ^2A)\rho \mathrm{tanh}\rho M\gamma _0\rho ^2$$

for some $`\gamma _0>0`$. If now $`|xx_0|M`$, then

$$F_\rho (x)\alpha \rho ^2+C\rho +\underset{|xx_0|M}{sup}c^{}(x)\gamma _1$$

provided $`\rho `$ is decreased further. Then (B.13) implies $`\frac{d}{dt}w^2u^2\gamma w^2u^2`$, and (B.12) follows. Theorem 2.1 and (B.5) now imply that for all $`\phi L^{2,k}`$ and $`0<m<k,t>0`$

$`S(t)\phi _{2,m}`$ $`C(1+t)^{mk}\phi _{2,k}`$ $`B.14a`$$`B.14b`$
$`[S(t)\phi ]_x_{2,m}`$ $`C(1+t)^{mk}N_0(t)\phi _{2,k}`$
with $`N_0(t)=\mathrm{max}(1,t^{1/2})`$. Finally note that (B.7) and (B.13) imply

$`S(t)\varphi _{4,m}`$ $`Ct^{mk}\varphi _{2,k}`$ $`B.15a`$$`B.15b`$
$`[S(t)\varphi ]_x_{4,m}`$ $`Ct^{mk}\varphi _{2,k}`$
for $`t1`$, for all $`0<m<k`$.

## Appendix C: An Integral Estimate

###### Lemma C.1

Let $`0<\alpha <\beta `$ with $`\beta >1`$. Let $`M:(0,\mathrm{})`$ be bounded on $`[1,\mathrm{})`$ and integrable on $`(0,1)`$. Then there exists a constant $`C=C_{\alpha ,\beta }`$ such that for all $`t0`$

$$_0^tM(ts)(1+ts)^\alpha (1+s)^\beta 𝑑sCt^\alpha .$$

###### Demonstration Proof

Split the integral:

$$_0^tM(ts)(1+ts)^\alpha (1+s)^\beta 𝑑s_0^t(1+ts)^\alpha (1+s)^\beta 𝑑s+_0^1M(s)(1+s)^\alpha (1+ts)^\beta 𝑑s.$$

The second integral is clearly bounded by $`C(1+t)^\beta `$. The first integral becomes

$$_0^t(1+ts)^\alpha (1+s)^\beta 𝑑s=_0^{t/2}(1+ts)^\alpha (1+s)^\beta 𝑑s+_{t/2}^t(1+ts)^\alpha (1+s)^\beta 𝑑s.$$

If $`\alpha ,\beta >1`$, one estimates further

$$\mathrm{}(\beta 1)^1(1+t/2)^\alpha +(\alpha 1)(1+t/2)^\beta C(1+t)^\alpha .$$

If $`\alpha <1<\beta `$, the same estimate yields

$$\mathrm{}(\beta 1)^1(1+t/2)^\alpha +(1\alpha )(1+t)^{1\alpha }(1+t/2)^\beta C(1+t)^\alpha .$$

Finally if $`\alpha =1<\beta `$, then the modified estimate holds

$$\mathrm{}(\beta 1)^1(1+t/2)^1+\mathrm{log}(1+t/2)(1+t/2)^\beta C(1+t)^1.$$

This proves the lemma. A closer look at the proof shows that the exponent $`\alpha `$ on the right hand side cannot be improved.


text:
# Vortices in a rotating Bose-Einstein condensate: critical velocities and energy diagrams in the Thomas-Fermi regime

## I Introduction

Since the first experimental achievement of Bose-Einstein condensates in atomic gases in 1995, many properties of these systems have been studied experimentally and theoretically, and particularly the existence of vortices BuR ; CD ; DGPS ; FCS ; F ; FS ; MCWD ; SF . A way to create vortices consists in rotating the trap confining the atoms: for sufficiently large velocities, it becomes energetically favorable to have vortices in the system. Theoretical studies of this type of experiments have often been made in the framework of the nonlinear Schrödinger equation or Gross-Pitaevskii equation, well known for superfluids, but which provides a very good description of Bose-Einstein condensates: it is assumed that the $`N`$ particles of the gas are condensed in the same state for which the wave function $`\varphi `$ minimizes the Gross-Pitaevskii energy. By introducing a rotating frame at the angular velocity $`\stackrel{~}{𝛀}=\stackrel{~}{\mathrm{\Omega }}𝐞_z`$, the trapping potential becomes time independent, and the wave function $`\varphi `$ minimizes the energy

$`_{3D}(\varphi )`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}|\varphi |^2}+{\displaystyle \frac{m}{2}}{\displaystyle \underset{\alpha }{}}\omega _\alpha ^2r_\alpha ^2|\varphi |^2`$ (1)
$`+{\displaystyle \frac{N}{2}}g_{3D}|\varphi |^4\mathrm{}\stackrel{~}{𝛀}(i\varphi ,\varphi \times 𝐱),`$
under the constraint $`|\varphi |^2=1`$. Here, for any complex quantities $`u`$ and $`v`$ and their complex conjugates $`\overline{u}`$ and $`\overline{v}`$, $`(u,v)=(u\overline{v}+\overline{u}v)/2`$. The terms in the energy correspond to the kinetic energy, the trapping potential energy, the interaction energy and the inertial due to the change of frame.

When the atoms are strongly confined in the $`z`$axis, the situation can be simplified to a two dimensional problem where the wave function $`\psi `$ depends only on $`𝐱=(x,y)`$ and it minimizes

$`_{2D}(\psi )`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}|\psi |^2}+{\displaystyle \frac{m}{2}}{\displaystyle \underset{\alpha =x,y}{}}\omega _\alpha ^2r_\alpha ^2|\psi |^2`$ (2)
$`+{\displaystyle \frac{N}{2}}g|\psi |^4\mathrm{}\stackrel{~}{𝛀}(i\psi ,\psi \times 𝐱),`$
where $`g=g_{3D}(m\omega _z/2\pi \mathrm{})`$. The constraint $`|\psi |^2=1`$ is also imposed. Our study was at first motivated by the work of Castin and Dum CD  who have studied the equilibrium configurations by looking for the minimizers in a reduced class of functions for the 2D case and have done numerical computations in 2 and 3D. Their analysis is in the Thomas-Fermi regime, where the mean interaction energy per particle is larger than $`\mathrm{}\omega _{x,y}`$.

Our aim is to provide a mathematical framework for a rigorous study of the energy $`_{2D}`$ and its minimizers in the Thomas Fermi limit. We first observe that this energy has a striking similarity with the high-kappa, high-field limit of the Ginzburg-Landau free energy used in the modeling of superconductors. Thus, we expect the energy $`_{2D}`$ will develop similar behavior as those for the Ginzburg-Landau energy studied in ASS ; S1 ; S . In particular, the results obtained in the context of Ginzburg-Landau energy may be applied to $`_{2D}`$, in the Thomas Fermi regime, to yield an asymptotic development of the energy as well as the critical velocities for the nucleation of vortices and the location of these vortices. Due to the close resemblance, we will not be concerned with the detailed derivations in this paper, but rather focus on the conclusions one can draw from the asymptotic developments. To our knowledge, some of our estimates to be presented later have not been given in the literature previously. Let us point out that our method to compute the energy is also new in this context and very different from the ones in, for instance, CD ; FCS ; SF .

We define the characteristic length $`d=(\mathrm{}/m\omega _x)^{1/2}`$ and set

$$\epsilon ^2=\frac{\mathrm{}^2}{2Ngm}.$$

In the Thomas-Fermi approximation, $`\epsilon `$ is small, which will be our asymptotic regime. We re-scale the distance by $`R=d/\sqrt{\epsilon }`$ and define $`u(𝐫)=R\psi (𝐱)`$ where $`𝐱=R𝐫`$. Assume that $`\omega =\omega _x`$ and $`\omega _y=\lambda \omega `$ with $`0\lambda 1`$ and we set $`\mathrm{\Omega }=\stackrel{~}{\mathrm{\Omega }}/\epsilon \omega `$. Since the trapping potential is stronger than the inertial potential, we have $`\mathrm{\Omega }<1/\epsilon `$. The energy can be rewritten as:

$`E_{2D}(u)`$ $`=`$ $`{\displaystyle \frac{1}{2}|u|^2}+{\displaystyle \frac{1}{2\epsilon ^2}}(x^2+\lambda ^2y^2)|u|^2+{\displaystyle \frac{1}{4\epsilon ^2}}|u|^4`$ (3)
$`+𝛀(iu,u\times 𝐫).`$
Due to the constraint $`|u|^2=1`$, we can add to $`E_{2D}`$ any multiple of $`|u|^2`$ so that it is equivalent to minimize

$$|u|^2+2𝛀(iu,u\times 𝐫)+\frac{1}{2\epsilon ^2}|u|^4\frac{1}{\epsilon ^2}a(𝐫)|u|^2$$

where $`a(𝐫)=\alpha (x^2+\lambda ^2y^2)`$ for some constant $`\alpha `$ to be determined. Let $`𝒟`$ be the ellipse $`\{a>0\}=\{x^2+\lambda ^2y^2<\alpha \}`$. We impose the following constraint on $`\alpha `$:

$$_𝒟a(𝐫)=1.$$
(4)
Indeed, as $`\epsilon `$ tends to 0, the minimizer will satisfy that $`|u|^2`$ will be close to $`a`$ so that the constraint will be satisfied automatically by $`u`$ if we impose (4). Equation (4) leads to $`\alpha ^2=2\lambda /\pi `$. If $`\lambda =1`$, that is $`\omega _x=\omega _y`$, then $`𝒟`$ is a disc of radius $`R_0`$ with $`R_0^4=2/\pi `$.

To study the problem analytically, it is reasonable to minimize the energy over the domain $`𝒟`$ with zero boundary data for $`u`$. Indeed, when $`a0`$, the energy is convex so that the minimizer $`u`$ goes to zero exponentially at infinity (see the numerical observation in CD  and the analysis on the behavior near the boundary of $`𝒟`$ as well as the decay at infinity of the order parameter in DPS ; FF ). Denote $`H^m(𝒟)`$ the space of square integrable functions defined on the domain $`𝒟`$ that have square integrable derivatives up to the order $`m`$, and $`H_0^1(𝒟)`$ the space of functions in $`H^1(𝒟)`$ satisfying the zero boundary condition. Denote the norm $`\left(_𝒟|v|^2\right)^{1/2}`$ by $`v`$ for any square integrable function $`v`$, we then consider the problem

$$\mathrm{min}E_\epsilon (u)\text{subject to }uH_0^1(𝒟),u=1$$
$`(P)`$
where

$`E_\epsilon (u)`$ $`=`$ $`{\displaystyle _𝒟}|u|^2+2𝛀(iu,u\times 𝐫)`$ (5)
$`+{\displaystyle \frac{1}{2\epsilon ^2}}(a(𝐫)|u|^2)^2.`$
Note that critical point $`u`$ of $`E_\epsilon `$ is a solution of

$$\mathrm{\Delta }u2i(𝛀\times 𝐫).u=\frac{1}{\epsilon ^2}u(a|u|^2)+\mu _\epsilon u\text{ in }𝒟,$$
(6)
with $`u=0`$ on $`𝒟`$ and $`\mu _\epsilon `$ is the Lagrange multiplier. The specific choice of $`\alpha `$ in (4) will imply that the term $`\mu _\epsilon u`$ is negligible in front of $`au/\epsilon ^2`$.

We want to study the behavior of $`\mathrm{min}E_\epsilon (u)`$ as $`\epsilon `$ goes to 0. In Section 2, we compute an asymptotic development of the energy and in Section 3, the critical velocities of nucleation of vortices and the location of the vortices. In Section 4, we study the evolution in imaginary time and construct some numerical algorithms. In Section 5, we present some computational results and the energy diagrams.

## II Asymptotic development of the energy

To study the behavior of the minimizer of the energy when $`\epsilon `$ goes to zero, we observe that the form of the energy (5) is close to the Ginzburg-Landau functional studied in ASS , S1  where the magnetic field has been replaced by a rotating term and similar to S  except for the trapping potential and the minimization over constraint. The main idea is to decouple the energy into three terms: a part coming from the solution without vortices, a vortex contribution and a term due to rotation. The estimate on the vortex contribution was developed in BBH ; BR ; S1 ; S .

### II.1 The solution without vortices

Firstly, we are interested in solutions without vortices, that is $`u`$ has no zero in the interior of $`𝒟`$. Thus we consider functions of the form $`\eta =fe^{iS}`$, where $`\eta `$ is in $`H_0^1(𝒟)`$ and $`f`$ is real and has no zero in the interior of $`𝒟`$. We consider first minimizing $`E_\epsilon `$ over such functions without imposing the constraint that the $`L^2`$ norm is 1, that is, $`f`$ and $`S`$ minimize

$`_\epsilon (f,S)`$ $`=`$ $`{\displaystyle _𝒟}|f|^2+{\displaystyle \frac{1}{2\epsilon ^2}}(af^2)^2`$ (7)
$`+{\displaystyle f^2|S𝛀\times 𝐫|^2}f^2\mathrm{\Omega }^2r^2.`$
We have $`f=0`$ on $`𝒟`$ and

$$\mathrm{\Delta }f+fS(S2𝛀\times 𝐫)=\frac{1}{\epsilon ^2}f(af^2)\text{ in }𝒟,$$
(8)
$$\mathrm{div}(f^2(S𝛀\times 𝐫))=0.$$
(9)
Equation (9) implies that there exists $`\xi `$ in $`H^2(𝒟)H_0^1(𝒟)`$ such that

$$f^2(S𝛀\times 𝐫)=\mathrm{\Omega }^{}\xi ,$$
(10)
where $`^{}\xi =(_y\xi ,_x\xi )`$. So $`\xi `$ is the unique solution of

$$\mathrm{div}(\frac{1}{f^2}\xi )=0\text{ in }𝒟,\xi =0\text{ on }𝒟.$$
(11)
Note that equation (10) is the equation for the velocity, but we prefer to write it as an orthogonal gradient for our later purposes of integration by part. In the special case where $`𝒟`$ is a disc, the minimum of (7) is reached for $`S=0`$ but this is not the case if $`𝒟`$ is an ellipse and there is a non trivial solution of (9).

#### The case of the disc

Assume that $`\lambda =1`$ so that $`𝒟`$ is a disc and $`a(r)=R_0^2r^2`$. As discussed earlier, $`S=0`$ in this case so that the energy becomes

$$_\epsilon (f)=_𝒟|f|^2+\frac{1}{2\epsilon ^2}(a(r)f^2)^2.$$
(12)
Let $`\eta _\epsilon `$ be the minimizer of (12) over real valued functions in $`H_0^1(𝒟)`$. Then, $`\eta _\epsilon `$ has no vortex, is independent of $`\mathrm{\Omega }`$ and satisfies

$$\mathrm{\Delta }\eta _\epsilon =\frac{1}{\epsilon ^2}\eta _\epsilon (\eta _\epsilon ^2a)\text{in}𝒟,\eta _\epsilon =0\text{on}𝒟.$$

When $`\epsilon `$ tends to 0, $`\eta _\epsilon ^2`$ is close to $`a`$ except on a boundary layer of size $`\epsilon ^{2/3}`$ close to $`𝒟`$. More precisely, using sub and super solutions, one can verify that

$$\sqrt{a(r)}\mathrm{tanh}(\delta (\sqrt{a(r)})^3|\mathrm{log}\epsilon |)\eta _\epsilon \sqrt{a(r)}$$

for $`|R_0^2r^2|C\epsilon ^{1/3}`$. In fact, one can construct a sub solution of the type above in any region $`|R_0^2r^2|C\epsilon ^{\beta /3}`$ with $`\beta <2`$. Then the value of $`\delta `$ is less than $`c(2\beta )`$.

The boundary layer can be analyzed using the change of variable $`x=(R_0r)/\epsilon ^{2/3}`$ and $`v_\epsilon (x)=\eta _\epsilon (r)/\epsilon ^{2/3}`$. Then $`v_\epsilon `$ satisfies the Painlevé equation

$$v^{\prime \prime }=v(v^22R_0x),v(0)=0,v(x)2R_0x\text{for}x\text{large}.$$

The boundary behavior has already been studied in FF  and DPS , but using matched asymptotics.

The energy of $`\eta _\epsilon `$ can also be estimated by a test function equal to $`\sqrt{a}`$ except on the boundary layer to get

$$E_\epsilon (\eta _\epsilon )\frac{2\pi }{3}|\mathrm{log}\epsilon |(1+o(1)).$$
(13)
#### The case of an ellipse

As discussed before, the minimum $`\eta _\epsilon =f_\epsilon e^{iS_\epsilon }`$ of (7) has a non-trivial phase. $`f_\epsilon ^2`$ tends to $`a`$ in every compact subset of $`𝒟`$ and the function $`\xi _\epsilon `$ given by (10) or (11) tends to the unique solution $`\xi `$ of

$$\mathrm{div}(\frac{1}{a}\xi )=0\text{ in }𝒟,\xi =0\text{ on }𝒟.$$
(14)
One can easily get that $`\xi (x,y)=a^2(x,y)/(2+2\lambda ^2)`$. Using (10), we can define $`S_0`$, the limit of $`S_\epsilon `$, to be the solution of $`a(S_0𝛀\times 𝐫)=\mathrm{\Omega }^{}\xi `$ with zero value at the origin. We have $`S_0=C\mathrm{\Omega }xy`$ with $`C=(\lambda ^21)/(\lambda ^2+1)`$. We see that $`S_0`$ cancels when $`\lambda =1`$ that is in the case of the disc. This computation is consistent with the one in SF , though it is derived in a different way.

### II.2 Decoupling the energy

Let $`\eta _\epsilon =f_\epsilon e^{iS_\epsilon }`$ be the vortex free minimizer of $`E_\epsilon `$ discussed previously without imposing the constraint on the norm of $`u`$. Let $`u_\epsilon `$ be a minimizer of $`E_\epsilon `$ under the constraint $`_𝒟|u|^2=1`$ and let $`v_\epsilon =u_\epsilon /\eta _\epsilon `$. Since $`\eta _\epsilon `$ satisfies the Gross Pitaevskii equation (8)-(9), we have

$`{\displaystyle _𝒟}(|v|^21)({\displaystyle \frac{1}{2}}\mathrm{\Delta }f_\epsilon ^2{\displaystyle \frac{1}{\epsilon ^2}}f_\epsilon ^2(af_\epsilon ^2)+|f_\epsilon e^{iS_\epsilon }|^2`$
$`2f_\epsilon ^2(S_\epsilon 𝛀\times 𝐫)=0.`$
Using this identity, one can get that the energy $`E_\epsilon (u_\epsilon )`$ decouples as follows

$`E_\epsilon (u_\epsilon )=E_\epsilon (\eta _\epsilon )+G_\epsilon (v_\epsilon )`$
$`+2{\displaystyle _𝒟}|\eta _\epsilon |^2(S_\epsilon 𝛀\times 𝐫)(iv_\epsilon ,v_\epsilon ),`$ (15)
where

$$G_\epsilon (v_\epsilon )=_𝒟|\eta _\epsilon |^2|v_\epsilon |^2+\frac{|\eta _\epsilon |^4}{2\epsilon ^2}(1|v_\epsilon |^2)^2.$$

This decoupling was used in S  in the case of a disc where $`S_\epsilon =0`$.

### II.3 Estimate on the energy

We now estimate the terms in (II.2). The first term $`E_\epsilon (\eta _\epsilon )`$ is a constant only depending on $`\epsilon `$, and not on the solution type, that is with or without vortices. The second term gives a contribution coming from the vortices and the third term is due to the vortices and rotation.

We use the analysis of vortices developed in BBH  and later in BR ; S1 ; S . Let $`𝒟_\epsilon ^{}=D\{x,\mathrm{dist}(x,𝒟)\epsilon ^\beta \}`$, with $`\beta <1`$. Then in $`𝒟_\epsilon ^{}`$ it is possible to define vortices for $`v_\epsilon `$ in the following way: there exist balls $`B_i=B(p_i,\epsilon ^\beta ^{})`$ where $`p_i`$ are points in $`𝒟_\epsilon ^{}`$ at mutual distance bigger than $`8\epsilon ^\beta ^{}`$ and $`\beta ^{}>\beta `$, such that $`|v_\epsilon |1/2`$ in $`𝒟_\epsilon ^{}_iB_i`$. Moreover, the degree $`d_i=deg(v_\epsilon /|v_\epsilon |,B_i)`$ is not zero and there is an estimate of the energy of $`v_\epsilon `$ in each ball $`B_i`$. This analysis means that vortices are in fact defined in the balls where $`v_\epsilon `$ is less than $`1/2`$ and has a non zero degree. This allows us to compute $`G_\epsilon (v_\epsilon )`$ for which only the gradient term in the vortex balls will give a contribution: each vortex gives a contribution in the amount of $`2\pi |\mathrm{log}\epsilon |`$ due to its degree and a contribution of lower order which comes from the interaction with the other vortices. Moreover, $`|\eta _\epsilon |^2`$ is almost $`a`$:

$`G_\epsilon (v)=2\pi |\mathrm{log}\epsilon |{\displaystyle \underset{i}{}}|d_i|a(p_i)`$
$`2\pi {\displaystyle \underset{ij}{}}d_id_ja(p_i)a(p_j)\mathrm{log}|p_ip_j|+O(1).`$ (16)
In order to estimate the third term in (II.2) we let $`X_\epsilon `$ be the solution of $`a(𝐫)(S_0𝛀\times 𝐫)=\mathrm{\Omega }^{}X_\epsilon `$ which is zero on the boundary of $`𝒟_\epsilon ^{}`$. That is $`X_\epsilon `$ solves (14) but with zero boundary data on $`𝒟_\epsilon ^{}`$ instead of $`𝒟`$. Hence $`X_\epsilon `$ converges to $`\xi `$ the solution of (14). An integration by part on the last term of (II.2) using the definition of $`X_\epsilon `$ and the definition of the degree of $`v_\epsilon `$ on vortex balls and the fact that the higher order term comes from an integration on the vortex balls yields

$`{\displaystyle _𝒟}|\eta _\epsilon |^2(S_\epsilon 𝛀\times 𝐫)(iv,v)`$
$`={\displaystyle _{𝒟_\epsilon ^{}_iB_i}}𝛀(iv,dX_\epsilon \times v)(1+o(1))`$
$`={\displaystyle \underset{i}{}}2\pi \mathrm{\Omega }d_iX_\epsilon (p_i)(1+o(1))`$
$`={\displaystyle \underset{i}{}}{\displaystyle \frac{\pi \mathrm{\Omega }d_i}{1+\lambda ^2}}(\alpha |x_i|^2\lambda ^2|y_i|^2)^2(1+o(1)).`$ (17)
Finally, one can derive from (II.2)-(II.3)-(17) an asymptotic development of the energy for a solution with vortices.

$`E_\epsilon (u_\epsilon )E_\epsilon (\eta _\epsilon )2\pi |\mathrm{log}\epsilon |{\displaystyle \underset{i}{}}|d_i|a(p_i)`$
$`2\pi {\displaystyle \underset{ij}{}}d_id_ja(p_i)a(p_j)\mathrm{log}|p_ip_j|`$
$`{\displaystyle \frac{\pi \mathrm{\Omega }}{1+\lambda ^2}}{\displaystyle \underset{i}{}}d_i\left(\alpha |x_i|^2\lambda ^2|y_i|^2\right)^2.`$ (18)
Note that the minimal energy for solutions without vortices in $`𝒟_\epsilon ^{}`$ is $`E_\epsilon (\eta _\epsilon )+O(\epsilon |\mathrm{log}\epsilon |)`$: it is not exactly $`E_\epsilon (\eta _\epsilon )`$ since $`\eta _\epsilon `$ is a minimizer without the constraint $`\eta _\epsilon =1`$, but it almost equals to $`E_\epsilon (\eta _\epsilon )`$ since $`_𝒟a=1`$ and $`|\eta _\epsilon |^2`$ approaches $`a`$ asymptotically.

## III Critical velocities

### III.1 Critical velocity for the existence of one vortex

Let $`u_\epsilon `$ be a minimizer of $`(P)`$ with one vortex at a point $`p`$ in $`𝒟`$ with coordinates $`(x,y)`$ and let $`\mathrm{\Delta }E_\epsilon `$ be the difference between $`E_\epsilon (u_\epsilon )`$ and the energy of a solution without vortex ($`E_\epsilon (\eta _\epsilon )+O(\epsilon |\mathrm{log}\epsilon |)`$):

$`\mathrm{\Delta }E_\epsilon =2\pi (\alpha x^2\lambda ^2y^2)`$
$`\left(|\mathrm{log}\epsilon |{\displaystyle \frac{\mathrm{\Omega }}{1+\lambda ^2}}(\alpha x^2\lambda ^2y^2)\right)(1+o(1)).`$ (19)
This expression has been obtained by Svidzinsky and Fetter SF  using a different method. The form of $`\mathrm{\Delta }E_\epsilon `$ allows the computation of two critical velocities $`\mathrm{\Omega }_s`$ and $`\mathrm{\Omega }_1`$ for the existence of vortices: $`\mathrm{\Omega }_s`$ is the velocity for which the solution with one vortex starts to be locally stable and $`\mathrm{\Omega }_1`$ for which it starts to be globally stable. For $`\mathrm{\Omega }<\mathrm{\Omega }_s`$, $`\mathrm{\Delta }E_\epsilon `$ is a decreasing function of $`|p|`$, the position of the vortex; $`|p|=0`$ is a local maximum of $`\mathrm{\Delta }E_\epsilon `$. For $`\mathrm{\Omega }_s<\mathrm{\Omega }<\mathrm{\Omega }_1`$, $`|p|=0`$ is a local minimum for $`\mathrm{\Delta }E_\epsilon `$. Note that $`\mathrm{\Delta }E_\epsilon (pD)=0`$ and $`\mathrm{\Delta }E_\epsilon (|p|=0)>0`$.

For $`\mathrm{\Omega }>\mathrm{\Omega }_1`$, $`|p|=0`$ is the global minimum for $`\mathrm{\Delta }E_\epsilon `$. :

$`\mathrm{\Omega }_s={\displaystyle \frac{1+\lambda ^2}{2\alpha }}|\mathrm{log}\epsilon |={\displaystyle \frac{1+\lambda ^2}{4\sqrt{\lambda }}}\sqrt{2\pi }|\mathrm{log}\epsilon |,`$
$`\mathrm{\Omega }_1={\displaystyle \frac{1+\lambda ^2}{\alpha }}|\mathrm{log}\epsilon |={\displaystyle \frac{1+\lambda ^2}{2\sqrt{\lambda }}}\sqrt{2\pi }|\mathrm{log}\epsilon |`$
that is

$$\stackrel{~}{\mathrm{\Omega }}_s=\frac{1+\lambda ^2}{4\sqrt{\lambda }}\omega \sqrt{\frac{\pi \mathrm{}^2}{Ngm}}\mathrm{log}\left(\frac{Ngm}{\mathrm{}^2}\right)^{1/2},$$

$$\stackrel{~}{\mathrm{\Omega }}_1=\frac{1+\lambda ^2}{2\sqrt{\lambda }}\omega \sqrt{\frac{\pi \mathrm{}^2}{Ngm}}\mathrm{log}\left(\frac{Ngm}{\mathrm{}^2}\right)^{1/2}.$$

Note that Castin-Dum CD  for the case $`\lambda =1`$ find $`\stackrel{~}{\mathrm{\Omega }}_1=\omega \sqrt{(\pi \mathrm{}^2)/(Ngm)}\mathrm{log}\left((C/\sqrt{\pi })(Ngm)/\mathrm{}^2\right)^{1/2}`$, with $`C1.8`$, hence $`C/\sqrt{\pi }1`$ which gives a value of $`\stackrel{~}{\mathrm{\Omega }}_1`$ very close to ours. They also have $`\mathrm{\Omega }_1=2\mathrm{\Omega }_s`$ for the case $`\lambda =1`$. Increases in anisotropy yield higher $`\stackrel{~}{\mathrm{\Omega }}_1`$ as already noticed in FCS , but as $`\lambda `$ tends to infinity, $`\mathrm{\Omega }_1`$ becomes bigger than $`1/\epsilon `$ so that vortices cannot be stabilized.

It can be proved that there exists $`k_\epsilon `$ which tends to zero with $`\epsilon `$ such that for $`\mathrm{\Omega }<\mathrm{\Omega }_1k_\epsilon `$, the minimizer of $`E_\epsilon `$ has no vortex and for $`\mathrm{\Omega }>\mathrm{\Omega }_1+k_\epsilon `$, there exists a minimizer with a vortex. Moreover, for $`\mathrm{\Omega }_1+k_\epsilon <\mathrm{\Omega }<\mathrm{\Omega }_1+O(1)`$, any minimizer has one vortex of degree 1 tending to the origin. The proof consists in constructing a test function with a vortex at the origin and computing the energy of this test function. This yields an upper bound for the energy. The lower bound relies on estimates for $`G_\epsilon (v)`$ from BBH  and S1 ; S .

### III.2 Critical velocity for $`n`$ vortices

Similarly, one can compute $`\mathrm{\Omega }_n`$, the critical velocity for the existence of $`n`$ vortices. For this purpose, one can prove that as $`\epsilon `$ goes to 0, the vortices tend to the origin and they are of degree 1 that is $`d_i=1`$. This is similar to S1 , S . The test function consists in putting the $`n`$ vortices on a polygon centered at the origin of size $`1/\sqrt{\mathrm{\Omega }}`$ in $`x`$ and $`1/\lambda \sqrt{\mathrm{\Omega }}`$ in $`y`$. More precisely, we let $`\stackrel{~}{p}_i`$ with coordinates $`(\stackrel{~}{x}_i,\stackrel{~}{y}_i)`$ be such that $`\stackrel{~}{x}_i=x_i\sqrt{\mathrm{\Omega }}`$ and $`\stackrel{~}{y}_i=\lambda y_i\sqrt{\mathrm{\Omega }}`$. This allows us to estimate the energy of a solution with $`n`$ vortices centered at $`\stackrel{~}{p}_i`$ from (II.2),(II.3),(17):

$`E_\epsilon (u)=E_\epsilon (\eta _\epsilon )+2\pi n\alpha (|\mathrm{log}\epsilon |{\displaystyle \frac{1}{1+\lambda ^2}}\mathrm{\Omega }\alpha )`$
$`+\pi (n^2n)\alpha ^2\mathrm{log}\mathrm{\Omega }+w(\stackrel{~}{p}_1,\mathrm{},\stackrel{~}{p}_n)+C_n+o(1)`$ (20)
where $`C_n`$ is a constant that depends on $`n`$ and $`\lambda `$ and

$`w(\stackrel{~}{p}_1,\mathrm{},\stackrel{~}{p}_n)=\pi \alpha ^2{\displaystyle \underset{ij}{}}\mathrm{log}\left(|\stackrel{~}{x}_i\stackrel{~}{x}_j|^2+{\displaystyle \frac{|\stackrel{~}{y}_i\stackrel{~}{y}_j|^2}{\lambda ^2}}\right)`$
$`+2\pi \alpha {\displaystyle \underset{i}{}}(\stackrel{~}{x}_i^2+\stackrel{~}{y}_i^2)\left({\displaystyle \frac{2}{1+\lambda ^2}}{\displaystyle \frac{|\mathrm{log}\epsilon |}{\mathrm{\Omega }\alpha }}\right).`$ (21)
Recall that $`\alpha ^2=2\lambda /\pi `$. For fixed $`\lambda `$, $`w`$ is of order 1, hence is of lower order than the previous terms. Then the critical velocity for the existence of $`n`$ vortices can be computed from (20)

$`\mathrm{\Omega }_n=(1+\lambda ^2)\sqrt{{\displaystyle \frac{\pi }{2\lambda }}}|\mathrm{log}\epsilon |`$
$`+{\displaystyle \frac{(1+\lambda ^2)}{2}}(n1)\mathrm{log}\left((1+\lambda ^2)\sqrt{{\displaystyle \frac{\pi }{2\lambda }}}|\mathrm{log}\epsilon |\right),`$ (22)
and the critical velocity in the original parameters is $`\stackrel{~}{\mathrm{\Omega }}_n`$

$`\stackrel{~}{\mathrm{\Omega }}_n={\displaystyle \frac{(1+\lambda ^2)}{2}}\omega (\sqrt{{\displaystyle \frac{\pi \mathrm{}^2}{Ngm\lambda }}}\mathrm{log}\left({\displaystyle \frac{Ngm}{\mathrm{}^2}}\right)^{1/2}`$
$`+(n1)\sqrt{{\displaystyle \frac{\mathrm{}^2}{2Ngm}}}\mathrm{log}\left({\displaystyle \frac{1+\lambda ^2}{2\sqrt{\lambda }}}\sqrt{2\pi }\mathrm{log}({\displaystyle \frac{Ngm}{\mathrm{}^2}})^{1/2}\right)).`$
### III.3 Location of vortices

Once $`\mathrm{\Omega }`$ is close to $`\mathrm{\Omega }_n`$, the location of the vortices is characterized by the configuration of points $`\{\stackrel{~}{p}_i\}`$ which minimizes the function $`w`$ given by (21). In non-dimensionalized variables, the points are given by $`Rp_i/\sqrt{2\pi |\mathrm{log}\epsilon |}`$. For convenience, we define

$$\rho =\sqrt{\frac{2\pi }{\lambda }}\left(\frac{2}{(1+\lambda ^2)}\frac{|\mathrm{log}\epsilon |}{\mathrm{\Omega }\alpha }\right).$$

Note that given the value of $`\mathrm{\Omega }_n`$ in (III.2), to leading order, $`\rho `$ is equal to $`\sqrt{2\pi /\lambda }/(1+\lambda ^2)`$. We use the value of $`\alpha `$ and $`\rho `$ to get a simplified expression for $`w`$:

$`w(\stackrel{~}{p}_1,\mathrm{},\stackrel{~}{p}_n)=2\lambda ({\displaystyle \underset{ij}{}}\mathrm{log}(|\stackrel{~}{x}_i\stackrel{~}{x}_j|^2+{\displaystyle \frac{|\stackrel{~}{y}_i\stackrel{~}{y}_j|^2}{\lambda ^2}})`$
$`+\rho {\displaystyle \underset{i}{}}(\stackrel{~}{x}_i^2+\stackrel{~}{y}_i^2)).`$ (23)
The critical points of $`w`$, and thus the vortex positions, satisfy:

$`\rho \stackrel{~}{x}_i={\displaystyle \underset{ji}{}}{\displaystyle \frac{\lambda ^2(\stackrel{~}{x}_i\stackrel{~}{x}_j)}{\lambda ^2|\stackrel{~}{x}_i\stackrel{~}{x}_j|^2+|\stackrel{~}{y}_i\stackrel{~}{y}_j|^2}},`$ (24)
$`\rho \stackrel{~}{y}_i={\displaystyle \underset{ji}{}}{\displaystyle \frac{\stackrel{~}{y}_i\stackrel{~}{y}_j}{\lambda ^2|\stackrel{~}{x}_i\stackrel{~}{x}_j|^2+|\stackrel{~}{y}_i\stackrel{~}{y}_j|^2}}.`$ (25)
An immediate observation is that

$$\underset{i}{}\stackrel{~}{x}_i=\underset{i}{}\stackrel{~}{y}_i=0.$$
(26)
By multiplying the equations with $`\stackrel{~}{x}_i`$ and $`\stackrel{~}{y}_i`$ respectively and adding the results together, one can obtain

$$\underset{i}{}(\stackrel{~}{x}_i^2+\stackrel{~}{y}_i^2)=n(n1)/(2\rho ).$$
(27)
Similarly, multiplying the equations with $`\stackrel{~}{y}_i`$ and $`\lambda ^2\stackrel{~}{x}_i`$ respectively and adding the results together, one gets

$$\rho (1\lambda ^2)\underset{i}{}\stackrel{~}{x}_i\stackrel{~}{y}_i=0.$$
(28)
Unlike equations (26,27) where the dependence on $`\lambda `$ is implicit, the equation (28) leads to a property

$$\underset{i}{}\stackrel{~}{x}_i\stackrel{~}{y}_i=0,\text{ for }\lambda 1.$$
(29)
The above observations lead to more precise predictions on the location of vortices.

For instance, in the case $`n=2`$ we get $`\stackrel{~}{x}_1=\stackrel{~}{x}_2`$ and $`\stackrel{~}{y}_1=\stackrel{~}{y}_2`$. For $`\lambda =1`$, we have an infinite set of solutions consisting in two points on the circle $`\stackrel{~}{x}_i^2+\stackrel{~}{y}_i^2=1/\rho `$, symmetric with respect to the origin. For $`\lambda 1`$, (29) leads to $`\stackrel{~}{x}_i\stackrel{~}{y}_i=0`$ for $`i=1,2`$, and we have a pair of solutions with either $`\stackrel{~}{x}_i=0`$, $`\stackrel{~}{y}_i=\pm \sqrt{1/2\rho }`$, or $`\stackrel{~}{y}_i=0`$, $`\stackrel{~}{x}_i=\pm \sqrt{1/2\rho }`$. Checking the corresponding values of $`w`$, we get that for $`\lambda 1`$, the minimizer of $`w`$ corresponds to having both vortices staying on the long-axis of the ellipse in the original scaling (that is on the $`x`$ axis if $`\lambda >1`$ and on the $`y`$ axis if $`\lambda <1`$). The estimate of the location is in agreement with the numerical solutions given later.

For the case $`n=3`$, we also get that the 3 vortices are on the long-axis of the ellipse: one centered at the origin whereas the other pair stays symmetrically on the long-axis with $`\stackrel{~}{x}_i=\pm \sqrt{3/2\rho }`$ if $`\lambda >1`$. Similar discussions can be carried out for other values of $`n`$. In general, for the case of $`n`$ vortices, we let $`R_1^2=n(n1)/2\rho `$ and $`\widehat{x}_i=R_1\stackrel{~}{x}_i`$, $`\widehat{y}_i=R_1\stackrel{~}{y}_i`$. It follows from (27) that the points $`\{\stackrel{~}{p}_i\}`$ are localized by the minimum of

$$\underset{ij}{}\mathrm{log}(|\widehat{x}_i\widehat{x}_j|^2+\frac{|\widehat{y}_i\widehat{y}_j|^2}{\lambda ^2}).$$

under the constraint $`_i\widehat{x}_i^2+\widehat{y}_i^2=1`$. It is reasonable to expect that this formula will lead to a vortex array as observed in FCS2 . Note that the minimizer $`\{\widehat{p}_i=R_1\stackrel{~}{p}_i\}`$ has no explicit dependence on $`\mathrm{\Omega }`$, thus, we expect for a given $`\lambda `$, the vortex configuration to be of similar structure for values of $`\mathrm{\Omega }`$ close to $`\mathrm{\Omega }_n`$. On the other hand, for a given $`n`$ and $`\mathrm{\Omega }`$, we expect, however, that for sufficiently large $`\lambda `$, that is for highly anisotropic traps, the minimizer of $`w`$ is given by the collinear solutions with vortices all located on the long axis of the elliptical trap.

To our knowledge, theoretical investigation had been restricted to the critical velocity for nucleation of one vortex. Here, we are able to deal with the case of multiple vortices and to precisely characterize the location of the vortices. Our estimates are also consistent with the numerical results given later.

## IV Evolution equation and numerical schemes

To numerically compute the energy minimizers of (5), we notice that the energy in (5) can be rewritten as

$$_𝒟\left\{|(i𝐀)u|^2+\frac{1}{2\epsilon ^2}(a_\epsilon (𝐫)|u|^2)^2\right\}+c_\epsilon $$
(30)
where $`a_\epsilon (𝐫)=a(𝐫)\epsilon ^2\mathrm{\Omega }^2r^2`$,

$$𝐀=\left(\begin{array}{c}y\\ x\end{array}\right)𝛀,\text{ and }c_\epsilon =_𝒟\{\frac{1}{2\epsilon ^2}(a^2(𝐫)a_\epsilon ^2(𝐫))\}.$$

The above formulation of the energy has a striking similarity with the high-kappa high-field Ginzburg-Landau energy CDGP  with a variable coefficient DGP2 .

### IV.1 Evolution in the imaginary time

To numerically compute the minimizers of (5), we consider the time-dependent equation in the imaginary time:

$$\frac{u}{t}\left(i𝐀\right)^2u+\frac{1}{\epsilon ^2}|u|^2u\frac{a_\epsilon (𝐫)}{\epsilon ^2}u=\mu _\epsilon (u)u$$
(31)
in $`𝒟`$ with initial condition $`u(𝐫,0)=u_0(𝐫)`$ in $`𝒟`$ and boundary condition $`u=0`$ on $`𝒟`$. Here, $`\mu _\epsilon (u)`$ denotes the Lagrange multiplier. Assume that $`u_0`$ satisfies the constraint $`u_0=1`$. Then, by taking

$$\mu _\epsilon (u)=_𝒟\left\{\left|\left(i𝐀\right)u\right|^2+\frac{1}{\epsilon ^2}|u|^4\frac{a_\epsilon (𝐫)}{\epsilon ^2}|u|^2\right\}𝑑𝒟,$$

we get

$$\frac{d}{dt}\left(_𝒟|u|^2𝑑𝒟1\right)\mu _\epsilon (u)\left(_𝒟|u|^2𝑑𝒟1\right)=0.$$

Thus, the constraint $`_𝒟|u|^2=1`$ is ensured at all time. Moreover, using $`(u,u_t)=0`$, we get the energy estimate:

$$\frac{1}{2}\frac{d}{dt}_\epsilon (u)+\frac{u}{t}^2=0.$$

Thus, we easily get that for any $`(0,T)`$, if $`u_0H_0^1(𝒟)`$ and $`u_0=1`$, there exists a unique strong solution $`u`$ of (31) satisfying the constraint $`u=1`$. Using argument similar to that in LS ; LD , we may also get that as $`t\mathrm{}`$, $`u`$ approaches to a steady state solution which is a critical point of the energy. Well-posedness for $`L^2`$ initial data may also be be obtained.

### IV.2 Numerical schemes

There are various ways to solve the time-dependent Gross-Pitaevskii equations, see for example CD  or FCS . We take the advantage of the striking similarity with the high-kappa high-field time-dependent Ginzburg-Landau equations DG , and adapt a code developed in DG ; DGP ; DGP2 . Spatially, we use a standard finite element approximation, see DGP ; DGP2  for details. Here, we focus on the time-discretization and the treatment of constraint. It has been observed that there are some steady states exhibiting meta-stability, thus it is important to get asymptotically stable schemes for large time which in general requiring the use of implicit schemes with no limitations on the time step size.

Let $`\{u_n\}`$ be approximate solutions of $`\{u(t_n)\}`$ at discrete time $`\{t_n\}`$ with time-step $`\mathrm{\Delta }t_n=t_nt_{n1}`$. We discuss two time-discretization schemes and also some results of numerical experiments.

#### A first order backward-Euler in time discretization

Given $`u_{n1}`$, we first solve for $`u^{}`$:

$`{\displaystyle \frac{u^{}u_{n1}}{\mathrm{\Delta }t_n}}\left(i𝐀\right)^2u^{}\mu (u_{n1})u^{}`$
$`+{\displaystyle \frac{1}{\epsilon ^2}}|u^{}|^2u^{}{\displaystyle \frac{1}{\epsilon ^2}}a_\epsilon u^{}=0`$ (32)
Then, we apply the projection $`u_n=u^{}/u^{}`$. Both the backward Euler step and the projection step gives only first order in time accuracy.

#### A norm-preserving, energy-decreasing second order scheme

For any $`u`$, $`v`$ and their complex conjugate $`\overline{u}`$, $`\overline{v}`$, we let $`f(u,v)=(|u|^2+|v|^2)(u+v)/2`$ which satisfies $`f(u,v)(\overline{u}\overline{v})=(|u|^4|v|^4)/2`$. Given $`u_{n1}`$, we first solve for $`u^{}`$:

$`{\displaystyle \frac{2(u^{}u_{n1})}{\mathrm{\Delta }t_n}}\left(i𝐀\right)^2u^{}\nu (u^{})u^{}`$
$`+{\displaystyle \frac{1}{\epsilon ^2}}f(2u^{}u_{n1},u_{n1}){\displaystyle \frac{1}{\epsilon ^2}}a_\epsilon u^{}=0`$ (33)
where $`\nu (u^{})`$ is given by

$`\nu (u^{}){\displaystyle _𝒟}|u^{}|^2={\displaystyle _𝒟}\left\{\left|\left(i𝐀\right)u^{}\right|^2\right\}`$
$`+{\displaystyle \left\{\frac{1}{\epsilon ^2}f(2u^{}u_{n1},u_{n1})\overline{u}^{}\frac{a_\epsilon }{\epsilon ^2}|u^{}|^2\right\}}.`$
Then we let $`u_n=2u^{}u_{n1}`$. Taking the inner product of the equation (33) with $`u^{}`$, we get $`(u^{}u_{n1},u^{})=0`$, which leads to $`u_n^2=u_{n1}^2`$. That is, the norm is preserved at each time step. Taking the inner product of the equation (33) with $`u^{}u_{n1}`$, it is easy to get

$$2\frac{u_nu_{n1}^2}{\mathrm{\Delta }t_n}+_\epsilon (u_n)_\epsilon (u_{n1})=0.$$

Thus, during the discrete time evolution, the energy decreases. This discrete scheme is second order in time and unconditionally stable. It also preserves some essential features of the continuous dynamic system, making it suitable for long time integration and for studies of meta-stabilities of the solutions.

### IV.3 Description of the numerical experiments

We have used the above schemes to calculate various numerical solutions for the parameter values $`\epsilon =0.02`$, $`\lambda =1`$ and $`\lambda =1.5`$. The spatial finite element space is taken to be $`C^0`$ piecewise quadratic elements on triangular meshes. As we are mostly interested in the minimizers of the Gross-Pitaevskii energy, the time-evolution is employed as a means of marching to the steady state solutions. For this reason, we have used variable time-steps in order to accelerate the convergence in time. The nonlinear systems are solved by a Newton like methods at each step. Though such a method is computational costly per step, this drawback is offset by its unconditional asymptotic stability for marching to the steady state. We have also computed the solutions using refined meshes to ascertain the numerical convergence.

To obtain solutions for various velocities, we have used a number of differential initial conditions. For example, we have used $`|u_0(𝐫)|^2=a(𝐫)`$ for $`𝐫𝒟`$ which serves as a good approximation to the steady state solution, especially in the case of vortex-free solutions. We note that for large values of $`\mathrm{\Omega }`$, this choice of initial condition can also lead to steady state solutions with multiple vortices. Detailed solution branches are described in the next section. In addition, we have also used other initial conditions that manually plant vortices in the domain in order to find different solution branches. Finally, a continuation in the parameter $`\mathrm{\Omega }`$ has often been employed to follow a particular solution branch and to compute the bifurcation diagrams. The continuation procedure also provides a test for the local stability of the numerical solutions: when one branch becomes unstable, the solution jumps onto a different branch.

## V Numerical results and bifurcation diagrams

We now present some pictures of numerical solutions and discuss the various solution branches.

### V.1 Description of solutions

#### The case of a disc

For any $`\mathrm{\Omega }`$, there is a vortex free solution, which is close to $`a(x,y)`$ except near the boundary. For $`\mathrm{\Omega }=0`$, in addition to this vortex free solution, there is also a one vortex solution, as illustrated in the first column of Figure 1.

For larger $`\mathrm{\Omega }`$, solutions with multiple vortices are shown in the other columns of Figure 1. For instance, solutions with 2 and 3 vortices for $`\mathrm{\Omega }=15`$ are shown in the second column, solutions with 4 and 5 vortices for $`\mathrm{\Omega }=17.5`$ are shown in the third column, and solutions with 6 and 7 vortices for $`\mathrm{\Omega }=20`$ are shown in the last column. The parameter values for which the single-vortex and multi-vortices solutions exist are to be presented in the next section. Figures 2 and 3 provide the surface plots (and better view) of $`|u|`$ for a vortex free solution at $`\mathrm{\Omega }=0`$ and a solution at $`\mathrm{\Omega }=20`$ with four vortices respectively. Each solution has a top and bottom view, the paraboloid shapes are easy to visualize from the picture.

#### The case of an ellipse

We now present some solutions for an ellipse corresponding to $`\lambda =1.5`$. In Figure 4, the contour plots of the magnitude of the solution $`|u|`$ with $`\mathrm{\Omega }=17.5`$ are drawn while in Figure 5, $`\mathrm{\Omega }=25`$. For $`\mathrm{\Omega }=17.5`$, it is interesting to compare the location of vortices with the analysis of section III.c: in this case $`\rho =0.79`$, the long axis is 0.9887 and we find for the location of the vortex $`x=0.19`$ for $`n=2`$ and $`x=0.33`$ for $`n=3`$. The picture shows that the vortex is about 1/5 of the long axis for 2 vortices and 1/3 for 3 vortices, which is consistent with the analysis. We believe that for a bigger number of vortices $`n`$, their location corresponds to the positions minimizing (23).

Note that there were two vortex configurations with five vortices for $`\mathrm{\Omega }=25`$, among which, the non-symmetric configuration corresponds to a solution with a lower free energy, though the difference in their energy values is very small.

### V.2 Branches of solutions

An issue that we have studied is the existence of branches of $`n`$ vortex solutions as $`\mathrm{\Omega }`$ is varied and especially which one is the minimizer.

#### The case of the disc

For $`\mathrm{\Omega }=0`$, the solution with lowest energy is the vortex free solution. We start with this solution as initial value for the time dependent problem for a slightly bigger $`\mathrm{\Omega }`$. This device allows us to continue the branch of the vortex free solution as $`\mathrm{\Omega }`$ is increased. We find that the vortex free solution is obtained as the limit when $`t`$ is large of the evolution equation up to $`\mathrm{\Omega }=19`$. For $`\mathrm{\Omega }=19`$, six vortices are nucleated from the boundary and eventually one vortex moves to the center and final configuration is similar to that in Figure 1 (top right). Now if $`\mathrm{\Omega }`$ is decreased from the value 19 using the 6 vortices solutions as initial value, we see that this 6 vortices solution branch exists down to $`\mathrm{\Omega }=16`$ when it jumps to 4 vortices. If we decrease $`\mathrm{\Omega }`$ further, then we stay on the 4 vortex branch down to $`\mathrm{\Omega }=13`$ when it drops to a 2 vortex branch. Similarly, if we increase $`\mathrm{\Omega }`$ from 13, the 2 vortices solution will persist up to $`\mathrm{\Omega }=21`$.

As for the one-vortex solution branch, it is computed by planting a vortex-like function at the center in the initial condition for $`\mathrm{\Omega }=10`$, then the branch is computed by continuation in $`\mathrm{\Omega }`$ (both increasing and decreasing). It is interesting to observe the fact that this branch extends all the way to $`\mathrm{\Omega }=0`$, such persistence of the one-vortex solution even for the zero velocity has been elaborated by various authors, see FCS  for instance. Since implicit integration is used, we are able to indeed march to to the steady state and to ascertain that this persistence is not due to the metastability. On the other hand if the vortex is planted away from the origin, it disappears for small $`\mathrm{\Omega }`$ as we will see later.

#### The case of an ellipse

The same kind of behaviour occurs for the vortex free solution branch with $`\lambda =1.5`$: when $`\mathrm{\Omega }`$ is increased from 0, we stay on the vortex free branch up to $`\mathrm{\Omega }=22.5`$ when the solution jumps to the 4 vortex branch of Figure 5. If we decrease $`\mathrm{\Omega }`$ starting from this 4 vortex branch, the solution will stay on it down to $`\mathrm{\Omega }=15`$ when it jumps to 2 vortices and at $`\mathrm{\Omega }=10`$ it jumps to the vortex free solution. Similarly to the case of the disc, if one vortex is planted at the center at time 0, it will persist in time even down to $`\mathrm{\Omega }=0`$.

For $`\mathrm{\Omega }`$ large enough, several vortices are nucleated at the same time. Basically, we are counting on initial conditions and the round-off errors to break the symmetry of the region and strong symmetry presence often makes symmetry breaking much harder to achieve. For the disc case, we expect equal chance of vortex nucleation from any point of the boundary. It turns out that for $`\mathrm{\Omega }=19`$, an unstable front produces oscillations of almost equal magnitude, and spins off 6 vortices at the same time. Had the disturbance being unevenly distributed, then it would be possible for some vortices to get spin off ahead of others, and due to the repulsion, the others may never have a chance to appear, thus we may see solutions with a smaller number of vortices.

### V.3 Energy diagrams

We now discuss the energy diagrams in relation to the discussion of the critical velocities given in the earlier sections. In Figures 6 and 7, we have plotted the energy given by (3) as a function of $`\mathrm{\Omega }`$ for the various branches of solutions (according to the number of vortices). Again, $`\epsilon =0.02`$ in our computation.

#### The case of the disc (Figure 6)

As discussed earlier, the vortex-less solution, in the case of the disc, exists for all values of $`\mathrm{\Omega }`$ and is independent of $`\mathrm{\Omega }`$, we thus have a constant line for its energy. These vortex-free solutions are the global energy minimizers for small $`\mathrm{\Omega }`$ ($`\mathrm{\Omega }<9.3`$) whereas for $`\mathrm{\Omega }>9.3`$, they have larger energy than the one-vortex solution. For multi-vortex solutions, we see that each becomes the global energy minimizer for a range of values of $`\mathrm{\Omega }`$. It is interesting to compare this result with the value found in section 2 where $`\mathrm{\Omega }_1=9.8`$. Similarly, we obtain from (III.2) that $`\mathrm{\Omega }_n\mathrm{\Omega }_{n1}1.77`$. Though the value of $`\mathrm{\Omega }_1`$ is slightly overestimated, the difference $`\mathrm{\Omega }_n\mathrm{\Omega }_{n1}`$ looks good for small $`n`$: the numerics indicate $`\mathrm{\Omega }_2=12.0`$, $`\mathrm{\Omega }_3=13.6`$, $`\mathrm{\Omega }_4=15.8`$ and our theoretical computations yield $`\mathrm{\Omega }_2=11.6`$, $`\mathrm{\Omega }_3=13.4`$, $`\mathrm{\Omega }_3=15.2`$.

However, when $`\mathrm{\Omega }`$ is increased from 0, we saw that we stay on the vortex free branch up to $`\mathrm{\Omega }=19`$. This means that the vortex free solution is a local minimizer up to this value. We do not have any theoretical estimate for this value of local minimum. Similarly, for multi-vortex solutions, hysteresis loops are present. For the solution with six vortices, there are two possible configurations, one with all six on a concentric circle, one with only five on a concentric circle while the remaining vortex at the center of the disc. This occurs for $`\mathrm{\Omega }=22.5`$. Though the difference in the energy is hardly noticeable, the solution with a center vortex does have a smaller energy value.

#### The case of an ellipse (Figure 7)

Now, we discuss the elliptical case with $`\lambda =1.5`$. The energy versus $`\mathrm{\Omega }`$ curves are given in Figure 7.

The vortex-less solution, in the case of the ellipse, is no longer independent of $`\mathrm{\Omega }`$ as in the case of the disc, as illustrated by the dependence its energy on $`\mathrm{\Omega }`$. Aside from that, the hysteresis phenomena also occurs just like for the case of a disc.

Our numerics indicate that $`\mathrm{\Omega }_1=12`$, $`\mathrm{\Omega }_2=15`$. It is interesting to compare this result with the value found in section 2: we obtain that $`\mathrm{\Omega }_113`$ from (III.2) that $`\mathrm{\Omega }_n\mathrm{\Omega }_{n1}2.88`$. Though the value of $`\mathrm{\Omega }_1`$ is slightly overestimated, the difference $`\mathrm{\Omega }_n\mathrm{\Omega }_{n1}`$ looks good again since we find $`\mathrm{\Omega }_2=14.9`$.

### V.4 Displacement of the vortex from the center

Based on the earlier estimate (III.1) for the one-vortex solution, we see that for small $`\mathrm{\Omega }`$, the displacement of the vortex away from the center leads to the drop of energy ($`\mathrm{\Omega }<4.9`$). For slightly bigger $`\mathrm{\Omega }`$ a vortex at the center is a local minimum ($`\mathrm{\Omega }<9.8`$). Let us analyze the time dependent equation using a vortex off center as initial condition and let us examine how it evolves. We find that in the case of the disk, for $`\mathrm{\Omega }<7`$, a displacement on the size of a tenth of the radius causes the displaced center vortex in the one-vortex solution to move towards the boundary. For $`\mathrm{\Omega }>10`$, the center vortex in the one-vortex solution moves back towards the center if under a displacement on the size of a third of the radius. For intermediate values of $`\mathrm{\Omega }`$, the vortex moves back to the center for small displacement but moves towards boundary for large displacement.

The following two pictures (Figures 8 and 9) show, in the case of a disc and a small displacement, the marching away of the center vortex for the solution with $`\mathrm{\Omega }=0`$ (starting from top row then to the bottom row) and the marching back to the center for the solution with $`\mathrm{\Omega }=7.9`$.

## VI Conclusion

We have presented a new framework for the study of the Gross Pitaevskii energy in the Thomas Fermi limit: we have defined a nondimensionnalized parameter $`\epsilon `$ and we have computed theoretically an asymptotic development of the energy, the critical thermodynamic velocities of nucleation of vortices and the location of vortices. This extends the results of CD . We have proposed and implemented time integration schemes which enjoy the norm-preserving and energy-descreasing features and thus useful for studying the stability and metastability of solutions. We have also presented energy diagrams computed numerically for the various vortex solutions. We have noticed that our theoretical predictions for critical thermodynamic velocities are quite consistent with the numerics which encourages us to think that our approximation $`\epsilon `$ small is correct though $`|\mathrm{log}\epsilon |`$ is not small.

In our computation, we have taken $`\epsilon =0.02`$ and the number of vortices we have observed ranges from 0 to 10 for $`\mathrm{\Omega }`$ between $`0`$ and $`25`$.

Our aim for a future work is to use the results of PR  for Ginzburg-Landau in dimension 3 to get results for rotating Bose Einstein condensates in dimension 3, and especially the critical velocities and an asymptotic development of the energy. In particular, it is obtained in SF  that $`\mathrm{\Omega }_s=3/5\mathrm{\Omega }_1`$ for a vortex line parallel to the $`z`$ axis, but we would like to get an expression of the energy for any type of vortex line. We will also carry out further numerical simulations in 3 dimension to compare with experimental data.

###### Acknowledgements.

The authors would like to thank Y.Castin for very interesting discussions and also T.Riviere, E.Sandier, S.Serfaty. This work is supported in part by a joint France-Hong Kong research grant.