Superfluids Passing an Obstacle and Vortex Nucleation

We consider a superfluid described by the Gross-Pitaevskii equation passing an obstacle \[\epsilon^2 \Delta u+ u(1-|u|^2)=0 \ \mbox{in} \ {\mathbb R}^d \backslash \Omega, \ \ \frac{\partial u}{\partial \nu}=0 \ \mbox{on}\ \partial \Omega \] where $ \Omega$ is a smooth bounded domain in $ {\mathbb R}^d$ ($d\geq 2$), which is referred as the obstacle and $ \epsilon>0$ is sufficiently small. We first construct a vortex free solution of the form $ u= \rho_\epsilon (x) e^{i \frac{\Phi_\epsilon}{\epsilon}}$ with $ \rho_\epsilon (x) \to 1-|\nabla \Phi^\delta(x)|^2, \Phi_\epsilon (x) \to \Phi^\delta (x) $ where $\Phi^\delta (x)$ is the unique solution for the subsonic irrotational flow equation \[ \nabla ( (1-|\nabla \Phi|^2)\nabla \Phi )=0 \ \mbox{in} \ {\mathbb R}^d \backslash \Omega, \ \frac{\partial \Phi}{\partial \nu} =0 \ \mbox{on} \ \partial \Omega, \ \nabla \Phi (x) \to \delta \vec{e}_d \ \mbox{as} \ |x| \to +\infty \] and $|\delta |<\delta_{*}$ (the sound speed). In dimension $d=2$, on the background of this vortex free solution we also construct solutions with single vortex close to the maximum or minimum points of the function $|\nabla \Phi^\delta (x)|^2$ (which are on the boundary of the obstacle). The latter verifies the vortex nucleation phenomena (for the steady states) in superfluids described by the Gross-Pitaevskii equations. Moreover, after some proper scalings, the limits of these vortex solutions are traveling wave solution of the Gross-Pitaevskii equation. These results also show rigorously the conclusions drawn from the numerical computations in \cite{huepe1, huepe2}. Extensions to Dirichlet boundary conditions, which may be more consistent with the situation in the physical experiments and numerical simulations (see \cite{ADP} and references therein) for the trapped Bose-Einstein condensates, are also discussed.


Introduction
This paper addresses the vortex nucleation and vortex shedding phenomena when a superfluid passes an obstacle. Of particular concern is the existence of associated steady state solutions of the following Gross-Pitaevskii equation passing an obstacle ǫ 2 ∆u + u(1 − |u| 2 ) = 0 in R d \Ω, ∂u ∂ν = 0 on ∂Ω (1.1) where ǫ > 0 is a sufficiently small constant, Ω is a smooth and bounded domain in R d , d ≥ 2 and ν denotes the unit outer normal.
The purpose of this paper is to construct two types of solutions to (1.1). This will be achieved by perturbation of two basic solution profiles. The first basic profile, which is a relatively clear one, is obtained through the so-called Madelung transformation, In other words, one is interested in solutions in the semiclassical regime (i.e. ǫ being sufficiently small), see [34] and references therein for discussions in evolutionary cases. Equation ( ( 1.4) There are classical works by L. Bers [8,9] in the two dimensional case, R.Finn-Gilbarg [19] and G. Dong [16] in higher dimensional cases. We summarize the basic results in the following theorem. Theorem 1.1. (i) There exists a δ * ∈ (0, 1) such that for |δ| < δ * , there exists a unique classical solution Φ = Φ δ (steady state solution of (1.4)). For |δ| > δ * , there are no classical solutions (the so-called shocks develop). (ii) The solution Φ δ has the property that for |δ| < δ * max x∈R d \Ω |∇Φ δ (x)| < 1 3 . (1.5) In the above theorem, δ * is called the sound speed of this problem. Thus for classical fluids passing an obstacle, the situation is relatively clear. For a fluid with a speed less than the sound speed, there is a smooth stationary flow. When the speed of fluids goes beyond a critical (sound) speed, there are no smooth stationary flows (shock develops). From the semiclassical limit formalism (see [34] and the references therein), one would expect a similar conclusion may be also true for superfluids passing an obstacle. Since superfluids are frictionless, there are no notions of "shock" in this case. Instead, there would be nucleation of vortices and the latter would introduce a dissipation mechanism that eventually destroy the superfluidity, see for examples [32,20,30,31,40,41,28,29]. However, there is no rigorous mathematical proof. Throughout this paper, we always assume that |δ| < δ * .
Though the conclusion of the above theorem may be expected (and it has been often assumed in many physics literature), it lacks of a rigorous mathematical proof. In fact, even at a low superfluid speed, we shall see that certain boundary layers near the obstacle may develop and this causes difficulties for analysis. See (2.5) below.
Let u ǫ = ρ ǫ e i Φǫ ǫ be the solution constructed in Theorem 1.2. It turns out that when d = 2 on the top of this solution a second solution exists. Interestingly the limiting profile of this second solution is the traveling wave solutions of Gross-Pitaevaskii equation. More precisely, set (1.7) Then v satisfies (coupled with homogeneous Neumann boundary condition) Let x 0 ∈ ∂Ω and perform a rescaling as follows: x = x 0 + ǫy. Formally letting ǫ → 0, (and after proper scaling), we obtain (assuming that ∇Φ δ (x 0 ) = |∇Φ δ (x 0 )| e 2 ) the following traveling wave equation coupled with the following boundary condition ∂U ∂y 1 (0, y 2 ) = 0. (1.10) We refer to Section 3.1 for more detailed derivations. Here (1.11) A simple computation shows that the subsonic condition (1.5) is equivalent to the following speed condition 0 < c < √ 2.
(1.12) The traveling wave problem (1.9)-(1.10) for Gross-Pitaevskii equation has been under study in many papers [2,7,21,22,23,24,33]. It has been proved that for c ≥ √ 2 there are no traveling waves ( [24]). When c < √ 2, variational method shows that there are traveling wave solutions to (1.9). For c small a perturbation argument can be used to show the existence ( [33]). The properties of solutions constructed in [33] would play an important role in the proof of our main theorem below. The asymptotic behavior and qualitative behavior of solutions are also studied in many papers [21,22,23,24]. What is remarkable and fascinating is the fact that the critical speed for existence of traveling wave solutions for the Gross-Pitaevskii equation on the entire plane is directly related to the critical (sound) speed for stationary flows of superfluids passing a smooth obstacle described by the Gross-Pitaevskii equation through an explicit but nonlinear algebra relation (1.11).
Our second result shows the traveling wave solutions to (1.9) persist for the superfluids problem (1.1), as long as |δ| is suitably small (see Section 3 below). Theorem 1.3. Let d = 2 and 0 < c < √ 2. Let U c be a solution of (1.9)-(1.10) satisfying a nondegeneracy condition. (See Key Assumption (3.16) below.) Then there exists c 0 > 0 and ǫ 0 such that for 0 < ǫ < ǫ 0 , |c| < c 0 problem (1.1) has at least two smooth solutions of the form where u ǫ is the solution given by Theorem 1.2 and U c is the traveling wave solution of problem (1.9)-(1.10). Here c is given by (1.11) Theorem 1.3 shows not only the phenomena of vortex nucleations in stationary flows of superfluids but also a somewhat surprising new phenomena that even before the "sonic" speed vortices can nucleate near the boundary of the obstacle. It does rigorously justify some seemly strange conclusions drawn from previous numerical studies in [26,27]. We believe that the second solution exists for all subsonic speed c < √ 2. But this remains to be an open problem and we wish to return to this issue later, see remarks in Section 3. Theorem 1.3 also shows the situation near the superfluid "sonic" speed may be much more complex than some formal studies done previously, [20,30,31,41], and the study of the latter situation would need a new set of tools and ideas.
As we mentioned earlier, the proofs of Theorems 1.2 and 1.3 are based on perturbations of two kinds of solutions. In the proof of Theorem 1.2 the primary ansatz is the solution to (1.4). The linearized operator is a system whose Fourier transforms are anisotropic. The additional difficulty is the existence of boundary layers. We use energy method and a priori estimates to prove Theorem 1.2. Theorem 1.3 is perturbed from a traveling wave solution to the Gross-Pitaevskii equation (1.9) in the whole space. Under a nondegenerate condition of solutions to (1.9), which we verify for |c| << 1, we prove Theorem 1.3 by finite dimensional Liapunov-Schmidt reduction method.
Throughout this paper, we always assume that d = 2 or 3 (which are the physical dimensions). (The result of Theorem 1.2 is likely to hold for d ≥ 4 as well.) The constant C is a positive generic constant independent of ǫ < ǫ 0 and c < c 0 . Denote B ρ (y) = {x | |x − y| < ρ}. We also use the following notation < y >:= 1 + |y| 2 , < f, g >= Re( fḡ).
Acknowledgments: The research of the first author is partially supported by the NSF grant DMS-1501000. The research of the second author is partially supported by NSERC of Canada.

Proof of Theorem 1.2
In this section, we prove the existence of vortex free solution, e.g. Theorem 1.2. First we introduce two nonlinear operators (2.1) Then equation (1.3) can be written in an operator form where C 2,α ν (R d \Ω) = C 2,α (R d \Ω) ∩ { ∂ρ ∂ν = 0 on ∂Ω} and α ∈ (0, 1) is the Hölder exponent. As mentioned in the introduction, formally letting ǫ = 0 in equation (2.2), we obtain the limiting problem (1.4). We first collect some additional properties of solutions to (1.4), which will be needed in later sections.
Proof. The asymptotic behavior can be found in [Theorem III, [8]] (in the case of d = 2) and [16] (in the case d = 3).

Approximation solution and boundary layer.
For the first approximation function, we take W 0 = (ρ δ , Φ δ ) where Φ δ is the solution given in Theorem 2.1 and ρ δ = 1 − |∇Φ δ | 2 . It is easy to see that We observe that for this initial approximate solution W 0 , the second component satisfies the Neumann boundary condition but the first component does not, which has to be corrected by a boundary layer. To this end, we now add a correction function to the first component: let ρ 1 be the unique solution of Observe that by the estimates for |∇Φ δ | in Lemma 2.1, it holds that On the other hand, using classical barrier function, we have the following estimates for ρ 1 : Let us choose the second approximation as follows We then compute . (2.10) Now we linearize around W 1 as follows and obtain ). (2.12) Observe that N 2 is of the form An important observation is that if ∂ρ ∂ν = ∂Φ ∂ν = 0 on ∂Ω, then it holds that g · ν = 0 on ∂Ω. (2.14) We aim to solve the following system of equations (2.16) The computations above provides basic estimates to proceed in the next steps.

A priori estimates.
To proceed with the perturbation, we need the following important a priori estimates.
2.5. Remarks on further estimates of (ρ 2 , Φ 2 ) near the boundary. For later purpose, we need to expand (ρ 2 , Φ 2 ) near a boundary point. It turns out that the estimate is better. Let d = 2 and x 0 ∈ ∂Ω. We may assume that the normal direction is ν x0 = e 1 , the tangential direction is τ x0 = e 2 and R d \Ω ⊂ {x 1 ≤ 0}. Rescale x = x 0 + ǫȳ. Then we see that For the first term, we have This implies that Since the remaining errors in (2.16) carries at least ǫ 2 order, we conclude that (2.32) Therefore we may take σ = 1 + σ 0 ∈ (1, 2) in the proof of Theorem 1.2 to obtain that for some σ 0 > 0. (In the computations below we may simply assume that ρ ǫ = ρ δ , Φ ǫ = Φ δ .) Now we look for another solution of the following form We see that v satisfies In the following, we explain the intuitive idea in the construction of the second solution. Let x 0 ∈ ∂Ω (to be determined later) and we perform a blow-up near (2.33)). (Without loss of generality we may assume that which is equivalent to the subsonic assumption (1.5) in Theorem 2.1. Problem (3.4) arises as the traveling wave solution u(ȳ − ct e d ) for the Gross-Pitaevskii equation It has been proved that a necessary condition for the existence of traveling waves is |c| < √ 2. In the case of small speed |c| << 1, the existence of traveling vortex (d = 2) and vortex rings (d ≥ 3) has been proved by Bethuel-Saut [2], Bethuel-Orlandi-Smets [3], and [33]. The method of [2,3] is variational while we used a perturbation approach which is related to what we will employ in this paper. The asymptotics of traveling wave solutions is studied in [21]- [24]. For general speed c, we refer to Bethuel-Gravajat-Saut [4], Gravajat [21,22,23,24] and references therein.

3.2.
Properties of traveling waves solutions. In this sub-section, we are concerned with the properties of the traveling wave solution to (3.4) in R 2 . So from now on, we assume that d = 2 and we consider In [33], we used a perturbation approach to prove the existence of a traveling wave solution to (3.7) with two opposite vortices traveling in the direction y 2 . We summarize the result in the following (ii) U c is even in y 1 , i.e., ∂U c ∂y 1 = 0 on y 1 = 0; (3.10) (iii) Writing U c (y) = S(y)e iϕ(y) , then it holds that Proof. Properties (i)-(iii) follow from the constructions given in [33]. (In [33], Schrodinger map is studied. But the same computations work for (3.7) as well.) For the asymptotics (3.11), it also follows from [21].
The following theorem give a complete characterization of the kernels of the linearized operator, at least when the speed c is small. The proof of it will be delayed to Section 6.
Theorem 3.1. Let U c be the solution constructed in [33]. Consider the linearized operator

12)
Then for c sufficiently small, the only bounded kernels satisfying where β 1 and β 2 are some constants.
In the general large c case, the following result is proved in several papers ([2]- [4]).
Proof. (1) and (2) are proved in [7] and (3) is proved in [21]. (Note that the exact asymptotics of U c − 1 is not precised in dimension two.) By Theorem 3.1, there exists a nondegenerate solution U c for c small. From now, we assume the following: Key Assumption: There exists a nondegenerate solution U c for c ∈ (0, c 0 ) for some c 0 ∈ (0, √ 2]. This means that there exists a solution to (3.7) such that the only solutions to the linearized problem (3.13) are where β 1 and β 2 are some constants.
Remark 3.1. For c small, there are higher energy solutions with l(l + 1) traveling vortices. These vortices are located at the roots of Adler-Moser polynomials. The existence and properties of these multiple vortices solutions are considered in [35]. We believe that these higher energy solutions are also nondegenerate, which may provide new solutions to equation (1.1).

3.3.
Perturbation of traveling wave solution and sketch of proof of Theorem 1.3. Under the Key Assumption, in the remaining sections, we use a finite dimensional reduction method to rigorously prove the existence of such solutions to (3.2) and hence give the proof of Theorem 1.3. The important step is the determination of x 0 ∈ ∂Ω. We solve (3.2) in three steps.
Step 1: Fixing each x 0 ∈ ∂Ω we show that there exists a unique solution to (3.2) and two Lagrange multipliers where Z 0 and Z 1 are defined at (4.27).
By placing x 0 near the maximum or the minimum of |∇Φ δ | 2 (x 0 ) on the ∂Ω, we can adjust x 0 such that λ 1 = 0. (Hence we always obtain at least two solutions, as stated in Theorem 1.3.) Step 3: In the last step we use the gauge invariance of the equation to show that λ 0 = 0.
In the following sections, we carry out this procedure. At first we need to understand the error and invertibility of the linearized operator.

Proof of Theorem 1.3
In this section we follow the steps outlined above to prove Theorem 1.3.

4.1.
Approximate solutions and error estimates. Let x 0 ∈ ∂Ω and x = x 0 +ǫy. (x 0 is a priori undetermined.) Without loss of generality we also assume that the outer normal direction ν x0 = e 1 and the tangential direction τ x0 = e 2 . In this stretched variable equation (3.2) is transformed to Here the new stretched domain becomes and ∂v ∂ν = 0 on ∂Ω ǫ,x0 .
We write W (y) := U c (ȳ) = S(ȳ)e iϕ(ȳ) where U c is given by Theorem 3.1 andȳ = 1 − |∇Φ δ (x 0 )| 2 y. (Note that S(ȳ) is not radially symmetric.) To avoid clumsy notations, we also drop the dependence on the speed c and the location x 0 . Then W satisfies By Theorem 3.2, we have the following decaying estimate for S and ϕ: We write (4.1) as a solution operator Using (4.2) we get . (4.5) In this section, we estimate the size of the error S[W ]. The first term in (4.5) is expanded as For the first and second terms in (4.6), we make use of the decaying estimate (4.3) and decaying estimate (2.3) . The first term has the following decay estimates 2ǫ|∇ log ρ ǫ ∇S| ǫ < y > −3 (4.7) while for the second term gives 2ǫ|∇ log ρ ǫ ∇ϕ| ǫ < ǫy > −2 < y > −2 ǫ 1−σ (1 + |y|) −2−σ (4.8) where σ ∈ (0, 1) is any given positive number. For the second term in (4.6), we have By the decaying estimates in (4.3) we get that (4.9) We note that for |y| >> 1 Using (4.3) again we obtain the following basic error estimates Concerning the boundary behavior, we can strengthen the boundary and use the fact that ∂W ∂y1 (0, y 2 ) = 0 as well as the decaying (4.3) to deduce that (4.12)

Equations in operator form.
We look for solutions of (4.1) in the form of a small perturbation of W . Let η be a smooth cut-off function such that η(y) = 1 for |y| > R and η(y) = 0 for |y| > 2R. (Here R is a large and fixed constant.) As in [15,18] or [33] we look for solutions of (4.1) of the form v = η(W + iW ψ) + (1 − η)W e iψ (4.13) where W ψ is small. We write ψ = ψ 1 + iψ 2 with ψ 1 and ψ 2 real-valued. Set (4.14) where For |y| > 2R, we have v = W e iψ and thus by simple computations we get iW e iψ which in terms of ψ can be written as Recalling that ψ = ψ 1 + iψ 2 and x = x 0 + ǫy, we have 1 L ǫ,2 contains linear terms which will be shown to be higher order, whileÑ ǫ contains nonlinear terms in ψ. Let us remark that the explicit form of all the linear and nonlinear terms will be very useful for later analysis.
Finally equation (4.15) has to be solved with the following boundary condition 4.3. Weighted norms and error estimates. We first introduce some norms. Let us fix two small positive numbers 0 < γ < 1, 0 < σ < 1. Recall that φ = iW ψ, ψ = ψ 1 + iψ 2 . Let R be a fixed but large number so that |W | ≥ 1 2 for |y| > R. Now we define the following norms for (complex) The forms of these norms are motivated by the expressions of (4.17)-(4.18). We refer to [15,18] for similar definitions.
To this end, we first need to consider the following linear problem We have the following key a priori estimates.
Lemma 4.1. There exists a constant C, depending on γ, σ only such that for all ǫ sufficiently small, and any solution of (4.29), it holds Proof. We proceed similar to the proof of Lemma 5.1 in [33]. See also similar arguments in [?, 18].
The key difference is now that we relax the decay of ψ 1 to be bounded only and that we don't have the (odd) symmetry assumption. To overcome this difficulty we use some ideas from [14]. We prove it by contradiction. Suppose that there exist a sequence of ǫ = ǫ n → 0, constants λ n 0 , λ n 1 , and functions φ n , h n , g n which satisfy (4.29) with φ n * = 1, h n * * ,i = o(1), g n * * ,b = o(1). (4.31) We first note that λ n 0 , λ n 1 = o(1). This follows by multiplying the equation (4.29) with (iW )z 0 , (iW )z 1 and integrating by parts. See Lemma 6.2 of [14].
Next we derive inner estimates first. Let R > 0 be fixed and large. We claim that φ n L ∞ (B4R) = o(1). In fact suppose not. Then by a limiting process we obtain a solution to the linear equation L 0 in R 2 + . Thanks to the Key Assumption (3.16), the kernels of L 0 consist of linear combinations of z 0 and z 1 only. But the limit of φ n is exactly orthogonal to the approximate kernels. This is impossible. (This argument is standard now so we omit the details. See [15,18,33].) Next we shall derive outer estimates: this is the more technical part. (To avoid the clumsy notation we drop the dependence on n.) We follow the proof in Section 6 of [14]. For |y| > 2R the system becomes (see (4.16)) (4.32) where f 1 , f 2 contain all remaining error terms which are of the order o(< y > 2+σ ), o(< y > 1+σ ) respectively. The first two terms in the second equation in (4.32) behavior like ∆ − 1, which implies by Maximum Principle It remains to consider the first equation in (4.32) only. Similar to the estimates of Lemma 6.1 of [14], for linear equation In the first equation in (4.32) we consider the linear term 2∇ x Φ ǫ (x 0 ) · ∇ψ 2 as a perturbation, since |∇ x Φ ǫ | δ. Therefore we get where the last term can be bounded by δ φ * .
All together we obtain the following outer estimates Combining both inner and outer estimates, we obtain that φ * = o(1), a contradiction to (4.31).
Remark 4.1. The condition that δ is small is only used in the outer estimate argument. We believe that this is just a technical condition.
We consider now the following projected linear problem  We state the following existence result for the projected linear problem.
Proposition 4.1. There exists ǫ 0 such that for all ǫ < ǫ 0 , the following holds: if (h, g) * * < +∞, then there exists a unique solution φ = T ǫ (h, g) to (4.33). Furthermore it holds that Proof. The proof is similar to that of [Prop. 4.1, [15]]. Instead of solving (4.33) in R 2 \Ω ǫ,x0 , we solve it in a bounded domain first: where M > 10R. By the same proof of a priori estimates, we also obtain the following estimates for any solution φ = φ M of (4.35): By working with the Sobolev space H 1 (R 2 \Ω ǫ,x0 ∩ B M (0)) ∩ {u = 0 on ∂B M (0)}, the existence then follows from Fredholm alternatives. (Recall that problem (4.1) admits a variational structure as follows 1 2 Here the fact that ρ ǫ e iΦǫ ǫ satisfies (2.2) is used.) Now letting M → +∞, we obtain a solution to (4.33) with the required properties.

4.5.
Step 1: projected nonlinear problem. Finally, we consider the full nonlinear projected problem (4.28). By (4.15) this is equivalent to (4.37) Using the operator T ǫ defined by Proposition (4.1), we can write (4.37) as which is equivalent to where G ǫ is the nonlinear operator at the right hand side of (4.38).

Dirichlet boundary condition
We discuss in this section how we can adjust the proofs to deal with the Dirichlet boundary conditions. We consider the Gross-Pitaevskii equation with Dirichlet boundary condition ǫ 2 ∆u + u(1 − |u| 2 ) = 0 in R 2 \Ω, u = 0 on ∂Ω. (5.1) We first discuss the existence of vortex free solutions. Same as before we let u = ρe Φ ǫ . Then we have For the boundary conditions of Φ we impose the usual Neumann boundary condition The first ansatz is W 0 = (ρ δ , Φ δ ). Similar to the Neumann boundary condition case, we need to add a boundary layer ρ 1 : The remaining proofs are similar to the Neumann boundary condition case. We omit the details.
To construct the second solution, we need to analyze the behavior of the first solution near the boundary and find the corresponding limiting traveling wave equation.
Let us rescale x = x 0 + ǫy where x 0 ∈ ∂Ω,Φ = Φ ǫ . (As before we also assume that ν x0 = e 1 , τ x0 = e 1 .) Then we have There exists a solution to (5.5) of the following form ρ = ρ 0 (y 1 ),Φ = by 2 (5.6) where ρ 0 is the unique solution of the following ordinary differential equation We claim that for b small we can construct a new solution U b to (5.7) with two opposing vortices. As in [33] we take the initial ansatz the same as before The only new error in the equation comes from the interaction with the boundary layer ρ 0 which is the following ρ ′ 0 ρ 0 ( y 2 y 1 d ((y 1 − d) 2 + y 2 2 )((y 1 + d) 2 + y 2 2 ) ) Note that ρ ′ 0 ∼ e −Cy1 , near the vortex y ∼ (d, 0) it is exponentially small. The L 1 norm of this error has the order O( 1 d ). The rest of the perturbation arguments in [33] goes through. We omit the details.

Proofs of Theorem 3.1
In this section, we prove the key nondegeneracy result Theorem 3.1. First it is easy to see that the following functions satisfy the equation (3.12) and the Neumann boundary condition ∂φ ∂y1 (0, y 2 ) = 0. Hence they belong to the kernel (3.13). To prove the converse statement, we note that z 2 = ∂Uc ∂y1 satisfies the equation (3.12), however does not satisfy the Neumann boundary condition. We now show that this function produces instead a nonzero eigenvalue. To this end, we go back to the construction process. The existence of a solution to (3.7) for c small is proved in the following steps. To align with the proofs in [33], we use the notation c = ǫ and assume that ǫ > 0 is small. We first introduce some definitions from [33].