Small-data global existence of solutions for the Pitaevskii model of superfluidity

We investigate a micro-scale model of superfluidity derived by Pitaevskii (1959 Sov. Phys. JETP 8 282–7) to describe the interacting dynamics between the superfluid and normal fluid phases of Helium-4. The model involves the nonlinear Schrödinger equation (NLS) and the Navier–Stokes equations, coupled to each other via a bidirectional nonlinear relaxation mechanism. Depending on the nature of the nonlinearity in the NLS, we prove global/almost global existence of solutions to this system in T2 —strong in wavefunction and velocity, and weak in density.


Introduction
Superfluids constitute a phase of matter that is achieved when certain substances are isobarically cooled, resulting in Bose-Einstein condensation.That Helium-4 (and also its isotope Helium-3) undergoes such a quantum mechanical phase transition was first experimentally discovered [Kap38,AM38] over 80 years ago and has been the subject of intense inquiry ever since.Despite this, a single theory that describes the phenomenon continues to elude us.
The general picture is that at non-zero temperatures, there is a mixture of two interacting phases: the normal fluid and the superfluid [PL11,Vin04,Vin06,BSS14,BDV01,BLR14].It is important to note that this is not like classical multiphase flow, where one can define a clear boundary between the two phases.Instead, some atoms are in the normal fluid phase, and some are in the superfluid phase, with both fluids occupying the entire volume.The normal fluid is well-modeled by the Navier-Stokes equations (NSE), while the description of the superfluid varies by the length scale that we are interested in (see [BBP14,Jay22] for a discussion).Briefly, the superfluid is described by the NSE at large scales [Hol01], a vortex model at intermediate scales [Sch78,Sch85,Sch88], and the nonlinear Schrödinger equation (NLS) at small scales [Kha69,Car96].The macro-scale, NSE-based description is a current topic of numerical research [VSBP19, RBL09, SRL11], and has also been rigorously analyzed [JT21].In this paper, we use the micro-scale, NLS-based model by Pitaevskii [Pit59], which has previously been considered in [JT22a,JT22b].
A missing piece of the physics puzzle here is the nature of the interaction mechanism.It is known that the interaction between the fluids is dissipative/retarding.Pitaevskii thus derived a micro-scale model that intertwines the NLS (for the superfluid) and the NSE (for the normal fluid).The coupling is nonlinear, bidirectional and transfers mass, momentum, and energy between the two fluids.For the combined system of both phases, the model respects the conservation of total mass and total momentum, while the total energy decreases in accordance with the dissipation.
The NLS, in its most popular form, is fundamentally a dispersive partial differential equation with a cubic nonlinearity that models systems with low-energy wave interactions, such as dipolar quantum gases [CMS08,Soh11].The well-posedness issues of NLS have been tackled in many situations [CKS + ], and its scattering solutions [Tao06,Dod16] have been of particular interest.The NLS can also be recast as a system of compressible Euler equations (referred to as quantum hydrodynamics or QHD) with an additional quantum pressure term [CDS12].This system is a special case of the more general Korteweg models, subject to much mathematical analysis.Hattori and Li [HL94] showed that the 2D QHD equations are locally well-posed for high-regularity data, and improved this to global well-posedness in the case of small data [HL96].Jüngel [JMR02] established local strong solutions to the QHD-Poisson system, formed by including a potential governed by the Poisson equation.The same model possesses local-in-time classical solutions in 1D when the data is highly regular [JL04].For initial conditions close to a stationary state, the solutions are globalin-time and converge exponentially fast to the stationary state.Blow-up criteria have also been derived for QHD [WG20,WG21].While the discussion so far has focused on strong solutions, there has also been rising interest in the weak formulation of QHD-like models.Antonelli and Marcati [AM09, AM12, AM15] introduced the novel fractional step method in the pursuit of finite-energy global weak solutions.The idea was to revert (from QHD) to NLS, which was easier to solve, and account for collision-induced momentum transfer via periodic updates to the wavefunction.In this process, the occurrence of quantum vortices could also be characterized by imposing irrotationality of the velocity field (away from vacuum regions).Using special test functions that permit better control of the quantum pressure term, Jüngel [Jün10] proved that the viscous QHD system admits weak solutions in 2D.For small values of viscosity, these solutions were global in time.The proof utilized a redefinition of the velocity that converts the hyperbolic continuity equation into a parabolic one, a technique that was pioneered by Bresch and Desjardins [BD04] for Korteweg systems in general.Vasseur and Yu [VY16b] expanded Jüngel's result to a wider class of test functions while adding some physically-motivated drag terms.Various forms of damping have appeared in the literature, primarily serving two different roles: (i) as an approximating scheme for both the compressible Navier-Stokes with degenerate viscosities [LX15,VY16a] as well as Korteweg-type systems [AS17, ACLS20, AS22], and (ii) as a means of proving global existence [Cha20] or relaxation to a steady state [BGLVV22,SYZ22].Most works involving Korteweg systems use the notion of κ-entropy that was first demonstrated in [BDZ15].Furthermore, even questions of non-uniqueness (and weak-strong uniqueness) of weak solutions have been addressed for the QHD-Poisson system with linear drag using convex integration [DFM15].
It is only at absolute zero temperature that superfluids can be well-approximated by the use of the NLS alone.For temperatures above zero and below about 2.17K, we have a mixture of both fluids.In this article, we consider Pitaevskii's model [Pit59] which couples the NLS and the NSE.The model was initially derived for a fully compressible normal fluid.While compressible fluids are more realistic in some scenarios, they are also much more challenging to both rigorously analyze and numerically simulate.[Fei04,Lio96a] contain several classical results on the compressible NSE.On the other hand, the incompressible NSE (no density equation) is arguably the most studied nonlinear partial differential equation in mathematics (see [Tem77,MB02,RRS16] for classical results).In this article, we approximate the normal fluid to be incompressible, but the density persists, varying from point to point in the flow domain.What results is an incompressible, inhomogeneous flow: compressible NSE appended with the condition of divergence-free velocity.This model of fluids was first investigated by Kazhikov for local weak solutions when the initial density is bounded from below [Kaz74], and vacuum states were allowed in an improvement by Kim [Kim87].Further advances for weak solutions were made by Simon [Sim90], who in particular analyzed their continuity at t = 0, and also proved the existence of global solutions in a less regular space.Meanwhile, Ladyzhenskaya and Solonnikov [LS78] presented the case for strong solutions: With the density bounded from below, it is possible to construct local (global) unique solutions in 3D (2D).Furthermore, if the data is small enough, one obtains global-in-time unique solutions.Results in the same spirit were proven by Danchin for small perturbations from the stationary state in critical Besov spaces [Dan03].He further established the inviscid limit of the incompressible inhomogeneous NSE in subcritical spaces [Dan06].The local existence theorem by Ladyzhenskaya and Solonnikov was shown to be valid for non-negative densities as long as the initial data satisfied a compatibility criterion [CK03].This work by Choe and Kim has since spurred on several other results that utilize such compatibility conditions on the initial data.
Given the immense interest in the NLS and NSE, the rigorous study of a coupled system should be a natural next step.Indeed, one such two-fluid model of superfluidity was analyzed by Antonelli and Marcati in [AM15].The superfluid was described by the NLS, and the normal fluid by the compressible NSE.This is similar to the system considered in this article, save for two key differences.Firstly, their model did not permit any mass transfer between the two fluids (which allows for global-in-time solutions).As we shall discuss, this is the biggest roadblock in Pitaevskii's model and essentially defines the strategy used.Secondly, the momentum transfer in their model is unidirectional and linear, affecting only the superfluid phase (as opposed to the bidirectional and nonlinear nature of the coupling in this work).
Thanks to the retarding interactions between the two phases, the NLS acquires a dissipative flavor and renders it parabolic.This lets us extract dissipative contributions to the energy estimates.To analyze the momentum equation of the NSE, we work with initial velocity in H 1 d .This yields appropriate regularity for the velocity, in order to adequately control the relaxation mechanism which contains quadratic terms in the velocity.Parting ways from [Kim87], we begin with an initial density field that is bounded from below.This is necessary since the continuity equation is unusual and is not a homogeneous transport equation.Our primary goal is to avoid the occurrence of zero or negative densities at any time.To this end, we must limit the effect of inhomogeneity, which is the relaxation mechanism that allows for mass and momentum transfer between the two fluids.As a serendipitous by-product of this non-zero density field, we also obtain control of x , which allows the use of compactness arguments to actually obtain strong continuity in time of the velocity field.
The crux of this work is to derive a priori estimates and carefully extract coercive terms that allow for norms to decay, while avoiding any derivatives on the density of the normal fluid.To engineer this decay, we include a linear drag term for the NSE.Additionally, we also present results for any polynomial-type nonlinearity in the NLS.We now mention the notation used in the article before describing the model and stating the results.
1.1.Notation.We denote by H s (T 2 ) the completion of C ∞ (T 2 ) under the Sobolev norm H s , while we use Ḣs (T 2 ) when referring to the homogeneous Sobolev spaces.Consider a 2D vector-valued function u ≡ (u 1 , u 2 ), where ) under the H s norm.The L 2 inner product, denoted by ⟨•, •⟩, is sesquilinear (the first argument is complex conjugated, indicated by an overbar) to accommodate the complex nature of the Schrödinger equation, i.e., ⟨ψ, φ⟩ := T 2 ψφ dx.Since the velocity and density are real-valued functions, we ignore the complex conjugation when they constitute the first argument of the inner product.
We use the subscript x to denote Banach spaces that are defined over T 2 .For instance, L p x := L p (T 2 ) and H s d,x := H s d (T 2 ).For spaces/norms over time, the subscript t denotes the time interval in consideration, such as L p t := L p [0,t] .The Bochner spaces L p (0, T ; X) and C([0, T ]; X) have their usual meanings, as L p and continuous maps (respectively) from [0, T ] to a Banach space X.
We also use the notation X ≲ Y and X ≳ Y to imply that there exists a positive constant C such that X ≤ CY and CX ≥ Y , respectively.When appropriate, the dependence of the constant on various parameters shall be denoted using a subscript as the article, C is used to denote a (possibly large) constant that depends on the system parameters listed in (2.4), while κ and ε are used to represent (small) positive numbers.The values of C, κ, and ε can vary across the different steps of calculations.
1.2.Organization of the paper.In Section 2, we present and discuss the mathematical model, along with statements of the main results.Several a priori estimates, at increasing levels of regularity, are derived in Section 3. The construction of the semi-Galerkin scheme and the renormalization of the density are discussed in Section 4.

Mathematical model and main results
The superfluid phase is described by a complex wavefunction, whose dynamics are governed by the nonlinear Schrödinger equation (NLS), while the normal fluid is modeled using the compressible Navier-Stokes equations (NSE).In all generality, the full set of equations can be found in [Pit59,Section 2].In what follows, we use a slightly simplified and modified version of the equations, arrived at by making the following assumptions.
(1) We consider a general power-law nonlinearity for the NLS.This is done by choosing the internal energy density of the system to be 2µ p+2 |ψ| p+2 , for 1 ≤ p < ∞ (see Remark 2.5).We also assume that the internal energy is independent of the density of the normal fluid.
(2) We work in the limit of a divergence-free normal fluid velocity.This means that the pressure is a Lagrange multiplier, rendering the equations of state and entropy unnecessary.Note that, due to the nature of the coupling between the two phases, the density of the normal fluid is not simply transported.(3) A linear drag term has been included in the momentum equation to account for the lack of coercive estimates for the velocity.(4) Planck's constant (h) and the mass of the Helium atom (m) have both been set to unity for simplicity.
We now state the equations used in this paper: Here, ψ is the wavefunction describing the superfluid phase, while ρ, u, and q are the density, velocity and pressure (respectively) of the normal fluid.The normal fluid has viscosity ν and drag coefficient α, while µ (positive constant) is the strength of the scattering interactions within the superfluid 1 .This scattering nonlinearity has an exponent p ∈ [1, ∞).Finally, λ is a positive constant that indicates the coupling strength between the two phases.The coupling is denoted by the nonlinear operator B.
The Schrödinger equation dictates the evolution of the wavefunction, generated via the action of the Hamiltonian (roughly, the energy) of the system.The coupling B resembles the relative kinetic energy 2 between the two phases.This is evident upon recalling that the quantum mechanical momentum operator (in the position basis) is −ih∇.The purpose of this coupling is to allow for mass/momentum transfer between the two phases as a means of relaxation or dissipation.
These equations are supplemented with the initial conditions We use periodic boundary conditions, i.e., we are working on the two-dimensional torus [0, 1] 2 .
2.1.Weak solutions and the existence theorems.Having stated the model, the notion of weak solutions to (NLS), (NSE), (CON), and (DIV) (with initial conditions (INI) and periodic boundary conditions), henceforth referred to as the Pitaevskii model, is as follows.
Definition 2.1 (Weak solutions 3 ).For a given time T > 0, a triplet (ψ, u, ρ) is called a weak solution to the Pitaevskii model if the following conditions hold.
2 There is also the nonlinear wavefunction term, so that the relaxation to equilibrium also depends on the potential energy of the superfluid.
Remark 2.2.We note that the last two terms in (NSE) are gradients, just like the pressure term, and thus vanish in the definition of the weak solution (since the test function is divergence-free).Henceforth, we absorb these two gradient terms into the pressure, relabeling the new pressure as q.
We are now ready to state our main results.
Theorem 2.3 (Global existence).Fix any p ∈ [1, 4), and let e. in T 2 .Then, there exist a global weak solution (ψ, u, ρ) to the Pitaevskii model such that the density is bounded between m f ∈ (0, m i ) and M f := M i + m i − m f , if the initial data satisfy the smallness criterion Also, the solution has the regularity for 1 ≤ r < ∞.Additionally, the solution also satisfies the energy equality 1 2 ρ(t)u(t) (2.8) For the case of higher-order nonlinearities, i.e., when p ≥ 4, we obtain "almost global" existence.
Theorem 2.4 (Almost global existence).In the case of p = 4, the solution to the Pitaevskii model has the same regularity properties as in Theorem 2.3, except that their existence is guaranteed on , where ε is the size of the (sufficiently small) initial data.
For p > 4, the existence time scales polynomially with the size of the data, as T ∼ ε − p p−4 .In both cases, these solutions also satisfy the energy equality on [0, T ].
While deriving the a priori estimates, we have to distinguish between the cases 1 ≤ p < 2, p = 2, 2 < p < 4, p = 4, and p > 4.This is due to the poor control we have on the superfluid mass.Given that we are on T 2 , and our equations do not preserve functions with vanishing mean, the L 2 norm becomes the limiting factor even in the decay of higher norms.In the case of the wavefunction, this corresponds to the mass of the superfluid.Similarly, for the velocity, we do not get coercive estimates from the viscosity term alone, at least at the level of the kinetic energy estimate.Thus, we introduce a linear drag term.
Remark 2.5.Since the self-interaction term in (NLS) involves a discontinuity due to the complex magnitude, evaluating the H 2 norm as in (3.51) requires p ≥ 1.In particular, points of superfluid vacuum (ψ = 0) may lead to problems.As an illustration, consider D 2 (|f | p f ) for a real-valued function f , which can be regularized as D 2 (f 2 + ε) p 2 f .Upon differentiation, the most problematic term is (f 2 + ε) p 2 −2 f 3 (Df ) 2 .To be able to handle this term in the limit ε → 0, at the points where f = 0, we require that 2 p 2 − 2 + 3 = p − 1 ≥ 0. This argument can be easily extended to a complex-valued function.
Remark 2.6.The regularity of the solutions seem to suggest that the wavefunction and velocity are strong solutions.Indeed this is true, as they are strongly continuous in their topologies.On the other hand, the density is truly a weak solution and is the reason for referring to the triplet as a weak solution.This low regularity of the density influences the nature of the calculations that are employed.
The proofs of both Theorems 2.3 and 2.4 follow from detailed a priori estimates, and a semi-Galerkin scheme to construct the solutions.The a priori estimates only differ slightly for various ranges of the values of p, as will be illustrated.The general approach to the problem is motivated by that of [Kim87], but we do not allow the density to vanish anywhere.This is because the presence of u in the nonlinear coupling means we are required to control it in L ∞ (T 2 ) to prevent the formation of vacuum (and regions of negative density).Beginning from the usual mass and energy estimates, we derive a hierarchy of several energies for the wavefunction and velocity.
2.2.Significance of the results.The holy grail of superfluid modeling is to find a unified description that works at all length scales, and rigorous validation of any proposed models is crucial to this process.The thrust of this paper is the analysis of Pitaevskii's description of superfluidity, the most important feature of which is to characterize the mass transfer between the two fluids.In the course of proving the main theorems, we quantify the conversion of superfluid into normal fluid (Lemma 3.1), confirming the interaction-induced relaxation mechanism.We establish the validity of the model in the limit t → ∞ even as the superfluid mass decreases (polynomially) quickly.The transition in the behavior of the solutions, from global to almost-global, as the self-interactions are increased in strength, is in accordance with the decreasing mass decay.However, the threshold p = 4 still begs for a physical explanation.Of the assumptions underlying our theorems, relaxing the demands of small data and positive normal fluid density would be important future advancements in the context of the Pitaevskii model.
The rigorous analysis of superfluid models is a fairly new topic, and we expect for this work to pave the way for further results in this direction.Some questions of interest, particularly of consequence to physicists and engineers, are the issues of stability and compressibility.For example, in [Pit59], Pitaevskii investigated the propagation of sound waves in superfluid Helium by studying the case when the superfluid has only small density gradients.It has to be noted that his derivation of the model accounted for the contributions to the internal energy of the system from both fluids.Thus, by utilizing appropriate self-interactions (for instance, non-local potentials, or including the normal fluid density), it would be important to test the model against experimental findings.A mathematical guarantee of the existence of solutions to the Pitaevskii model is essential to complement the efforts to numerically simulate such complicated systems [BSZ + 23].It is worth mentioning that a better understanding of superfluidity could be revolutionary to most modern experiments in physics (including the Large Hadron Collider [Leb94,RM18]), and also to the fields of quantum computing [HDT21], gravitational wave astronomy [SDLPS17], and dark matter [vKEE + 23].All of these use helium as a cryogen, often as a superfluid-normal fluid mixture due to the superfluid's excellent thermal conductivity [Vin04].
2.3.The strategy.The nonlinear coupling terms in (NLS) and (NSE) may be the most obvious differences between this model and other standard fluid dynamics models, but the source term in (CON) is the most troublesome.The backbone of our approach towards proving global existence is ensuring a positive lower bound for the density at all times.This involves a meticulous handling of the a priori estimates so as to obtain coercive terms that lead to global-in-time bounds.Throughout the calculations, we ensure that the density norms are only in Lebesgue spaces: ρ is not smooth enough to be differentiated (even weakly).Before we outline the strategy, we discuss some properties of the coupling operator B. Henceforth, we refer to the linear (in ψ) part of B as B L .Thus, (2.9) Proof.Both calculations follow using integration by parts.
(1) By (2.9) and incompressibility of u, we have (2) Similarly, In the last inequality, we used Hölder's and Young's inequalities to cancel the third term with the first two terms: x .□ Remark 2.8.Given that B provides a relaxation mechanism, it is tempting to treat it, or at least its linear part B L , as a dissipative second-order elliptic operator whose eigenfunctions can be used as a basis for the semi-Galerkin scheme.Even though B L is symmetric and has a non-negative real part, this cannot work since it has time-dependent coefficients, and so its eigenvalues and eigenfunctions also depend on time.Moreover, B L does not have a spectral gap at 0. Its eigenvalues are not known to be bounded from below by a positive number.
Thus, by integrating (CON) over T 2 , the advective term vanishes and using Lemma 2.7, we have This implies that the overall mass of the normal fluid does not decrease with time.Put differently, the coupling causes superfluid to be converted into normal fluid, on average.However, the RHS of (CON) need not be non-negative pointwise in T 2 .So it is not inconceivable that the density of the normal fluid may locally vanish, or even take negative values!To prevent physically unrealistic density fields, and because our estimates require a strictly positive density, we fix a positive lower bound for ρ.Based on this, we define our existence time T , so that ρ does not drop below the lower bound until time T .Our goal is to show that this lower bound can be maintained for arbitrarily long, provided we begin from sufficiently small data.
Definition 2.9 (Existence time).Start with an initial density field 0 Given 0 < m f < m i , we define the existence time for the solution as (2.11) A formal solution to the continuity equation can be written using the method of characteristics.Let X α (t) be the characteristic starting at α ∈ T 2 .To wit, the characteristic solves the differential equation (2.12) Here, u is the velocity of the normal fluid.So, along such characteristics, (2.13) From (2.11) and (2.13), it is clear that a sufficient condition to ensure the density is bounded from below by m f is 2λ for all T ≤ T * .This can in turn be ensured through the sufficiency So, we are looking to show that (2.15) − actually a stronger version of it − holds irrespective of T , so that we can conclude that the density is always greater than m f .This is achieved by selecting small enough data, and allows us to deduce the global existence of solutions.Since Bψ involves a second-order derivative, its L ∞ x boundedness leads us to high-regularity spaces.The momentum equation (NSE) is used to estimate ∥u∥ L 2 t H 2 x and ∥u∥ L 2 t H 1 x , which are useful in handling parts of ∥Bψ∥ L ∞ x .As a by-product of these calculations, we are also able to bound x , which plays a part in the compactness arguments for the strong time-continuity of u.The Schrödinger equation (NLS) is used to derive increasingly higher-order a priori estimates of ψ.In all these calculations, we work with density that is only in L ∞ x .

A priori estimates
Throughout this section, we derive the required a priori estimates, using formal calculations.We assume the wavefunction and velocity are smooth functions and that the density is bounded from below by m f > 0 in [0, T ].Here, T is any time less than the local existence time T * , and is extended to global existence in Section 3.5.Proof.Multiplying (NLS) by ψ, taking the real part, and integrating over T 2 gives 1 2 The Laplacian term on the RHS of (NLS) vanishes using integration by parts.By Lemma 2.7, the second term in (3.1) is bounded from below by the L p+2 x norm, so we get 1 2 Since we are in a domain of unit volume, Hölder's inequality leads to d dt It is now easy to conclude that the mass of superfluid (using the quantum mechanical interpretation of the wavefunction) decays algebraically in time.Namely, where x is the initial mass of the superfluid.□ 3.2.Energy estimate.In this subsection (Section 3.2), we derive the governing equations for the energy In Section 3.3, we work with a higher-order energy X(t), combining it with E(t) in Section 3.3.3.We begin by acting with the gradient operator on (NLS), multiplying by ∇ ψ, and taking the real part.This gives Integrating over T 2 , we observe that the first term on the RHS vanishes upon integration by parts due to the periodic boundary conditions.The second term on the RHS is similarly integrated by parts to yield 1 2 Now, we rewrite the first term on the RHS by expressing the Laplacian in terms of the operator B, giving us a dissipative contribution to the energy estimate.Namely, |ψ| p Re( ψBψ). (3.7) We also have to account for the potential (self-interaction) energy of the wavefunction.To obtain this, we multiply (NLS) by 2 ψ and take the real part to obtain Multiplying the above equation with µ|ψ| p and integrating over T 2 leads to 2µ p + 2 The terms on the RHS are canceled once we include the energy of the normal fluid.We first rewrite (NSE) in the non-conservative form, and apply the Leray projector (see Remark 2.2) to get Here, P is the Leray projector, which projects a Hilbert space into its divergence-free subspace, thus removing any purely gradient terms.We also apply the Leray projector to (NSE) to obtain Taking the inner product of both (NSE') and (NSE-L) with u, using incompressibility, and adding them, we arrive at the energy equation for the normal fluid, 1 2 Therefore, by adding (3.9) and (3.10), we obtain the energy equation Thus, the energy is bounded from above as with denoting the initial energy of the system.Next, we wish to show that the energy actually decays algebraically in time, under a certain smallness condition on the initial data.First, note that where we used an argument similar to the one from the proof of Lemma 2.7 to get the last inequality.We now use (2.9) to see that We then bound the first term on the RHS using Hölder inequality and Gagliardo-Nirenberg (GN) interpolation as For the second term in (3.15), we interpolate the L 3 x norm, while also applying the Hölder, Poincaré, and Young inequalities, as well as the GN interpolation inequality, to get (3.17) For sufficiently small values of κ and E 0 , the RHS of (3.17) can be absorbed into the LHS of (3.15).We also use the Poincaré inequality to convert the last term on the LHS of (3.15) into a coercive term for the internal energy term 2µ p+2 ∥ψ∥ p+2 L p+2 x in E(t).To this end, we observe that (3.18) In the last inequality, we interpolated between the L p+2 x and L 2 x norms, which may be done when p > 2. By choosing κ sufficiently small, we can absorb the second term on the RHS into the LHS.For p ≤ 2, we can simply replace ∥ψ∥ p+2 While we have the required coercive terms on the LHS, we cannot yet obtain a decay estimate for E(t), since the second term on the RHS is out of reach using E only.In order to control it, we set up an analogous inequality for a higher-order energy.
3.3.Higher-order energy estimate.In this subsection, we obtain further bounds for ψ and u, this time with one more derivative than the energy E.
Once again, the first term on the RHS of (NLS) vanishes due to the boundary conditions.We now estimate the terms on the RHS of (3.20).For the first term, x , which gives a dissipative term for ψ.For the term I 4 , we again integrate by parts, followed by Hölder's inequality to obtain x .
Thus, (3.20) becomes (3.21) The first of these terms is bounded as using the Poincaré and GN interpolation inequalities.We have also applied Young's inequality to extract out dissipative terms in the last step.We again use κ to denote a small number whose value shall be fixed later on, and C κ is a constant whose value depends on κ and the system parameters.
Similarly, for the second term on the RHS of (3.21), we have Finally, we apply the Sobolev embedding and Poincaré inequalities to bound I 7 .This leads to Combining all these inequalities into (3.21)results in where we have absorbed κ D 3 ψ 2 L 2 x into the LHS with a sufficiently small κ.

3.3.2.
The Navier-Stokes equations.We shall now derive a higher order estimate for the velocity field, which shall be combined with (3.25).Starting with (NSE'), we first multiply it by ∂ t u and integrate over the domain to obtain , we control the RHS.For the first term, x .In going to the last inequality, we used the GN interpolation and Poincaré inequalities.Finally, Young's inequality lets us extract the required dissipative term.For the second integral in (3.26), where the Bψ term is handled via the GN interpolation and Young's inequalities.In the third integral in (3.26), x ∥ψ∥ 2 L 6 x ∥Bψ∥ 2 x , where the term Bψ is handled just like in I 9 .Finally, for the last term in (3.26), Re(ψBψ)|u| 2 . (3.27) We estimate the second term on the RHS of (3.27) using the Hölder and GN interpolation inequalities.This gives Similarly, for the third term in (3.27), αλ Substituting the above estimates into (3.26),we arrive at where C κ depends on κ and the system parameters.
So far, we have obtained equations for ∥∇u∥ L 2 x and ∥∆ψ∥ L 2 x , while including the higher-order dissipation corresponding to the wavefunction, ∥∇(Bψ)∥ 2 x .What remains is to consider the higher-order velocity dissipation ∥∆u∥ 2 L 2 x .To this end, we multiply (NSE') by −θ∆u, with θ to be determined, and integrate over the domain.This gives x with a small coefficient, so it can be absorbed into the LHS.Thus, we have The second integral is manipulated just as I 8 and yields x .The bound for the integral I 14 follows from the GN interpolation, Poincaré, and Young inequalities, as x .In a similar manner, we have Finally, for the last integral in (3.29), Thus, (3.29) becomes x . (3.30) We now add (3.25), (3.28), and (3.30).We also observe that where the last three terms on the RHS are the same as I 5 , I 6 , and I 7 in (3.21).We bound them just as in (3.22)-(3.24).Choosing θ sufficiently small, and subsequently κ also small enough, we absorb x and ∥∆u∥ 2 L 2 x on the RHS into the corresponding terms on the LHS.Finally, what remains is where we absorbed x with an appropriate choice of κ.This is the higher-order energy estimate.

The Grönwall inequality step.
Having derived the equations for the higher-order norms of u and ψ, and while accounting for the relevant dissipative terms, the goal now is to use a Grönwalltype argument.Lemma 3.2 (Algebraic decay rate for energies).The sum of the energy E(t) and the higher-order energy X(t) := ∥∆ψ(t)∥ 2 x decays algebraically in time as (1 + t) Proof.We begin by denoting x , so we can rewrite (3.31), after updating θ, κ, E 0 , and S 0 to be sufficiently small, as where Q 1 (X + E) is a strictly super-linear polynomial, while Q 2 (X + E) contains both linear and super-linear terms.To arrive at (3.32), we have also expanded the Sobolev norms as for the velocity, and ∥ψ∥ 2 for the wavefunction.Next, we add (3.11) and (3.32) to end up with We use the Poincaré inequality to rewrite Y in order to get decaying norms.Indeed, x ≳ X.Additionally, we also use the analysis in (3.14 x , which in turn can be downgraded to ∥∇ψ∥ 2 L 2 x using the Poincaré inequality.One can also represent ∥Bψ∥ 2 L 2 x on the RHS of (3.35) by where we have used the estimates (3.16) and (3.17), and the GN inequality.After all of the above manipulations, (3.35) now reads where β depends on the system parameters, and the polynomials Q 1 and Q2 are strictly superlinear.The first term on the RHS results from the estimates in (3.18).As for the second term on the RHS, we note that this can be absorbed into the LHS by tweaking S 0 .
For notational convenience, we write Z := X + E and use Q := Q 1 + Q2 to denote the strictly super-linear polynomial in the RHS of (3.36), leaving us with The Duhamel solution for Z(t) obeys We set the size of the initial data as Z(0) =: Z 0 ≤ ε.We need to use a bootstrap argument to show that Z(t) ≤ 3ε for t ∈ [0, T ].Specifically, we prove that the hypothesis Z(t) ≤ 4ε for t ≤ t 1 leads to the stronger conclusion Z(t) ≤ 3ε for t ≤ t 1 , where t 1 ∈ [0, T ].To this end, we estimate each integral on the RHS of (3.38).The first integral is split into two parts to take advantage of the exponential decay factor.Thus, The last inequality is a result of the exponential decay of the first term, compared to the algebraic decay of the second.The second integral in (3.38) is more straightforward, with Now we choose ε small enough (call it ε 0 ) so that the RHS ≤ ε.Similarly, the contribution from the first integral term in (3.38) is made less than ε for all S 0 ≤ ε 0 < 1 small enough.This completes the bootstrap argument, and we can see that indeed Z(t) ≤ 3ε.For ε 0 sufficiently small, the linear dissipation in (3.37) dominates the nonlinearities, and we may write the equation as whose solution, following (3.37)-(3.39),obeys Returning to (3.35), we now absorb the last term on the RHS into the LHS, which is possible for small enough data since sup 0≤t≤T Z(t) ≤ 3ε 0 .Furthermore, in the regime of small data, the super-linear polynomial Q 1 can be dominated by the linear term on the RHS, which leaves us with Employing the bound for Z from (3.41) in the RHS of (3.42) and integrating over [0, T ], we estimate the dissipation as (3.43)The last inequality holds because S 0 < 1.Thus, we can achieve small values for the RHS of (3.43) by selecting appropriate Z 0 and S 0 .
Another useful estimate for the dissipative terms results from integrating (3.42) over the time interval [t, 2t], where t ≥ 1.This gives (3.44)This time-decaying bound on the dissipation is necessary to obtain a sharp control of the dynamics at large times.□

3.4.
The highest-order a priori estimate for ψ.From the previous analysis, we have obtained Bψ ∈ L 2 [0,T ] H 1 x .However, as pointed out in the discussion following Definition 2.9, we seek To this end, we would like to obtain an even higher order a priori estimate, only for ψ.Lemma 3.3 (Algebraic decay rate for highest-order norm of ψ).For S 0 , E 0 , and Z 0 small enough, and with s = 5 4 , the homogeneous Sobolev norm ∥ψ(t)∥ Ḣ2s x decays algebraically in time as (1+t) x is sufficiently small, the higher-order dissipation ∥ψ∥ L 2 [0,T ]

Ḣ2s+1
x may be made as small as required, independent of the time T .
x .This inclusion is precisely what we need to control the term |u| 2 ψ in the coupling using the a priori estimates up to this point.
Proof.With s = 5 4 , we apply (−∆) s to (NLS), to get Just as in Sections 3.2 and 3.3.1,we multiply by (−∆) s ψ, and integrate the real part over T 2 .As a result, the second term on the LHS yields Re Using the self-adjoint property of the Laplacian allows us to conclude that the first term on the RHS of (3.45) vanishes, since For the second term on the RHS of (3.45), we have In the last expression of (3.46), we expand the operator B and use Hölder's inequality to arrive at (3.47) Rewriting the LHS in terms of the homogeneous Sobolev norms and the RHS in terms of the usual Sobolev norms, we have Since 2s − 1 = 3 2 , the algebra property of Sobolev norms is applicable.Using this, (3.4) and (3.41), we estimate the RHS of (3.48).The first term requires interpolation4 and yields where we have retained only the terms that decay the slowest.In arriving at the last inequality in (3.49), we use the fact that Z 0 < 1.For the second term in the RHS of (3.48), we similarly obtain While the H 3 2 x norm could have been interpolated between H 1 x and H 2 x , it does not provide an improved estimate since ∥u∥ x are both bounded by (3.43).In the last term of (3.48), in view of Remark 2.5, we have where the penultimate inequality is obtained using (3.4) and (3.41).Therefore, (3.48) becomes With the Poincaré inequality, we replace the dissipative term on the LHS with W (t) := ∥ψ(t)∥ 2

Ḣ2s
x .We employ calculations similar to (3.39) to estimate the integrals, i.e., splitting them over 0, t 2 and t 2 , t .We also use (3.43), in particular ∥u∥ 2 ≲ 1, to simplify the exponential factors outside the integrals.In all, we end up with for all t ∈ [0, T ].We simplify further by making use of (3.43) and (3.44) for ∥u∥ 2 x , respectively.This leads to We use (3.55) in (3.52) and integrate over [0, T ] to obtain the final dissipative estimate This shows that with small enough data one can achieve an arbitrarily small value (independent of T ) for this highest-order dissipation.
Similarly to (3.44), it is possible to also get a time-decaying estimate by integrating (3.52) over the time interval [t, 2t] for t ≥ 1.This leads to where we have used (3.44) and (3.55), and retained only the slowest decaying terms.□ The high-norm control in (3.56) and (3.57) is important because the inequalities can be translated into the desired bounds (on two fewer derivatives) for Bψ.Indeed, where for the last three terms, we replaced the homogeneous Sobolev norms by the larger inhomogeneous norms.Combining the analysis in (3.49)-(3.51)with (3.4), (3.43), (3.56), and (3.58), we get the sought-after dissipation bound The estimates in (3.59) and (3.60) are used to ensure that the density remains bounded from below.
3.5.Ensuring positive density.We now have all the a priori estimates to return to (2.15).For it to hold true, a sufficient condition is Depending on the value of p, we now divide the analysis into several cases: 1 ≤ p < 2, p = 2, 2 < p < 4, p = 4, and p > 4.
3.5.1.The case 1 ≤ p < 2. Owing to the Poincaré inequality and (3.43), we have and this bound holds for all p ≥ 1.For the first term of (3.61), we integrate (3.4), yielding since 2 p > 1.From (3.59), (3.62), and (3.63), we conclude that the condition in (3.61) can be achieved if is sufficiently small.Thus, the density satisfies m f ≤ ρ ≤ M i + m i − m f for all T > 0, as long as the initial data are small enough.
For p ≥ 2, the integral of the superfluid mass, i.e., ∥ψ(t)∥ 2 L 2 x , cannot be bounded uniformly in [0, T ].This is where the decaying estimates in (3.44) and (3.60) prove to be useful.
3.5.2.The case p = 2.We split the time integral in (3.61) over the ranges 0 ≤ t ≤ 1 (short-time) and t ≥ 1 (long-time).We start with the long-time estimate the LHS of (3.61) with p = 2.For the first term, we have Using the Poincaré inequality and (3.44) gives . This leads us to which is the long-time contribution (independent of N ) of the constraint in (3.61).It can be made as small as required with an appropriate choice of W 0 + Z 0 + S 0 .Finally, we verify the short-time control as well.The superfluid mass bound in (3.4) means that Similarly, using (3.43), we get From (3.59), (3.67), and (3.68), we have which can be made small enough to satisfy (3.61).This lets us conclude that the density is bounded from below uniformly in time, for the case p = 2. Thus, we have the necessary global bound.
3.5.3.The case 2 < p < 4. We begin, once again, with the long-time analysis, i.e, for t ≥ 1.From (3.4), we have Using the Poincaré inequality and (3.44), .71), we have Once again, the slowest decaying term is the dominant one.Therefore, we have (3.73)The sum converges (uniformly in N ) because p < 4. Hence, we obtain good long-time control of the LHS of (3.61) for 2 < p < 4.
What remains is to check that we also maintain short-time control.To this end, we have from (3.4), and from (3.43), which is the short-time control we are seeking.This implies global solutions, since the density is bounded from below uniformly in time.
3.5.4.The case p ≥ 4. The arguments for short-time control in Section 3.5.3remain valid even for p ≥ 4.However, the long-time estimates breaks down.Specifically, the geometric series in (3.73) diverges.We see that for , p > 4.
(3.76) Therefore, in this scenario, global-in-time estimates elude us due to the logarithmic/polynomial dependence on T .We can, however, guarantee almost global existence of solutions.Given a set of system parameters, we can ensure that ρ ≥ m f for any finite time T > 0 as long as we start from small enough initial data (depending on T ).In other words, if the size of the data is ε, then we have T ∼ e Having derived the required a priori estimates, we now establish the existence of a weak solution for a truncated form of the governing equations, and then pass to the limit.4.1.Constructing the semi-Galerkin scheme.The finite-dimensional wavefunction and velocity are constructed using eigenfunctions of the Laplacian and the Leray-projected Laplacian, respectively.4.1.1.The approximate wavefunction.Consider the negative Laplacian −∆ on the torus T 2 , with the domain D(−∆) = H 2 .It has a discrete set of non-negative and non-decreasing eigenvalues {β j }, and the corresponding eigenfunctions {b j } ∈ C ∞ (T 2 ) can be chosen to be orthonormal in L 2 x and orthogonal in H 1 x .We define the approximate wavefunction as for N ∈ N ∪ {0} and d N k (t) ∈ C.
4.1.2.The approximate velocity.We consider the Leray-projected Laplacian (or Stokes operator) A = −P∆ with the domain D(A) = L 2 d ∩ H 2 (see [RRS16, Chapter 2], for instance).The Stokes operator (like the Laplacian) has a discrete set of non-negative and non-decreasing eigenvalues {α j }, and the corresponding divergence-free, vector-valued eigenfunctions {a j } ∈ C ∞ (T 2 ) can be chosen to be orthonormal in L 2 d,x and orthogonal in H 1 x .We define the approximate velocity as for N ∈ N ∪ {0} and c N k (t) ∈ R. 4.2.The initial conditions.4.2.1.The initial wavefunction and initial velocity.We begin by defining P N (respectively, Q N ) to be the projections onto the space spanned by the first N + 1 eigenfunctions of A (respectively, −∆).Then, we truncate the initial conditions for the velocity and wavefunction accordingly: ), it is necessary to establish that the truncated initial conditions converge to the actual ones in the relevant norms.
ψ, and (2) The proof utilizes the equivalence of norms between Sobolev spaces and fractional powers of the negative Laplacian/Stokes operator (see Theorem 2.27 in [RRS16]).Given the regularity of ψ 0 and u 0 , we deduce the convergence of the approximate initial conditions by applying Lemma 4.1.

4.
3.1.The continuity equation.Having described the (approximate) initial conditions and the semi-Galerkin scheme, we now establish the existence of solutions to the "approximate" equations, starting with the continuity equation.It is given by where Just as in (2.15), we see that the constraint that fixes the local existence time Since the norms in (4.5) are bounded by the size of the initial data, the time T N is independent of N .Hence, we use T to denote the time of existence, with T arbitrarily large for 1 ≤ p < 4 and T bounded for p ≥ 4 (as specified in Theorem 2.4).We now establish the analogs of Lemmas 2.2 and 2.3 from [Kim87].These constitute the existence of a unique solution to (4.4) and a Picard iteration scheme for the same, respectively.
Proof.Consider the evolution equation for the characteristics of the flow, x N (0) = y N ∈ T 2 .
x N is well-defined.We now write the solution to (4.4) along characteristics as That (4.7) uniquely solves (4.4) can be verified using the property of the "inverse-characteristics" y(t, x).For any τ ∈ R, where the last equality is due to Euler's chain rule.□ Now, we consider a convergent sequence of velocities and wavefunctions that belong to the finitedimensional subspaces spanned by the truncated Galerkin scheme.Given such a convergent sequence, we show that the sequence of density fields satisfying (4.4) is also convergent, and this shall be used to complete a contraction mapping argument below.
x the unique solution to the system , where ρ N solves (4.4).
Proof.We begin by defining Ψ Given the convergence of y N n derived above, and because ρ N 0 ∈ C 1 x , the first term on the RHS vanishes.The second and third terms vanish on account of the following argument.Note that Ψ N n has its highest order term of the form ψ N n ∆ψ N n (second derivative), and so the assumed convergence of ψ N n in the C 0 x is finite, uniformly in n. □ 4.3.2.The Navier-Stokes equation.We now consider an "approximate momentum equation", composed of the approximate wavefunction and velocity fields defined by (4.1) and (4.2), respectively.Namely, .11) Recall that the incompressiblity condition is built-in, because the eigenfunction basis used to construct the velocity fields are divergence-free.Now, taking the L 2 inner product of (4.11) with a j (x) for 0 ≤ j ≤ N , we arrive at a system of equations for the coefficients describing the time-dependence of the approximate velocity fields, as where Since we have both lower and upper bounds on the density in the chosen interval of time, we can use Lemma 2.5 in [Kim87] to show that the matrix R N (t) is invertible.Therefore, we arrive at which is the desired evolution equation (written vectorially) for the coefficients c N j (t).

4.3.3.
The nonlinear Schrödinger equation.As in the previous section, we derive an evolution equation for the coefficients of the approximate wavefunction, by considering an "approximate NLS".Namely, Recall that B L = B − µ|ψ| p , i.e., the linear (in ψ) part of the coupling operator.Performing an L 2 inner product with b j (x), we get where Written vectorially, the evolution equation for the coefficients d N j (t) becomes where 4.4.Compactness arguments.We now extract convergent subsequences from the a priori estimates in Section 3. Beginning with the density, we know that ρ Moreover, from (4.4),

.18)
The second inequality is due to the (compact) embedding L 2 x ⊂ H −1 x for T 2 .All the terms in the last line are finite (uniformly in N ) by virtue of the a priori estimates.Therefore, using the Aubin-Lions-Simon lemma, we conclude the strong convergence of a subsequence of the density as Consider a relabeled subsequence ρ N that strongly converges to ρ in C([0, T ]; H −1 x ), so that (4.1) and (4.2) are also appropriately relabeled.For a.e.s, t ∈ [0, T ] and any ω ∈ H showing that ⟨ρ N (t), ω⟩ H −1 ×H 1 is uniformly continuous in [0, T ], uniformly in N due to (4.18).Due to the embedding H 1 x ⊂ L r x for all 1 ≤ r < ∞, we conclude, using the Arzela-Ascoli theorem, that ρ N is relatively compact in C w ([0, T ]; L r x ).We move on to the velocity.Based on the a priori estimates, we extract a subsequence of u N that weakly converges to u Applying the Lions-Magenes lemma (see [Tem77, Chapter 3]), we deduce that u ∈ C([0, T ]; H 1 d,x ).Based on the L ∞ t L ∞ x bound on the density, and the above strong convergences, it is easy to see that ρ N u N and ρ N u N ⊗ u N converge in C([0, T ]; L 2 x ) to ρu and ρu ⊗ u, respectively.Next, we consider the wavefunction.Again, we extract a subsequence that converges weakly to x ∩ L 2 [0,T ] H 7 2 x .From this and (NLS), we have ∂ t ψ ∈ L 2 [0,T ] H 3 2 x .Thus, the Lions-Magenes lemma yields ψ ∈ C([0, T ]; H x and H 1 d,x , respectively.For the momentum, we have where 1 r + 1 r ′ = 1 2 .Using the embedding H 1 x ⊂ L r ′ x to handle the velocity in the first term of the RHS, it is easy to see that the initial momentum converges in the L 2 x norm.The approximate solutions (ψ N , u N , ρ N ) are smooth enough to satisfy (2.1)-(2.3).The aforementioned compactness results allow us to pass to the limit of N → ∞ and arrive at the weak solutions (ψ, u, ρ).4.5.Renormalizing the density.At this point, we know that ρ N * − ⇀ ρ in L ∞ t L ∞ x .We wish to use the technique of renormalization to extend this to ρ N → ρ in C 0 t L r x , for 1 ≤ r < ∞.To achieve this, we will adapt a classical argument (see, for instance, Theorem 2.4 in [Lio96b]).We begin by defining a sequence of unit-mass mollifiers ζ h (x) = 1 h 2 ζ x h , where h will eventually be taken to 0. Next, for a given weak solution ρ ∈ L ∞ t L ∞ x , we mollify (CON) to obtain where g h := g * ζ h , Ψ := 2λ Re(ψBψ), and R h := u • ∇ρ h − (u • ∇ρ) h is a commutator.We multiply this by η ′ (ρ h ), for a C 1 function η : R → R.This yields The Sobolev embedding H 2 x ⊂ W 1,r 1 x for any r 1 ∈ [1, ∞) implies that u ∈ L 2 t W 1,r 1 x .From Lemma 2.3 in [Lio96b], we note that R h vanishes in L 2 t L r 1 x (and also in L ∞ t L 2 x ) as h → 0, by choosing r 1 > 2. Similarly, Ψ h converges to Ψ in C 0 t L 2 x .Finally, note that η ′ (ρ h ) is uniformly continuous since ρ (and ρ h ) take values in a compact subset of R. Therefore, using a test function σ, we may pass to the limit h → 0 in (4.22).In other words, if ρ is a weak solution of the continuity equation, then η(ρ) solves (in a weak sense) ∂ t η(ρ) + u • ∇η(ρ) = η ′ (ρ)Ψ.(4.23)This is the renormalized continuity equation.Taking the difference of (4.21) for h 1 , h 2 > 0, we write the analog of (4.22) for η(ρ h 1 − ρ h 2 ), with η(x) = x 2n where n ∈ N. Integrating over T 2 leads to x , it follows from the Sobolev embedding and Hölder's inequalities that Ψ = ψBψ ∈ L 1 t L r 1 x for any r 1 ∈ [1, ∞).Between this, the commutator estimate in Lemma 2.3 of [Lio96b], and the boundedness of ρ 0 , we find that all of the terms on the RHS of (4.24) vanish as h 1 , h 2 → 0, giving us a Cauchy sequence in C([0, T ]; L 2n x ).Hence, ρ h converges to ρ in C([0, T ]; L 2n x ).We have, so far, proved that our "original approximations" of the continuity equation ρ N converge in C w ([0, T ]; L r x ) to ρ, and that ρ also belongs to C([0, T ]; L 2n x ).To achieve what we set out to prove, i.e., that ρ N converges strongly in C([0, T ]; L r x ) to ρ, it remains to show that the L r x norms are continuous in time.It is sufficient to illustrate this for r = 2 (or n = 1), in order to deduce it for the other values of r.Explicitly, if there is a sequence of times t N → t, then we need ρ N (t N ) to converge in L 2 x to ρ(t).Returning to (4.4), we look at its renormalized version with η(x) = x 2 , and integrate over T 2 (and then from 0 to t N ) to get Since we know that ρ ∈ C([0, T ]; L 2 x ), we can do the same calculation with (CON), except over the time interval 0 to t.This yields x , we can use the strong convergence in (4.19) to handle the first term on the RHS.The second and third terms follow from simple Hölder's inequalities, and the strong convergence of ψ N of B N ψ N .Finally, the last term is integrable on [0, T ], so as t N → t, it vanishes.In summary, which, along with the weak-in-time continuity deduced earlier, implies strong convergence of ρ N to ρ in C 0 t L 2n x for all n ∈ N. Interpolating between Lebesgue norms extends this result to C 0 t L r x for all r ∈ [1, ∞).
4.6.The energy equality.The smooth approximations to the weak solutions satisfy an energy equation, given by (2.8), i.e., 1 2 ρ N (t)u N (t) for a.e.t ∈ [0, T ].From our choice of the initial conditions and their approximations (see Section 4.2), we can ensure that as N → ∞, the RHS converges to the initial energy E 0 defined in (3.13).Indeed, for the first term, Moreover, based on the results of Section 4.4, we can conclude that all the terms on the LHS of (4.26) converge strongly to the corresponding terms with the approximate solutions replaced by the weak solution.The first term on the LHS can be dealt with the same way as the first term on the RHS in (4.27), by simply including a sup t outside the absolute values.□ This completes the construction of the solutions.Together with the global/almost global estimates from Section 3, we can conclude the results of Theorems 2.3 and 2.4.
L p 2 +1 by ∥ψ∥ p+2 L 2 x since we are on a finite-size domain.Thus, irrespective of the value of p, (3.15) becomes

ε 4 .
for p = 4 and T ∼ ε − p p−4 for p > 4.This is the scaling expressed in Theorem 2.4.Existence of weak solutions (Proof of Theorems 2.3 and 2.4)

Lemma 4 .
1 (The projections Q N and P N are convergent).If ψ ∈ H r x and u ∈ H s d,x for any 0

(4. 6 )
Since u N ∈ C 0 [0,T ] C 1 x , there exists a unique solution x N (t, y N ) ∈ C 1 [0,T ] C 1 x .Owing to the incompressibility of the flow u N , it follows that det allowing us to conclude that the characteristics are C 1 diffeomorphisms and therefore, invertible.This means

4. 3 . 4 .
Fixed point argument for the coefficients.For a fixed N , a standard contraction mapping argument shows that (4.13) and (4.16) have unique solutions that are continuous in [0, T ].For a pair (u N n , ψ N n ), equivalently (c N n , d N n ), using Lemma 4.2, we can find a solution ρ N n .Owing to the smoothness (in space) of the eigenfunctions used in the approximate velocity and wavefunction, we conclude that u Therefore, performing an iteration on the triplet (c N n , d N n , ρ N n ) and using Lemma 4.3, we conclude that the sequence ρ N n converges to ρ N ∈ C 0 [0,T ] C 0 x .

x
).Additionally, we also haveB N ψ N C 0 t L 2 x − −−→ Bψ,due to the regularity of u and ψ.As for the initial conditions, by construction itself (Section 4.2.2),we have ρ N 0 L r x −→ ρ 0 for 1 ≤ r < ∞.Also, Lemma 4.1 states that ψ N 0 and u N 0 converge to ψ 0 and u 0 in H 5 2

ρ
Re ψBψ .Subtracting the last two equations, and taking the limit N → ∞, we observe that the first terms on the RHS cancel (recall from Section 4.2.2 that ρ N Thanks to the uniform boundedness of ψN B N ψ N in L 1 [0,T ] H 3 2 this is just the inverse of the characteristic, i.e., if the flow were reversed.Due to the flow being incompressible, we know that the matrix ∂y N n ∂x is invertible.Also, as shown in the proof of the previous lemma, ∂ ∂t y N n = −u N n • ∇ x y N n .This implies that the derivatives of y N n with respect to both space and time are bounded uniformly in n, t and x.Thus, by the Arzela-Ascoli theorem, we can extract a subsequence that converges N .Just as in (4.8), we can show that the solution to (4.9) is n→∞x N .Consider the map y → x N n (t, y) and define its inverse y N n (t, x); n→∞ y