On the mass transfer in the 3D Pitaevskii model

We examine a micro-scale model of superfluidity derived by Pitaevskii in 1959 which describes the interacting dynamics between superfluid He-4 and its normal fluid phase. This system consists of the nonlinear Schr\"odinger equation and the incompressible, inhomogeneous Navier-Stokes equations, coupled to each other via a bidirectional nonlinear relaxation mechanism. The coupling permits mass/momentum/energy transfer between the phases, and accounts for the conversion of superfluid into normal fluid. We prove the existence of weak solutions in $\mathbb{T}^3$ for a power-type nonlinearity, beginning from small initial data. The main challenge is to control the inter-phase mass transfer in order to ensure the strict positivity of the normal fluid density, while obtaining time-independent a priori estimates.


Introduction and mathematical model
In this article, we present a rigorous analysis of the Pitaevskii model (a micro-scale description) of superfluidity [Pit59] in three dimensions.The system consists of a superfluid phase and a normal fluid phase, described by modified versions of the nonlinear Schrödinger equation (NLS) and the Navier-Stokes equations (NSE), respectively.This is one of many different theories proposed to explain and quantify the underlying mechanisms of superfluidity.For more details, see [BDV01, Vin06, PL11, BBP14, JT21, Jay22] and references therein.The Pitaevskii model, specifically, works well in the context of small length scales (≪ inter-vortex spacing), and has previously been explored in [JT22a,JT22b,JJK23].A similar model has also been numerically simulated in [BSZ + 23].The superfluid phase is represented by a complex-valued wavefunction ψ, while the normal fluid is characterized by its density ρ, velocity u, and pressure q.The form of the Pitaevskii model used here is as follows (with the prefix "c" in the equation labels signifying that the equations are coupled): The strength of the superfluid's scattering self-interactions is measured by µ > 0, while ν > 0 and α > 0 denote the viscosity and drag coefficient of the normal fluid.The quantum scattering is a power-type nonlinearity, with exponent p ∈ [1, ∞).The interactions between the two phases are mediated by the nonlinear coupling operator B, and the strength of the coupling is quantified by the positive constant λ.The structure of this coupling permits bidirectional mass/momentum/energy transfer between the two phases, which results in a relaxation mechanism for the otherwise nondissipative NLS.Indeed, the coupling gives (c-NLS) a parabolic flavor.These equations are supplemented with the initial conditions ψ(0, x) = ψ 0 (x), u(0, x) = u 0 (x), ρ(0, x) = ρ 0 (x) a.e.x ∈ T 3 .

(INI)
We use periodic boundary conditions, i.e., we are working in a 3-dimensional torus [0, 1] 3 .The above equations are a slight modification of Pitaevskii's original work [Pit59], which were valid for any type of scattering interactions and for a compressible, heat-conducting normal fluid.The simplifying assumptions used here are detailed in [JJK23].While (DIV) implies that the velocity is divergence-free, this does not mean that the density is constant along particle trajectories.In fact, the density is governed by (c-CON), a transport equation with a complicated source term.This inhomogeneity on the RHS is the principal limiting factor of the analysis, since it forces us to ensure that the normal fluid density does not become zero (vacuum) or negative (unphysical).By integrating (c-CON) over T 3 , the advective term vanishes and using the positivity of the operator B [JJK23, Lemma 2.7], we have This implies that the net mass of the normal fluid does not decrease with time.Due to conservation of total mass (of both fluids), we have In other words, the coupling results in a conversion of superfluid into normal fluid, on average.However, the RHS of (c-CON) is not necessarily non-negative for all x ∈ T 3 , and we have to control the local mass transfer between the two phases.Indeed, this means bounding the L ∞ norm (in space) of the source term, Re(ψBψ), as described in Section 2.2.Since the coupling is a secondorder differential operator, this essentially translates into finding an estimate of ∥ψ∥ L ∞ x ∥∆ψ∥ L ∞ x , which dictates the required regularity of the wavefunction.In the case of no source term for the continuity equation, there are several results pertaining to the standard version of the incompressible, inhomogeneous NSE -see [Kaz74, LS78, Kim87, Sim90, Lio96, Dan03, CK03, Dan06] among many other references.
The NLS, an archetypal dispersive PDE, is often used to study quantum systems with low-energy wave interactions, like dipolar gases [CMS08,Soh11] and quantum ferrofluids [BSG + 18].There is a rich literature on the mathematical analysis of the NLS (see [Tao06] for a collection of results).One particular aspect of NLS that has attracted much attention is its quantum hydrodynamic (QHD) reformulation [JMR02, Jün10, AS17, AMZ21, WG21, AMZ23].Another connection with the compressible fluid dynamics community has been the study of Korteweg models [HL94, HL96, BD04, CDS12], of which QHD and even capillary flows [AS22] are special cases.The linear drag included in (c-NSE) is not unheard of: similar terms have been used previously to either prove the global existence of solutions [Cha22], or to show relaxation to a steady state [BGLVV22,SYZ22].For QHD with a combination of linear drag and electrostatic forces, see [JL04,AM09,AM12].
Given the vast mathematical literature that exists on the NSE and NLS independently, it is surprising that there are hardly any on a combination of the two, like the Pitaevskii model or others [Kha69,Car96].The first attempt to study a combined model was by Antonelli and Marcati [AM15], in which a fractional time-step method pioneered by the same authors [AM09, AMS12] was utilized.In this approach, the standard, uncoupled NLS is solved over a small time interval at the end of which the wavefunction is "updated" to account for the interactions with the normal fluid.However, this was still a uni-directionally coupled system, i.e., the wavefunction ψ was dependent on the density ρ and velocity u, but not the other way around.This yields the standard form of the continuity equation: without a source term.
The small-data global existence of solutions for the Pitaevskii model on T 2 was established in [JJK23] for the nonlinearity exponent p ∈ [1, 4).Moreover, for p ≥ 4, an almost-global existence was shown, wherein the existence time grows with decreasing data size, exponentially for p = 4 and polynomially for p > 4. The path to the required a priori bounds of ψ is by energy-type estimates, up to L ∞ t H 5 2 x .In the current work, we obtain the global existence of weak solutions in T 3 , for the entire range of power nonlinearities (characterized by 1 ≤ p < ∞), thus proving that the inter-phase mass transfer may be controlled for small initial data.We emphasize that we deal with density that is simply in L ∞ x , and thus not necessarily differentiable.This limits the highest regularity that one can subject the momentum equation (and the velocity) to.In turn, through the coupling operator B, this also restricts the regularity we can achieve for the wavefunction.The primary objective is to derive time-independent a priori bounds on ψ that ensure that the RHS of (c-CON) will not lead to non-positive densities.
The approach used for 2D is not applicable in 3D, owing to the insufficiency of the x bound for the velocity.This regularity of u is not enough to derive an energy estimate of the required order for ψ in 3D.Any higher regularity of u would mean taking the derivatives of ρ, which is off limits.Another roadblock is the extremely slow mass decay for high values of p, i.e., when the superfluid interacts much more strongly with itself rather than the normal fluid.In such a scenario, a purely energy-based method [JJK23] led to a bifurcation in the existence time of solutions (global or almost-global) for different ranges of the nonlinear index p.The mass decay rate (independent of dimension) pervades all levels of energy estimates, dominating the decay of higher norms as well.
We overcome these challenges by a hybrid approach: combining the decay of superfluid mass with a maximal regularity estimate for parabolic equations (Section 2.3).The time-control of the mass conversion and of the higher order energy norm are presented in Lemmas 3.1 and 3.3, and state that the L 2 x and Ḣ2 x norms of the wavefunction decrease as (1 + t) p , respectively.These allow us to evaluate the integral in (2.13), independent of the final time T , leading to global control of the solution.In order to apply maximal regularity, the initial wavefunction ψ 0 must belong to an interpolation (Besov) space, which is carefully chosen to be marginally larger than H 2 , so that the assumption ψ 0 ∈ H 2 is sufficient for our purposes.The resulting solution ψ is shown to belong to C([0, ∞); H 2 ), among other spaces.To summarize, we demonstrate that wielding the power of parabolic regularity allows us to guarantee global solutions, even when the mass decay is exceedingly small.We now briefly outline the notation used in this article.Following this, we state and discuss the main result in Section 2. Several a priori estimates are derived in Section 3, which ends with an argument on ensuring a positive lower bound for the density.In this paper, we only present the required a priori estimates, as the general construction of solutions (and the density renormalization) follows as in [JJK23, Section 4].
We use the subscript x on a Banach space to denote that the Banach space is defined over T 3 .For instance, L r x := L r (T 3 ) and H s d,x := H s d (T 3 ).For spaces/norms over time, the subscript t is used, such as L r t .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 κ is used to represent a (small) positive number.The values of C and κ can vary across the different steps of calculations.

Main result and discussion
2.1.Weak solutions and the existence theorem.First, we define the notion of a weak solution used here.
Definition 2.1 (Weak solutions 1 ).For a given time ), and (ii) ψ, u, and ρ satisfy the governing equations in the sense of distributions for all test functions, i.e., and and Remark 2.2.The last two terms in (c-NSE) are pure gradients, and thus we can absorb them into the pressure, relabeling the latter as q.Due to the use of the divergence-free test functions, all gradient terms in the definition of the weak solution disappear.Now, we state the main result.
e. in T 3 .Then, there exists a global 1 See Remark 2.5.
weak solution (ψ, u, ρ) to the Pitaevskii model such that the density is always bounded between m f ∈ (0, m i ) and M f := M i + m i − m f , provided the initial data satisfy the smallness condition (2.4) The solution has the regularity for a sufficiently small δ 1 > 0, and 1 ≤ s ≤ 6.Additionally, the solution satisfies the energy equality (2.8) The proof of Theorem 2.3 is based on a priori estimates, a semi-Galerkin scheme to construct solutions, and an adaptation of the classical renormalization procedure for the density [Lio96, Theorem 2.4].Since the coupling operator B contains the velocity u by itself (and not in combination with ρ), we limit the calculations to when the density has a positive lower bound.This gives us an indication to the level of regularity expected of the RHS of (c-CON), which in turn defines the spaces in which ψ and u belong to.In order to achieve this, we derive the required a priori control for the wavefunction and velocity, while ensuring that the density is neither differentiated nor does it become zero anywhere in the domain.
Before a more specific discussion on the method of proof, a few remarks about the result are warranted.
Remark 2.4.Once we have ρ ∈ L ∞ t,x , the renormalization procedure from [DL89] can be used to show that ρ indeed belongs to C t L s x .In a finite 3D domain, the Sobolev embedding H 1 ⊂ L s for 1 ≤ s ≤ 6 accounts for the integrability in Theorem 2.3.It is worth mentioning that in the analogous result in 2D, we get 1 ≤ s < ∞ due to the more favorable embedding.
Remark 2.5.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 follow, and in fact also prevent us from concluding uniqueness of the weak solutions.See [JT22b] for results akin to weak-strong uniqueness for the Pitaevskii model.
Remark 2.6.In Lemma 3.1, we establish that the mass of superfluid decreases with time (algebraically) and goes to 0 as t → ∞.Due to overall mass conservation, this means an increase in normal fluid mass, and an eventual conversion of all the superfluid into normal fluid.This inter-phase mass transfer is one of the underlying physical phenomena that the Pitaevskii model was designed to explain.In the macro-scale models of superfluidity (like the HVBK equations [Hol01,JT21], the coupling between the two fluids is suggestively called mutual friction, as it dissipates the overall energy of the system.In the micro-scale model that we are concerned with, such an energy sink can be interpreted as "heating up" the Bose-Einstein condensate (superfluid particles) into excited states (normal fluid particles).
Remark 2.7.We point out that Theorem 2.3 is also valid for the 2D Pitaevskii model (with the same regularity, except for the density renormalization argument holding for all 1 ≤ s < ∞).
Remark 2.8.It would be interesting to consider the Pitaevskii model in R 3 .It may help to localize the dynamics by using an external confining potential (in the NLS) that rises in strength with increasing distance from the origin.This would be akin to trapped-ion quantum systems in condensed matter physics.In such a scenario, it is plausible to expect that the current results from T 3 would continue to be valid in R 3 as well.We thank one of the anonymous referees for posing this question.
2.2.The strategy.As indicated above, the main difficulty is to guarantee that if we begin from ρ 0 that is bounded below, then the time evolution does not result in a degeneration where ρ vanishes at certain points in the domain.Hence, we define our existence time so that ρ does not go below a fixed lower bound until t = T * .So we aim to show that the chosen lower bound can be maintained for an arbitrarily long time.
Definition 2.9 (Local existence time).Start with an initial density field 0 < m i ≤ ρ 0 (x) ≤ M i < ∞.Given 0 < m f < m i , we define the existence time for the solution as (2.9) Consider the Lagrangian path of a particle starting at y ∈ T 3 , as it is advected by the local velocity.These characteristic curves are denoted by X y (t) and solve the ODE given by where u is the velocity of the normal fluid.Traveling along such a curve, we observe that is a (formal) solution to the continuity equation.From (2.9) and (2.11), it is clear that a sufficient condition to ensure the density is bounded from below by (2.12)This is, in turn, guaranteed by the sufficiency By selecting small enough data so that all the arguments may be bootstrapped, it is possible to achieve (2.13) independently of T > 0. Since Bψ involves a second-order derivative, its L ∞ x norm translates into high-regularity Sobolev spaces.For u, this means showing that it belongs to x , which also proves useful in establishing strong continuity in time for the solution.As for the wavefunction, we make use of the parabolic nature of (c-NLS) to derive the necessary regularity (see Lemma 2.11 below).It is important to remember that throughout these calculations, we handle the density only in L ∞ x , and not in any derivative spaces.2.3.Elliptic operators and maximal parabolic regularity.We now define uniform ellipticity in the context of complex-valued Banach spaces, before stating the maximal parabolic regularity result that will be utilized in Lemma 3.4.Definition 2.10 (Uniform (K, ζ)-ellipticity [PS01]).For a complex-valued Banach space X, consider the differential operator with domain D(A(t, x)) ⊂ X, where α is a multi-index and ∂ α denotes spatial derivatives.The coefficients a α are bounded and uniformly continuous functions from [0, T ] × R n to C N ×N for some n, N ∈ N. The principal symbol associated with this operator is (2.15) The operator A(t, x) is said to be uniformly It is possible to show that a uniformly (K, ζ)-elliptic operator generates an analytic semigroup of negative type, leading to the maximal regularity below.

A priori estimates
We now 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.The derivations of some estimates are identical to the 2D case [JJK23], and are not repeated here.
3.1.Superfluid mass estimate.Lemma 3.1 (Algebraic decay rate of superfluid mass).The mass S(t) of the superfluid decays algebraically in time as (1 + t) where x is the initial mass of the superfluid.Proof.See Lemma 3.1 of [JJK23].□ 3.2.Energy estimate.The energy of the system is defined as Lemma 3.2.The energy balance of the system is given by Proof.See Section 3.2 of [JJK23].□ Integrating (3.3) in time, we observe that the energy is bounded from above as where denotes the initial energy of the system.Next, we show that the energy is not just bounded, but also decays (algebraically) with time.In order to achieve this, we observe that ∥Bψ∥ L 2 x can be rewritten as D 2 ψ L 2 x , at the expense of some nonlinear terms on the RHS.More precisely (see equation (3.14) in [JJK23] for the exact derivation), Thus, (3.3) now becomes (3.6) The first term on the RHS is estimated as x using Hölder's inequality and Sobolev embedding.For the second term in (3.6), we interpolate the L 3 x norm, and apply the Poincaré, Hölder's, and Young's inequalities, as well as Sobolev embedding to get We also use the Poincaré inequality to convert the last term on the LHS of (3.6) 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.7) In the last inequality, we interpolated between the L p+2 x and L 2 x norms, which is valid for p > 2. For a sufficiently small κ, the second term on the RHS is absorbed into the LHS.On the other hand, when 1 ≤ p ≤ 2, the finite size of the domain implies L 2 x ⊆ L p 2 +1 x , which again leads to (3.7).Thus, for any p ≥ 1, (3.6) becomes x . (3.8) In order to show a decaying norm, we need coercive terms on the LHS, which have been achieved.
To control the RHS (particularly the second term), we derive a balance equation for a higher-order energy X(t), defined in (3.19).Combining E(t) with X(t) allows us to close the estimates.
3.3.Higher-order energy estimate.We now present more a priori bounds for ψ and u, involving one more derivative than the energy E. (3.9) The first term on the RHS of (c-NLS) leads to a term which vanishes due to the periodic boundary conditions.We now estimate the RHS of (3.9).For the first term, we have x .
The first term on the RHS acts as the dissipative term for ψ at this higher-order energy level.For I 4 , we integrate by parts and use Hölder's inequality to obtain x .
Thus, (3.9) becomes (3.10) The first term is bounded using the Poincaré inequality, Sobolev embedding, and Lebesgue interpolation as x .We applied Young's inequality to extract out dissipative terms in the last step.Again, κ denotes a small number whose value shall be fixed later on, and C κ is a constant whose value depends on κ and the system parameters.For the second term on the RHS of (3.10), we have x Finally, we apply the Poincaré inequality and Sobolev embedding to bound I 7 .This leads to Combining all these inequalities into (3.10), and absorbing κ D 3 ψ 2 L 2 x into the LHS, we end up with x ∥ψ∥ 4 x . (3.12) This constitutes the higher-order estimate for the wavefunction.Next, we combine this with corresponding estimates for the velocity.Here, P is the Leray projector, which projects a Hilbert space into its divergence-free subspace, thus removing any purely gradient terms.Next, we multiply (c-NSE') by ∂ t u and integrate over the domain.This leads to (3.13) Henceforth, we repeatedly use the fact that the density is bounded both above and below (m , respectively.Thus, for the first term, In going from the second line to the third, we use the Sobolev embeddings and Lebesgue interpolation.Finally, Young's inequality lets us extract the required dissipative term.For the second integral of (3.13), we have where the Bψ term is handled via interpolation and Young's inequality, while the term ∇ψ is bounded using Sobolev embedding.For the third integral, where the Bψ term is handled just like in I 9 .Finally, for the last term, we integrate by parts and use (c-CON), which results in We estimate the second term in (3.14) via interpolation and a Sobolev embedding.This gives x .Similarly, for the third term in (3.14), αλ Putting together the above estimates into (3.13),we end up with where C κ depends on κ as well as the system parameters.In order to obtain the higher-order velocity dissipation ∥∆u∥ 2 L 2 x , we multiply (c-NSE') by −θ∆u, with θ > 0 to be fixed shortly, and integrate over the domain.This gives For the first term, we have The second integral is manipulated just as I 8 , namely, x .The integral I 14 requires Sobolev embedding, the Poincaré inequality, and Lebesgue norm interpolation, and reads x .In a similar manner, we have x .The last integral in (3.16) requires, like I 12 , only Hölder's and Young's inequalities.This provides In the end, (3.16) becomes θν 2 ∥∆u∥ 2 x . (3.17) We now add (3.12), (3.15) and (3.17).Then, we note from the definition of Bψ that x , where the last three terms on the RHS are exactly I 5 , I 6 , and I 7 .Using sufficiently small values for θ and κ, we also absorb x and ∥∆u∥ 2 L 2 x on the RHS into the LHS.Finally, we are left with d dt ∥∆ψ∥ 2 x .
(3.18)This is the higher-order energy estimate.Using similar arguments to Section 3.3.3 of [JJK23], we can use the Grönwall inequality to control the higher-order energy and dissipation.The results are summarized in the following lemma.
Lemma 3.3 (Algebraic decay rate for energies).We label the higher-order energy as Then, the sum Z := X + E decays as where Z 0 := Z(0).Moreover, the time-integral of the corresponding dissipation terms is also bounded.Specifically, (3.21)In addition, by integrating the higher-order energy estimate over [t, 2t] (for t ≥ 1), we end up with a time-decaying estimate for the dissipation, given by The last inequality in (3.21) is valid because Z 0 , S 0 ≤ 1.
3.4.Maximal parabolic regularity for ψ.From the previous analysis, we have obtained ψ ∈ L 2 [0,T ] H 3 x .However, as pointed out in the discussion following Definition 2.9, we seek Bψ ∈ L 2 [0,T ] L ∞ x , which follows from ψ ∈ L 2 [0,T ] H 7 2 + x .In the 2D case [JJK23], this was achieved by taking advantage of the Sobolev embedding x and deriving a "highest-order energy estimate" for ψ.This approach does not work here due to the embedding H x , which would require higher-order estimates on u and ρ.Instead, we exploit the parabolic nature of (c-NLS) and apply the method of maximal regularity to gain the necessary control of ψ.Lemma 3.4 (Maximal regularity for ψ).For δ as defined in Theorem 2.3, and a sufficiently small δ 1 > 0, we have the maximal regularity bound uniformly in time T , where f : (R + ) 2 → R + is a continuous polynomial, with f (0, 0) = 0.
Proof.We begin by rewriting (c-NLS) in a parabolic form as The differential operator in this case is A = − λ+i 2 ∆.Comparing this to (2.14), we see that m = 1, and a α = − λ+i 2 when α ∈ {(2, 0, 0), (0, 2, 0), (0, 0, 2)} and a α = 0 otherwise.Thus, the first condition in Definition 2.10 is satisfied.The principal symbol of the operator is from which it is clear that Ã(ξ) −1 is also bounded for |ξ| = 1.Finally, the spectrum σ λ+i 2 |ξ| 2 belongs to the sector Σ ζ 0 for tan ζ 0 > λ −1 > 0. Thus, the operator A is uniformly ( , 2 √ 1+λ 2 }, and the maximal parabolic regularity estimate is applicable.We now act upon (3.24) by (−∆) 3 4 + δ 1 2 for some δ 1 > 0 that will be determined shortly.We then apply Lemma 2.11 with X = L 2 (T 3 ) and X 1 = H 2 (T 3 ), which results in for r = 1 + δ with δ ∈ (0, 1 3 ), and this restriction will become clear when estimating I 18 in (3.27).It is important to note that L r t is calculated on [0, T ], where T > 0 is the local existence time defined in (2.9).We now estimate each of the terms on the RHS.The norm of the initial condition is found to belong to a Besov space as a result of the interpolation (see [AF03, Chapter 7]).Indeed, we have for sufficiently small values for δ, δ 1 , and δ 2 .The first inequality is due to the embedding H s+δ 2 ⊂ B s 2,q for any δ 2 > 0 (see [Lu21, Lemma 2.2]).Due to the restriction δ < 1 3 , we have 2δ 1+δ < 1 2 , implying that we may bound I 17 with the initial data, as in the last inequality of (3.26).Next, we deal with the second term on the RHS of (3.25) as The second inequality follows from the product rule for Sobolev norms [MB02, Lemma 3.4], and the third inequality from Agmon's inequality and Sobolev embedding.The fourth inequality is due to Hölder's inequality, while the final step follows from (3.20) and (3.21).In a similar way, we can analyze the third term on the RHS of (3.25), yielding , uniformly in T , and that their norms can be made small by an appropriate choice of S 0 and Z 0 .
3.5.Ensuring global-in-time positive density.We have now obtained all the a priori estimates needed to return to (2.13).Proof of Theorem 2.3: Using (CPL), we have (3.31) We now show that each of these four terms can be made as small as required, by choosing sufficiently small values of the data, i.e., S 0 and Z 0 .For the first term (I 21 ), we use Hölder's inequality, (3.1), and (3.20) to write (3.32) The above calculation assumes that p < 1 + 1 δ (which is consistent with the definition of δ in Theorem 2.3), so that a uniform-in-T bound can be obtained.From Lemma 3.4, we conclude the smallness (in terms of the initial data) of the last factor in the RHS, i.e., D 2 ψ L 1+δ [0,T ] H 3 2 +δ 1 x .Thus, I 21 is independent of T , and can be made sufficiently small by appropriate initial data.
Moving on to the remaining terms in (3.31), we have x . (3.33) All the terms are, once again, bounded in terms of the data according to (3.20) and (3.21).In the same way, we have where all the terms are controlled by the data.Finally, using (3.1) and (3.20), we arrive at and this calculation holds for all values of p ≥ 1.From the analysis in (3.32)-(3.35),we conclude that the LHS of (3.31) can be made sufficiently small to satisfy the constraint given by (2.13), for all values of time T > 0. This implies that the density always remains bounded below, and thus solutions are global in time for p < 1 + 1 δ .No matter how large (but finite) p is, it is possible to choose δ small enough so that the constraint in (2.13) is met.□ and ξ ∈ R n with |ξ| = 1.Here, σ(B) refers to the spectrum of the operator B, and Σ ζ := {z ∈ C : |arg z| ≤ ζ} is a sector in the right half of the complex plane.