Gross-Pitaevskii vortex motion with critically-scaled inhomogeneities

We study the dynamics of vortices in an inhomogeneous Gross-Pitaevskii equation $i \partial_t u = \Delta u + {1\over \varepsilon^2} (p_\varepsilon^2(x) - |u|^2)$. For a unique scaling regime $|p_\varepsilon(x) - 1 | = O(|\log \varepsilon|^{-1})$, it is shown that vortices can interact both with the background perturbation and with each other. Results for associated parabolic and elliptic problems are discussed.


Introduction
We consider the behavior of vortices in an inhomogeneous Gross-Pitaevskii (IGP) equation where Ω is a bounded, simply connected domain in R 2 . An important feature of solutions of (1.1) is the potential appearance of vortices, points where u(x, t) : R 2 → C has a zero with nontrivial topological degree. We concentrate on the case where there are n vortices of degree d j ∈ {−1, 1} at positions a j (t) ∈ R 2 for j = 1, . . . , n. We will show below that solutions take the form, Here x−a j (t) |x−a j (t)| is interpreted in the complex sense as z−a j (t) |z−a j (t)| for z −a j = x 1 −a j,1 +i(x 2 −a j,2 ) with ρ ε (x) → 0 as |x − a j | → 0 and |u(x, t)| 2 ≈ p 2 (x) away from the a j (t)'s. In a sense identifying the location and degrees of vortices is sufficient to describe the solution u(x, t).
The inhomogenous Gross-Pitaevskii equation (1.1) arises in several contexts, including the behavior of superfluid 4 He near a boundary, Bose-Einstein Condensation (BEC) under a very weak confinement potential, and nonlinear optics. In the first case it has been observed experimentally that anomalous behavior of superfluid 4 He vortex filament motion indicates J.L.M. was support in part by NSF Applied Math Grant DMS-1312874 and NSF CAREER Grant DMS-1352353. D.S. was supported in part by NSF CAREER grant DMS-0955687 and NSF grant DMS-1516565. The authors wish to thank Sylvia Serfaty and Rafe Mazzeo for helpful conversations throughout the development of this work and wish to thank the anonymous referees for their comments that helped to improve the paper. that vortex tubes can become pinned to sites on the material boundary, see [46,38,45]. Such pinning can be incorporated heuristically by the introduction of a potential p(x) in the Gross-Pitaevskii equation, see [34] or complex Ginzburg-Landau equations, see [15]. In the second case, IGP is a useful model in a BEC in which the vortices interact both with each other and with the trap potential. In the third case, IGP is useful in modeling optical vortices. See for instance [3,2,34], which are related to motion of vortices with non-vanishing total charge in inhomogeneous potentials and confinement.
The dynamics of vortices in Gross-Pitaevksii equations have been the subject of much study in past few decades. The connection between the vortex motion in superfluids, where p(x) ≡ 1, and simplified ODE's was made by Fetter in [11], and it was shown that vortices interact with each through a Coulomb potential. On the other hand in Bose-Einstein condensates with nontrivial trapping potentials, p(x) = O(1), it was shown by Fetter-Svidzinsky in [12] that vortices interact solely with the background potential and are carried along level sets of the Thomas-Fermi profile.
On the other hand it has been observed, both experimentally [14,32] and numerically [29], that vortices do interact with other vortices and with the background potential. In one particularly interesting direction, experimentalists have generated dipoles by dragging non radially symmetric BEC's through laser obstacles in [32] and by the Kibble-Zurek mechanism in [14]. These bound dipoles interact in nontrivial ways with each other and the background potentials. Reduced ODE models for the dynamics of these vortex configurations can be found in [30,41,42,44] and are in good agreement with numerical simulations of the full equation (1.1) and with physical experiments, [30]. The objective of our work is to place these widely-studied models on a rigorous footing.
We are interested in the critical asymptotic regime where vortices interact with both the background potential p(x) and each other. In the following, we will take p(x) = p ǫ (x) with (1.2) |1 − p 2 ε (x)| 1 |log ε| and show that asymptotic scalings of p ε (x) outside of the range (1.2) induce dynamics that are dominated by only the background potential or by only vortex-vortex interactions.
1.1. Background potential. Following Lassoued-Mironescu [26] we write where η ε is a minimizer of an energy (A.2) and consequently a nontrivial solution to the following elliptic PDE, (1.3) 0 = ∆η ε + 1 ε 2 η ε (p 2 ε − η 2 ε ) in Ω ν · ∇η ε = 0 on ∂Ω , then (1.1) is equivalent to the Schrödinger equation The local and global well-posedness of (1.4) for fixed ε can be established using a variation of the argument in Brezis-Gallouet [7]. Associated to (1.4) is a weighted Hamiltonian, consequently, e 1 ε (w) is the standard gauge-less Ginzburg-Landau energy density as in [4]. We set to be the total energy of the mapping w associated to background potential η ε . We will be interested in an asymptotic range of p ε (x)'s such that the resulting energy feels both the background potential and the other vortices at the same order. As will be elicited below, our regime corresponds to having the background η ε expanded in the following way: There are a set of conservation laws for (1.4), and to write them down efficiently, we introduce some notation. For b, c ∈ C set (·, ·) = 1 2 bc + bc and let us define then solutions satisfy the following differential identities: A third conservation law holds for the Jacobian, we which we now define, tensorially. If the matrix J is given by J ij with J 21 = −J 12 = 1 and J 11 = J 22 = 0 then curl F = J kℓ ∂ ℓ F k and if J(w) = 1 2 curl j(w), then following Jerrard-Smets [16] we have after a lengthy calculation Multiplying by φ ∈ C 2 0 (Ω) yields a conservation law observed in [16] (1.8) The dynamics of a vortex can be inferred by choosing suitable test functions φ in (1.8) which elicit the vortex positions. Since we expect the Jacobian J(w) to be roughly equal to a delta function, localized at the site of each vortex α, we can choose φ = x m χ(x − α), m ∈ {1, 2} for smooth cutoff function χ(x) = 1 in B r and compactly supported in B 2r . Integrating the left hand side over the support yields, formally, πα. On the right hand side, the first term has support in an annulus about each vortex, due to the structure of the test function, and since the energy concentrates outside of the support of the annulus, one sees that it generates an O(1) vortex-vortex interaction term. The second term on the right hand side has support inside of a ball about each vortex, and since the energy tensor is large in this region [24], it generates an O(|log ε||∇ log η ε |) vortex-potential interaction term. Therefore, for a vortex to interact with both the background potential and other vortices, it requires that |∇ log η ε | = O(|log ε| −1 ) for both terms to be of the same order. 1 is negligible.
In fact if |∇ log η ε | ≪ |log ε| −1 then the vortex dynamics will be dominated by the interaction with the background potential, as in [16], whereas if |∇ log η ε | ≫ |log ε| −1 then vortex dynamics will be driven by its interaction with the current generated by fellow vortices, see [9].
1.3. Weak Topology. In order to measure the distance of the vortices to the site of the expected location of the vortices given by the ODE (1.11), we use the flat normẆ −1,1 (Ω). We denote the normẆ −1,1 unless this topology is used on a subdomain, such as We use a helpful estimate for concentrations in the flat norm, see Brezis-Coron-Lieb [6]. Define then the following holds: The third term on the right hand side of (1.8) Proof. For a proof see [16] for example.
Consequently, the flat norm provides a good way to measure the location of singularities. Since the Jacobian, J(w), converges to a sum of weighted delta functions, we will useẆ −1,1 topology as a convenient metric to track vortex trajectories.
Associated to a set of locations α j ∈ Ω and degrees d j ∈ {±1} we can define a canonical harmonic map, w * that describes the limiting harmonic map with prescribed singularities, where harmonic ψ * is chosen so that ∂ ν w * (α, d) = 0 on ∂Ω. We will use w * when there is no ambiguity.

Energy Expansion.
Given a collection of n vortices with centers α = {α 1 , . . . , α n }, the Hamiltonian E ηε ε (w) can be expanded out to second order with E ηε which will be established in the sense of Γ-convergence in Proposition 4.3. Here n will be used throughout to represent the number of vortices, Q 0 (x) is the limiting rescaled background perturbation, and W (α, d) is the renormalized energy that arises in the work of Bethuel-Brezis-Hélein [4] 2 , given by where Ω r (α) ≡ Ω\∪ j B r (α j ). In order to express the boundary terms, we define G(·, α, d) = n j=1 G(·, α j , d j ) where ∆G(·, α j , d j ) = 2πd j δ α j in Ω with G = 0 on ∂Ω.
1.5. Discussion. Rigorous results on vortex dynamics for Gross-Pitaevskii were established when p(x) ≡ 1 by Colliander-Jerrard [9] and also Lin-Xin [28] which showed that vortices satisfy the Kirchoff-Onsager ODE for Euler point vortices. When there is a recent result of Jerrard-Smets [16] that showed that vortices travel along level sets of the Thomas-Fermi profile. In the parabolic setting there were rigorous results by Lin [27] and Jerrard-Soner [18] when p(x) ≡ 1 and by Jian-Song [22] when p = O(1). Mixed dynamics with p = O(1) and further forcing terms were discussed by Serfaty-Tice [39]. A variational proof of the parabolic dynamics for p ≡ 1 was given by Sandier-Serfaty [35].
1.6. Results. We wish to describe the dynamics of vortices in (1.4). If the vortices are located at positions α j (t) with degree d j ∈ {−1, 1}, then the vortices move via the ODE, Given the ODE generated by (1.11), we can define the collision time, Then for any 0 ≤ T < T col we can define (1.13) r min := min which is the minimum vortex-vortex or vortex-boundary distance until time T . Clearly r min > 0. Our main result is the following that captures both vortex-vortex and vortex-potential interactions. An important hypothesis will be that the initial data has control of the excess energy (1.14) D ε (w, α) ≡ E ηε ε (w) − H ε (α, Q 0 ). We will shorten the notation to D ε when the w and α are readily apparent. Theorem 1. Let p 2 ε = 1 + ρε(x) |log ε| where the ρ ε satisfies hypotheses of Proposition 1.2 with k ≥ 4. Let {α 0 j , d j } is a configuration of vortices such that r α 0 > 0, d j ∈ {−1, 1} and suppose w 0 ε , 0 < ε < 1 satisfies the well-preparedness hypotheses If w ε (t) is a solution to (1.1) with initial data w 0 ε then there exists a time T > 0 such that for all t ∈ [0, T ], as ε → 0, where the α j (t) are defined by (1.11), and T is independent of ε.
1. An important function satisfying the conditions of Proposition 1.2 below is simply ρ ε = ρ 0 ∈ H k+1 (Ω) is any fixed function such that ∇ρ 0 has compact support inside Ω.
One important example that arises from the Theorem 1 are dipoles that are scattered by inhomogeneities, see Remark 1.2 below and Section 5 for some illustrative numerical simulations.
Under slightly less regularity assumptions, we have a similar result for the gradient flow |log ε| where the ρ ε satisfies hypotheses of Proposition 1.2 with k ≥ 3. If {α 0 j , d j } be a configuration of vortices such that d j ∈ {−1, 1} and if u 0 ε , 0 < ε < 1 satisfies (1.16) and w ε = uε ηε satisfies the well-preparedness hypotheses then there exists a time T > 0 such that for all t ∈ [0, T ], as ε → 0, where the α j (t) are defined by (1.17).
Remark 1.2. There are straightforward adaptations of Theorem 1 to other contexts, including: • The Gross-Pitaevskii equation (1.1) on Ω ≡ R 2 . This can be achieved by combining the methods here with the arguments in [5,21]. In this case vortices move according to the same ODE (1.11); however, the renormalized energy W (a, d) is the classical Coulomb potential, as studied in [23,31] which results in a modified ODE In Section 4 we present first and second order Gamma-convergence results for the energy E ηε ε (w ε ), see Proposition 4.1 and Proposition 4.3, respectively. These results are similar to those found in [1], albeit ported to the case with critical inhomogeneities.

1.7.
Properties of η ε = η ε (p ε ). We collect here some information required on the behavior of η ε , as related to the background fluctuations p ε . The discussion of the proof will be given in Appendix A and is based off of standard elliptic theory estimates, so we simply state the results that we require for the vortex dynamics here.
We write We will assume that ρ ε → ρ 0 (x) in H k (Ω) with k made precise below. Using results such as [10] for the Dirichlet problem, we expect that p ε and η ε should be quite close. We make precise the nature in which that is true in the following Proposition, which will be proved in the Appendix. We remark here that the elliptic analysis required for the convergence that will implement in the Appendix is somewhat non-trivial, as the overall ellipticity of the underlying nonlinear elliptic problem is going to 0 in the limit as ε → 0. Hence, estimates must be done with some care, especially to understand regularity up to the boundary of our domain.
In particular, we have Q ε → ρ 0 in C 2,δ (Ω) with δ > 0 for integer k ≥ 4 as required for the Schrödinger dynamics below and Q ε → ρ 0 in C 1,δ (Ω), δ > 0 for integer k ≥ 3 as required for the gradient flow dynamics. 3 for convenience. In fact, the sharp estimate would include Schauder theory estimates directly on Hölder norms, but that would require uniform convergence in C 1 , which does not seem obvious given that the ∇Q term appears to vanish in the limit preventing the proof of uniform bounds of higher regularity. Remark 1.3. In the appendix, we will actually decompose η ε = 1 + Qε | log ε| , in which case this convention slightly simplifies the resulting elliptic analysis.
Remark 1.4. Similar results hold assuming Dirichlet boundary conditions, but we work here with Neumann as they are the most physically relevant.

Excess Energy Control
In this section we present an excess energy identity similar to [9,20]. The influence of the background requires some modifications, so we give a full proof. The form of the estimate will be similar to that found in [5].
Remark 2.1. Note that (2.2) and Rellich-Kondrachov implies Proof. Although it is possible to get explicit control on the error o ε (1) in (2.1) by a careful analysis as in [20,25], the weaker estimate (2.1) is sufficient to prove the vortex motion law. 1. We first decompose the excess energy into where σ is appropriately chosen and σ ≤ r ≤ r α . We will control terms I − III by variants of estimates found in [20,25,16]. We will choose in the following. 2. By a simple variation of a standard calculation, see for example [20], and so our primary concern will be with estimating the integral of (2.8) on a suitable subdomain Ω σ (α). This last term can be rewritten as We now fix our µ 0 and ε 0 in order to be able to cite some estimates from [17,20]. The tool we wish to use is in the proof of Lemma 4 in [20], which yields an estimate given some assumptions on the Jacobian and ε. Let K 1 be the constant in the assumptions of Lemma 4 from [20], which depends only on Ω.
We follow an argument in the proof of Theorem 6.1 of [5] for fixing these constants; however, there are several new constraints including the fact that the vortex ball radius, σ, cannot be too small due to the need to control terms like . We further restrict µ 0 = min{µ 1 , rα 8K 2 n 5 } where K 2 comes from the assumptions in Theorem 3 of [20], needed for the localization result (2.4).
We now choose an ε 0 to allow us to use the necessary array of estimates. First, we choose ε 2 so that for all ε < ε 2 , We then set ε 1 ≤ ε 2 so that εE 1 [20]. In particular (2.10) and (2.11) imply that the assumptions in Lemma 4 hold, and we find If |log ε| − 2 3 ≤ µ then again the assumptions hold for Lemma 4 of [20], and choose s ε = µ and σ = √ µ so Taking the sum over the two errors yields sufficient control on Finally, we note Combining these estimates together controls term I. 3. We control II using the explicit control in Lemma 12 [20], therefore, we can invoke Theorem 1.3 and Lemma 6.8 of [17]: Therefore, and so (2.1) follows. 5. We next turn to the proof of (2.2) and (2.3). We first choose vortex balls B rα (α j ); by (2.9), and arguing as in Step 4, Therefore, Q 0 (α j ) + W (α) + C(n, r α , Q 0 ), (2.16) and so, Since r ≤ r α , we can combine (2.1), η 2 ε ≥ 1 2 , and (2.18) to get bound (2.2). To prove (2.3) we use an L p bound on each vortex ball B rα (α j ) in Theorem 3.2.1 of [9]. Outside the vortex balls, we use use the decomposition argument; in particular, since . From Remark 2.1 and the above bound, one finds ≤ C(r α , Q 0 , p), due to Remark 2.1. Combining the vortex ball bound with the bound in the excised domain yields (2.3). 5. The restriction on µ 0 allows us to use Theorem 3 of [20], and estimate (2.4) follows directly. To establish (2.5) we can follow the proof of Theorem 11 in [25], along with the equipartioning result of [24], which lets us localize the stress-energy tensor ∇w ⊗ ∇w up to an error of order O( |log ε|).

Proof of Theorem 1
The proof of the vortex motion law entails first establishing the smooth evolution of the vortex paths ξ j (t) via the compactness theorems in Section 2. The second step will compare the vortex paths with the ODE (1.11).
Proof. 1. We first claim that for any T > 0, there exists an ε such that J(w ε ) is a continuous function on [0, T ] into L 1 (Ω). Since J(u) = −∇ ⊥ Ru · ∇Iu then since the solution map is a continuous function into H 1 (Ω). 2. We use (1.8) along with the results in Proposition 2.1. Set Choosing 0 ≤ r ≤ s and using (1.8) and (1.15) we generate the following bound, and so We can now employ a diagonalization argument. Taking subsequence ε k ⊂ ε, we generate a dense, countable collection of times t ∈ [0, T ], in which T depends solely on r a , such that the ξ j (t) satisfy ξ j (0) = a 0 j with (3.1) as k → ∞. Since the collection is dense, we can take the limit as k → ∞. The Lipschitz extension of the ξ j (t)'s on [0, T ] satisfy the same conditions.

Vortex Dynamics.
In this subsection we complete the proof of Theorem 1. We use the Lipschitz paths ξ j (t) generated along the subsequence ε k from Proposition 3.1. We will examine the differences between vortex paths and the ODE, which will lead eventually to a Gronwall argument.
The proof requires two technical calculations. The first describes the gradient of W (α) in terms of the canonical harmonic map w * . Lemma 3.2. Let α ∈ Ω then w * = w * (·, α) and the renormalized energy W (α) satisfy where ϕ ∈ C 2 (Ω) and ∇ 2 ϕ has support in a neighborhood of the α j 's.
Proof. See for example [9,28] for a proof.
The second technical result proves weak compactness of the supercurrent away from the vortex cores. This will be used to bound certain cross terms in the Gronwall argument. The argument is along the lines of similar proofs in [9,28], among other places, and we sketch most of the argument.
We set φ to be the average of φ over Ω; and consequently from the Neumann boundary conditions on w ε , where we used a Hodge decomposition of ∂ t X that arises from differentiating φ and ψ in time. In order to control the term with h ε 2 , we use elliptic estimates and bound ∇h ε This implies t 0 +τ t 0 Ω X · J∇h ε 2 dxdt = o ε (1), and the three bounds combine together to yield (3.4).
With Lemmas 3.2 and 3.3 in hand, we can now turn to the proof of the main result.
Proof of Theorem 1. Define We will estimate the difference between the vortex paths ξ j , generated by (1.4), and the positions a j , extracted by the ODE (1.11). We also let D ε (w ε (t), b(t)) = E ηε ε (w ε (t)) − W ε (b(t)) denote the time-varying excess energy.
Recall from (1.11) that we will estimate the differences between vortex paths and the ODE positions: We can control A and B using our Lipschitz bounds; in particular, where we used Lemma 10 of [20] to control the Hessian of W (a), with a bound dependent on r a . Likewise, we can use Appendix A on the regularity of η ε k to bound b a |∇Q(ξ j (t)) − ∇Q(α j (t))|dt ≤ Q 0 W 2,∞ b a µ(t)dt. (3.9) We now turn to the control on C. For shorthand we define components j k = ( j(wε k ) |wε k | ) k and j * k = (j(w * (ξ))) k of the supercurrents, and since ∂ k w ε k = (∂ k |w ε k | + ij k ) , then from Lemma 3.3, we have that lim k→∞ C 3 + C 4 = 0; note this is where we needed to consider the time-integral on (a, b).
Next we control D.
We can now bound and Finally, we bound E using (2.2): Combining together estimates generates a differential inequality, and choosing b = a + δ and dividing by δ implies the classical Gronwall inequality with µ(0) = 0. Therefore, µ(t) = 0 for all 0 ≤ t ≤ T . We can then restart the problem at time T and continue the dynamics until r a → 0. So far we have only shown the theorem for a subsequence. However, our argument can be applied to subsequences of any sequence of ε k → 0, and as the limit is independent of the chosen subsequence, we obtain convergence for ε → 0 without taking subsequences.

Second order Gamma Convergence and Gradient Flow
In this section we collect and restate some of our results on the static energy in the spirit of Γ-convergence, compare [1] for the corresponding results for p ε ≡ 1. For dynamics, we show that the approach of [35] can be adapted to treat the gradient flow as opposed to the Schrödinger dynamics studied in the previous sections.
All of our arguments work with either Neumann or Dirichlet boundary conditions, just as those in [1] and [35]. For simplicity, we assume Neumann boundary conditions throughout, and sometimes comment on changes necessary for the Dirichlet case.
Proposition 4.1. Let w ε ∈ H 1 (Ω; C) be a sequence of functions with E ηε ε (w ε ) ≤ K 1 |log ε|. Then the Jacobians J(w ε ) are compact inẆ −1,1 (Ω) and for a subsequence, where a i ∈ Ω are distinct points and d i ∈ Z \ {0}. Furthermore, Proof. Comparing with the standard Ginzburg-Landau energy E 1 ε , we clearly have (4.2) min As both the minimum and the maximum can be estimated as 1 + O( 1 |log ε| ), we obtain that the standard Ginzburg-Landau energy is bounded as This implies the compactness of J(w ε ), see [19,37].
The lower bound for the energy follows from the standard Ginzburg-Landau lower bounds and (4.2).
The following bound is the heart of the second order Γ-convergence argument.  Furthermore, (w ε ) are then bounded in W 1,p (Ω) for all p ∈ [1, 2) and satisfy the energy lower bound where W is the renormalized energy defined in (1.10). Also, there is K 3 such that Proof. We clearly have E 1 ε (w ε ) ≤ J * |log ε| + K 4 with K 4 depending on K 2 and J * as well as on Q 0 . Hence we can use Ginzburg-Landau arguments of Colliander-Jerrard [9] or Alicandro-Ponsiglione [1] to obtain the claimed structure of J * . The L p bound follows from Proposition 2.1. Let now r > 0 be so small that B r (a i ) are disjoint and do not intersect ∂Ω. Then we apply Proposition 4.2 to obtain In Ω r (α), we let w r be the optimal S 1 valued map as in the definition of the renormalized energy W , and letŵ * be the subsequential weak limit of w ε . Then |ŵ * | = 1 a.e. and hence by lower semicontinuity of the Dirichlet integral, By (4.2), we deduce that also lim inf ε→0 E ηε ε (w ε ; Ω r (α)) ≥ 1 2 Ωr(α) |∇w r | 2 dx.
Using that we obtain the claimed lower bound. The upper bound for the penalty term has been shown in Proposition 2.1.

4.2.
The gradient flow. Now we study the gradient flow analog of (1.1), Instead of deriving an equation of motion by studying the limits of the energy density, we are going to use the method of Γ-convergence of gradient flows. This has the advantage that we will require slightly less regularity and convergence properties for the potential term p(x) than in the Schrödinger part of this article. We will follow closely the approach of Sandier-Serfaty [35] to obtain the convergence of the gradient flow of E ηε ε to that of F . In fact, the proof goes through with almost no essential changes, and we will therefore assume the reader is familiar with the argument and some of the notation of [35] and mostly highlight the differences.
Again, we use the approach of dividing by the profile η ε , leading to an equation for w ε = uε ηε . As we assume Neumann boundary conditions ∂ ν η ε = 0 for η ε , then Neumann boundary conditions ∂ ν u = 0 imply Neumann boundary conditions ∂ ν w ε = 0. As before, we obtain in Ω ∂wε ∂ν = 0 on ∂Ω This is the gradient flow of E ηε ε with respect to the scalar product on X ε = L 2 (Ω; C) given by the quadratic form In fact, we have for φ ∈ C ∞ (Ω; C) the Gâteaux derivative w ε · φdx and on the other hand, ∂ t w ε , φ Xε = 1 |log ε| Ω η 2 ε ∂ t w ε φdx so integrating by parts we see that (4.4) is a gradient flow, more precisely ∂ t w ε = −∇ Xε E ηε ε (w). We will show that we can relate this gradient flow to that of F (α), where for α = (α 1 , . . . , α n ) ∈ Ω n , α i = α j for i = j, we set Q 0 (α j ) and we will use the following metric on Y = R 2n (the tangential space to Ω n ): Theorem 3. Assume that the profiles η ε solving (1. Let (w ε ) be a family of solutions of (4.4) with initial data w 0 ε . Assume that w 0 ε satisfy J(w 0 ε ) → π d j δ α 0 j inẆ −1,1 with d j ∈ {±1}. Let α(t) : [0, T m ) → R 2n be a solution of the system of ODEs where T m = T (α 0 ) ∈ (0, ∞] is the maximal time of existence for (4.5), i.e. the smallest time T such that as as t → T , we have r α(t) → 0.
If the initial data are very well-prepared, i.e. D ε (0) = E ηε ε (w 0 ε ) − H ε (α 0 ) → 0 as ε → 0, then for 0 ≤ t < T m , J(w ε (t)) → π d j δ α j (t) so we have convergence of the gradient flow of E ηε ε to that of F . Furthermore, D ε (t) = E ηε ε (w ε (t)) − H ε (α(t)) → 0 so the data continue to be well-prepared. Finally, we have for T < T m and any B i (t) that are disjoint open balls contained in Ω and centered at α i (t) that Remark 4.1. See Proposition 1.2 for conditions on the potential p ε that imply the convergence of Q ε assumed in the theorem. Note that we require less regularity and convergence than for the Schrödinger case.
We can now apply Proposition 3.2 and Proposition 3.3 of [35] to w ε and obtain that (for a subsequence) J(w ε (t)) → π d j δ b j (t)) , where b j ∈ H 1 (0, T 0 ; R 2 ) with b j (0) = α 0 j , i.e. the d j are constant. Corollary 7 of [36] now applies to show lim inf ε→0 1 |log ε| As |η 2 ε − 1| ≤ C |log ε| , this also implies The last equation is the lower bound part needed for the Γ-convergence of gradient flows argument.
To construct this path in function space, we use the same pushing as in [35]. Let B i = B ρ (α i ) with ρ < r α be pairwise disjoint balls and define χ t : Ω → Ω to be a oneparameter family of diffeomorphisms that satisfies With the phase corrector function ψ t defined as in (3.24) of [35], set v ε (χ t (x), t) = w ε (x)e iψt(x) .
To calculate the energy change along this path, we change variables y = χ t (x) and obtain with Jac(χ t ) denoting the Jacobian determinant, We now differentiate in time and set t = 0. When the derivative does not hit η ε • χ t , the resulting terms can be dealt with exactly as in [35] by our observation (4.7). The remaining terms are . Using the uniform convergence 1 |log ε| ∇(η 2 ε ) → ∇Q 0 and the fact that the logarithmic part of the energy is concentrated in the B i (recall Proposition 2.1 or (3.17) of [35]), we deduce that A 1 → π V i · ∇Q 0 (α i ) as ε → 0. To show that A 2 → 0, we recall (2.2) and note that |∇(η 4 ε )| → 0 uniformly in Ω. Together with the results of [35] for the terms not involving derivatives of η ε , we deduce (4.10).
On the interval [0, T 0 ), we can now apply Theorem 1.4 of [35] and obtain the claim of Theorem 3 up to the time T 0 . The global statement up to T m follows as in Section 3.3 of [35].

Numerical Simulations of the Vortex Dynamics
Armed with the dynamical equation (1.11), we can simulate the flow of vortices and observe the impact of the background potential at the scales we have studied in Theorem 1. For our simulations, we have used the renormalized energy which is the correct expression for Ω = R 2 , see Remark 1.2. This is an approximation for the renormalized energy for Ω = B R (0) for R ≫ 1 sufficiently large compared to |α|. We will focus here on the case of the dipole (a pair of vortices of opposite charge) interacting with potentials Q 0 = V i of the form Step Function to different material background, (5.2) While these potentials are not compactly supported, we can apply Theorem 1 by truncating the potentials outside a suitably large domain without affecting the dynamics we are plotting.
Recall that in the absence of the background potential, dipole dynamics simply move in in a straight line perpendicular to that connecting the vortex centers at a speed correlated to the vortex spacing.
The equations (1.11) for each choice of background are then plugged into the ode15s ODE solver in Matlab and integrated over time scales long enough to observe the impact of the background potential on the dipole dynamics. The results are recorded graphically in Figure 1.
Since we are working on a bounded set Ω, such a nontrivial η ε exists in H 1 for ε sufficiently small and p ε sufficiently bounded in H 1 using a direct method and Rellich-Kondrachov. We can improve the regularity away from the boundary. In particular ifρ ε ∈ H k then η ε ∈ H k+2 , even though the H k+2 norm blows up as ε → 0. This follows from standard elliptic theory results involving nonlinear Sobolev embeddings to bootstrap regularity, see for instance [8], Section 2.4. The existence of Euler-Lagrange equations, strong solutions, and classical solutions follow. As noted in (1.19) this is a slight shift of notation, as the Q 0 defined in the main text is not quite the limit of this family of functions, rather it takes the form The Euler-Lagrange perturbation equation is written as Then we can write this as 2. We next establish ε-dependent estimates on a model problem.
The required estimates follow by, for instance, carefully modifying [13][Theorem 7.32], which draws upon ideas originally put forth by Nirenberg [33]. The proof relies upon an analysis of regularity up to the boundary by flattening locally to the half-plane and using the Green's function there, application of difference operators to gain regularity and induction on regularity estimates for elliptic problems with Neumann boundary problems. Essentially, the argument boils down to the fact that the operator is coercive for each ε however. The equation also has a unique set of solutions by the same property. 4 To establish the bounds we need more precisely, we follow [13][Theorems 7.32 and 7.29], where the author studies Sobolev estimates on coercive elliptic equations. In our setting, these equations take the form L ε w := 1 + 1 |log ε| A − ε 2 2 ∆ w = f, ∂ ν w = 0 on Ω smooth enough, for A a bounded function of the same regularity as p.
We will first consider this equation on half-balls of radius s with flat boundary, say N(s) where y parametrizes the boundary and x parametrizes the interior. We claim we have an 4 Perhaps the closest argument of this type for Neumann boundary conditions can be found in [43], Chapter 5.7, Propositions 7.4 and 7.5 where it states that for a smooth enough domain Ω ⊂ R 2 , there exists a unique solution v to the elliptic equation and for all k = 0, 1, 2, . . . , given f ∈ H k , we have v 2 H k+2 (Ω) ≤ C ∆v 2 H k (Ω) + C v 2 H k+1 (Ω) . Since we need precise control on the constants of our estimate with respect to ε for a perturbation of this equation, which are not available by rescaling, we have included a proof for completeness. estimate of the form To see this, we first claim that (A.6) ε 2 ∂ γ x w H 1 (N (r)) ≤ f H k (N (r)) withr < r, with γ ≤ k + 1. For γ = 0, this is exactly the coercivity estimate. Otherwise, we use an inductive argument constructed via difference operators inside the Dirichlet form D(∆ γ h ζw, ∆ γ h ζw) for ζ a cut-off to N that vanishes on the curved part of the boundary of N. Commuting with the cut-off function and integrating by parts when necessary, we observe where then by induction we have Similarly, we claim that (A.7) ∂ γ x w L 2 (N (r)) ≤ f H k (N (r)) for 0 ≤ γ < k. This follows by instead putting all the derivatives onto f and using L 2 instead of H 1 norms for the f term on the right hand side. At the boundary, we can then establish (A.5), again proven using induction. To see this, we recognize that if γ is a multi-index and γ 2 = 0 or 1, then the estimate follows by (A.6) and (A.7). For use the fact that If the number of derivatives on y is 0 or 1, we use the above estimate. For more than that, we pull off the first two derivatives in y, make the substitution, then use the inductive hypothesis.
Then, once such estimates are established on half-balls, we can create a sequence of cutoff functions U 0 = B(0, R), . . . , U k = V 0 , U j+1 ⊂ U j , where V 0 ∪ V 1 ∪ · · · ∪ V M = B(0, R) and V j can be mapped to a half-ball for j > 0 and W 0 is an interior region. would require further work controlling the error terms relating to boundary condition of ρ ε and potentially restrict us to local convergence estimates. 4. To finish the a posteriori convergence, we need to get to H k convergence. Recall the equation (1 − ε 2 ∆)q = g ∂ ν q = 0, which we rewrite as (A.10) (q − g) − ε 2 ∆q = 0 ∂ ν q = 0.
We will replace q with Q ε and g by the right hand side of our elliptic equation.
For higher regularity, we have the formal calculation, where the boundary terms can be controlled by careful use of the Neumann boundary condition and the equations.
To establish the result rigorously up to the boundary, we must repeat the argument as in Step 2, but acting ∆ γ h on both sides of (A.10) and then multiplying both sides by ∆ γ h (q − g) and integrating. Specifically, we can consider the coercive Dirichlet form both in neighborhoods of the boundary, as well as in the interior to get elliptic estimates as in Step 2.