Global mass-preserving solutions for a two-dimensional chemotaxis system with rotational flux components coupled with a full Navier-Stokes equation

We study the chemotaxis-Navier-Stokes system \[\left\{\; \begin{aligned} n_t + u\cdot\nabla n&=\Delta n - \nabla\cdot (nS(x,n,c)\nabla c),&&x\in\Omega, t>0, \\ c_t + u\cdot\nabla c&=\Delta c - n f(c),&&x\in \Omega, t>0, \\ u_t + (u\cdot\nabla) u&=\Delta u + \nabla P + n \nabla \phi, \;\;\;\; \nabla\cdot u = 0,&&x\in\Omega, t>0 \end{aligned}\right.\tag{$\star$} \] with no-flux boundary conditions for $n$, $c$ in a bounded, convex domain $\Omega\subseteq\mathbb{R}^2$ with a smooth boundary, which is motivated by recent modeling approaches from biology for aerobic bacteria suspended in a sessile water drop. We further do not assume the chemotactic sensitivity $S$ to be scalar as is common, but to be able to attain values in $\mathbb{R}^{2\times 2}$, which allows for more complex modeling of bacterial behavior near the boundary. This is seen as a potential source of the structure formation observed in experiments. While there have been various results for scalar chemotactic sensitivities $S$ due to a convenient energy-type inequality and some for the non-scalar case with only a Stokes fluid equation (or other strong restrictions) simplifying the analysis of the third equation in ($\star$) significantly, we consider the full combined case under little restrictions for the system giving us very little to go on in terms of a priori estimates. We nonetheless manage to still achieve sufficient estimates using Trudinger-Moser type inequalities to extend the existence results seen in a recent work by Winkler for the Stokes case with non-scalar $S$ to the full Navier-Stokes case. Namely, we construct a similar global mass-preserving generalized solution for ($\star$) in planar convex domains for sufficiently smooth parameter functions $S$, $f$ and $\phi$ and only under the fairly weak assumptions that $S$ is appropriately bounded, $f$ is non-negative and $f(0) = 0$.


Introduction
We will study a mathematical model from biology describing chemotaxis, the directed movement of cells along a chemical gradient towards an attractant, and its interaction with the fluid said cells are suspended in. While pure chemotaxis models, which generally consist of two coupled partial differential equations, without added fluid interaction have been quite heavily studied for some time now (initiated by the seminal work of Keller and Segel in [9]), recent findings by Dombrowski et al. (cf. [7]) have shown that in populations of Bacillus subtilis much stronger liquid movement than could originate from independent bacteria can be observed after cell aggregates have been formed. This calls into question the prior modeling assumption that the liquid-cell interaction can basically be disregarded because each single organism has negligible influence on the liquid. This observation then led Tuval et al. (cf. [15]) to an extended modeling approach, which adds a full Navier-Stokes type fluid equation to the standard chemotaxis model. The central new interactions are that the attractant and the cells suspended in the solution are affected by convective forces exerted by the fluid, while the fluid is affected by buoyant forces exerted by the cells because of their comparatively large size.
In mathematical terms and after some normalization of parameters, this modeling approach by Tuval et al. leads to a system of coupled partial differential equations of the following form: n t + u · ∇n = ∆n − ∇ · (nS(x, n, c)∇c), x ∈ Ω, t > 0, c t + u · ∇c = ∆c − nf (c), x ∈ Ω, t > 0, u t + (u · ∇)u = ∆u + ∇P + n∇φ, ∇ · u = 0, x ∈ Ω, t > 0. (1.1) Here Ω ⊆ R N with N ∈ N is some space domain and the functions n = n(x, t), c = c(x, t), u = u(x, t), P = P (x, t) model the bacteria population, the attractant (e.g. oxygen) concentration and the fluid velocity field and associated pressure respectively. The system is further parameterized by the given functions S, which models the chemotactic sensitivity of the bacteria, f , which models the attractant consumption rate, and φ, which represents the gravitational potential. Convection is modeled by the terms u · ∇n and u · ∇c, while buoyant forces are modeled by the term n∇φ.
This system is fairly well understood if we assume the chemotactic sensitivity S to be a scalar function as this brings the system quite close to the classical Keller-Segel model in terms of the chemotactic interaction (cf. [9]). Amongst many other things, global well-posedness of this system in two or three dimensions with different assumptions on the initial data and parameter functions f , S, φ has been extensively studied (cf. e.g. [4], [8], [19], [11]). In some two dimensional settings, there are global, unique classical existence results available (cf. [19]). While in three dimensional cases there are still existence results available, they are somewhat less ambitious and deal with weaker notions of solutions and only consider eventual smoothness properties (cf. [19], [22], [23]). For a more broad survey of mathematical chemotaxis models and recent results about them, consult for instance [1].
For the functions f, S and φ that parametrize (1.1), we will throughout this paper assume that and that, for S = (S ij ) i,j∈{1,2} , and that and that φ ∈ W 2,∞ (Ω Complications: Loss of energy structure and nonlinear convection. The above setting presents us with two key complications, one inherited from the Stokes case already discussed in [24] and one we reintroduce from the original model by Tuval et al. The complication we inherited is of course the non-scalar chemotactic sensitivity S. Such S are especially interesting from a biological point of view because scalar S have been shown to lead to long time homogenization for some common parameters, which does not agree with the structure formation seen in experiments (cf. [20]). As these observations further suggest that spatial inhomogeneities often originate at the boundary (cf. [7]), modern modeling approaches introduce rotational flux components near said boundary, which necessitate sensitivity functions that look somewhat like with significant non-diagonal entries near the boundary (cf. [26], [27]). The most devastating consequence of this for the analysis of the model is that this leads to the apparent loss of an energy type functional, which is present in the scalar case and often key to proving global existence of solutions (cf. [8], [28], [19], [13]) or understanding their long time behavior (cf. [20], [23]). Results dealing with non-scalar sensitivities therefore often analyze the system under some very strong restriction on either S or the initial data (cf. [3], [16], [17]).
Even in the fluid-free version of (1.1), this loss of structure has led to a significant lack of knowledge and, e.g. for N = 2, global smooth solutions have thus far only been constructed under significant smallness conditions for c 0 (cf. [10] or [2] for an extension to the fluid case). For arbitrarily large data and dimension N , it is possible to construct global generalized solutions (not unlike those constructed in this paper) as seen in [21] at the very least.
The complication we reintroduced is adding the nonlinear convection term, which is disregarded in [24], back into the third equation making it a full Navier-Stokes fluid equation. This means that especially the semigroup methods used in the Stokes case lose some of their effectiveness reducing our immediate access to regularity information for the third equation and making it in general much tougher to handle.

Main ideas.
A lot of ideas from the Stokes case [24] translate fairly immediately for the first two equations in (1.1), especially concerning the handling of the somewhat problematic chemotactic sensitivity S as long as we still manage to provide similar bounds for u as in the reference. Therefore it is in establishing these bounds where our key ideas come in. As already mentioned, semigroup techniques lose some of their fruitfulness in the third equation due to the newly introduced nonlinear convection term, but similarly to the Stokes case, the rather weak regularity information of the form ∞ 0 Ω |∇n| 2 (n + 1) 2 ≤ C in Lemma 3.1 is of central importance yet again (In the Stokes case, it was mostly used to tease out integrability properties of the time derivative of ln(n + 1) and for some additional compactness properties). It in fact allows us in conjunction with some functional inequalities derived from the Trudinger-Moser inequality seen in Lemma 3.2 to conclude that terms of the form Ω (n + 1) ln n + 1 n + 1 and Ω n∇φ · u have similar time integrability properties. Time integrability of the former term proves useful to simplify some compactness results for n that have already been shown in [24], while integrability of the latter term will be the central keystone to showing sufficient L 2 type bounds for u and ∇u (cf. Lemma 3.4) by testing the third equation with u itself. While not as strong a set of bounds for u as in the Stokes case, this then proves to be enough for the construction of global mass-preserving generalized solutions in the sense of Definition 2.1 below via a similar approximation approach as the one seen in [24].

Generalized solution concept and approximate solutions
As regularity information is rather hard to come by due to the complications outlined in the introduction, we will not endeavor to construct a global classical solution for (1.1), but instead confine ourselves to a very generalized notion of solution. Because of the similarities to the pure Stokes case in [24] and our desire to not unnecessarily duplicate effort, we let ourselves be guided by the generalized solution concept seen in said reference, which we of course slightly adapt to the full Navier-Stokes case. This reads as follows: We call a triple of functions a global mass-preserving generalized solution of (1.  Furthermore at least in the Stokes case, some similar generalized solutions have been shown to eventually (from some large t > 0 onwards) attain such a level of regularity (cf. [25]). Similar to Reference [24], the key to our existence results will be approximate problems defined in the following way: We first fix families of functions (ρ ε ) ε∈(0,1) and (χ ε ) ε∈(0,1) with ρ ε ∈ C ∞ 0 (Ω) such that 0 ≤ ρ ε ≤ 1 in Ω and ρ ε ր 1 pointwise in Ω as ε ց 0 and constructed by standard methods. For ε ∈ (0, 1), we then define and consider the following initial boundary value problem: x ∈ Ω, t > 0, x ∈ Ω, t > 0, ∇n ε · ν = ∇c ε · ν = 0, u ε = 0, x ∈ ∂Ω, t > 0, n ε (x, 0) = n 0 (x), c ε (x, 0) = c 0 (x), u ε (x, 0) = u 0 (x), x ∈ Ω. (2.7) This kind of regularized version of (1.1) with (1.2) and (1.3) then easily admits a global classical solution because it substitutes standard Neumann boundary conditions for the more complex no-flux boundary condition seen in (1.2) and makes the first equation accessible to comparison arguments with a non-zero constant to gain a global upper bound for n ε , which would be much harder to achieve otherwise. This of course only works under similar assumptions on the parameter functions f , S, φ as proposed in the introduction. As the techniques to achieve an existence result for the approximated system above are fairly well-documented and do not appreciably differ for our case in comparison to e.g. the case studied in [19], we will only give a short sketch to justify the following existence theorem: Proof. Standard contraction mapping approaches in an appropriate setting (as e.g. seen in [19, Lemma 2.1] for a similar system) provide us with a classical solution for (2.7) on a space time cylinder Ω × [0, T max,ε ) with some maximal T max,ε ∈ (0, ∞] and a blow-up criterion of the following type: Here q is some number greater than 2. Non-negativity and positivity on Ω × (0, ∞) of c ε and n ε respectively are immediately ensured by maximum principle. Note now further that, because we defined the S ε to be zero for sufficiently large values of n, well-known comparison arguments can be used to already gain a global upper bound for n ε on the whole cylinder. This already rules out blowup regarding n ε . As the second equation in (2.7) is generally fairly unproblematic, similar boundedness results can be achieved for c ε (cf. Lemma 3.1 later in this paper). Regarding the possible blowup of c ε or u ε , we can look at the prior work done in Section 4.2 of Reference [19], where it is proven for much weaker prerequisites as already established here that the two norms in the blowup criteria regarding c ε and u ε respectively are bounded as well. Note that this is mostly done using the second and third equation of the system studied in said reference, which are apart from some slight generalizations the same as the second and third equation in (2.7). Only one step in the reference uses an energy inequality not available to us to establish a bound for Ω |∇c ε | 2 , which is also easily gained by a straightforward testing procedure for the second equation in (2.7) without using said energy inequality.

Results for n ε and c ε reusable from the Stokes case
Let us now start by revisiting some key results for the approximate solutions n ε , c ε , which can be derived in a very similar way to the Stokes case (cf. [24]) as they stem from just considering the first two equations in (2.7).
which for fitting choices of p yields (3.2) for finite p as well as (3.3). The case p = ∞ in (3.2) then follows via the limit process p → ∞. To now derive (3.4), we test the first equation in (2.7) with 1 nε+1 to obtain which we then further improve to 1 2 Ω |∇n ε | 2 (n ε + 1) 2 ≤

Estimates for u ε based on the Trudinger-Moser inequality
While semigroup methods proved very fruitful when the third equation in (2.7) is of Stokes type, they are in our case thoroughly thwarted by the nonlinear convection term (u ε · ∇)u ε . As such, we will use very different tools to at least partially recover some of the L p boundedness results for u ε and its derivatives seen in Lemma 3.3 of Reference [24] for the Stokes case. To fill the role of these tools, we therefore start by deriving two functional inequalities based on the Trudinger-Moser inequality (pioneered in [14], [12]) and inspired by [18]: and with ψ := 1 |Ω| Ω ψ for all ϕ ∈ C 1 (Ω), positive ψ ∈ C 1 (Ω) and a > 0.
Proof. By using the Trudinger-Moser inequality seen in Proposition 2.3 of Reference [5], which is applicable because Ω is convex and therefore finitely connected, we gain a constant K 1 ≥ |Ω| with for all ξ ∈ C 1 (Ω) with Ω ξ = 0 and Ω |∇ξ| 2 ≤ 1. Note that we can choose the constant β seen in Proposition 2.3 from the reference equal to 2π because Ω has a smooth boundary. Using Young's inequality to see that for all ξ ∈ C 1 (Ω). Now fix ϕ ∈ C 1 (Ω), positive ψ ∈ C 1 (Ω) and a > 0. We then observe that the estimate holds with m := Ω ψ because of Jensen's inequality. If we now apply (3.8) with ξ := aϕ to this and multiply by m a , we get that after some rearranging. This gives us our first result. Now only fix a positive ψ ∈ C 1 (Ω). We can then choose ϕ := ln ψ ψ and a := 2 in (3.9) to get that Because by Jensen's inequality Ω ln ψ ψ ≤ 0, this directly implies our second result.
Employing the second functional inequality (3.6), we can now use the fairly weak regularity information in (3.4) to derive the following preliminary integrability property for the family (n ε ) ε∈(0,1) , which will not only prove useful to derive bounds for the family (u ε ) ε∈(0,1) , but also to later simplify a compactness argument, which was already used in the Stokes case. (n ε + 1) ln n ε + 1 for all ε ∈ (0, 1).
This in combination with another application of our functional inequalities above then serves as the basis to prove the central result of this section, namely the L 2 type bounds for the functions (u ε ) ε∈(0,1) and their gradients, which replace the results gained via semigroup arguments in the Stokes case: for all ε ∈ (0, 1).
Proof. As our first step, we test the third equation in (2.7) with u ε itself to gain that and ε ∈ (0, 1).
We now apply (3.5) from Lemma 3.2 (with ϕ := ∇φ · u ε and ψ := n ε + 1) to the rightmost term in the previous inequality to gain a constant K 1 > 0 such that for any a > 0 and each t ∈ [0, T ], ε ∈ (0, 1). Further note that Here H φ denotes the Hessian of φ and C p is the Poincaré constant for Ω. If we now apply this to (3.14) and set for all t ∈ [0, T ] and ε ∈ (0, 1). Because of Lemma 3.3, there further exists a constant K 4 (T ) > 0 such that for all ε ∈ (0, 1). If we now integrate (3.15), this property of g ε then directly gives us that for all t ∈ [0, T ], which immediately implies (3.12) and (3.13).

Construction of limit functions
Having now established weaker (though still suitable) uniform bounds for u ε than those seen in the Stokes case, we will make the final preparations for the construction of limit function for our family of approximate solutions as ε ց 0. We do this by proving some additional, albeit fairly weak, boundedness results for the time derivatives of all three families (n ε ) ε∈(0,1) , (c ε ) ε∈(0,1) and (u ε ) ε∈(0,1) as this will provide the last prerequisite for some key applications of the Aubin-Lions lemma (cf. [13]): for all ε ∈ (0, 1).
Proof. We will only give this proof in full detail for (3.18) and then just provide a sketch for (3.16) and (3.17) as all three inequalities can be proven in quite similar a fashion and more detailed proofs for (3.16) and (3.17) can be found in [24].
This then allows us to use essentially three applications of the Aubin-Lions lemma (cf. [13]) to prove the following sequence selection and convergence result: There exists a sequence (ε j ) j∈N ∈ (0, 1) with ε j ց 0 as j → ∞ such that 2 and a.e. in Ω × (0, ∞) and as ε = ε j ց 0 and a triple of limit functions (n, c, u) defined on Ω × (0, ∞) and satisfying n, c ≥ 0 and ∇ · u = 0 almost everywhere.
Note that the non-negativity properties directly transfer from the approximate functions because of the pointwise convergence, while ∇·u = 0 is directly ensured by u being an element of L 2 loc ([0, ∞); W 1,2 0,σ (Ω)). and use of the boundary conditions in (2.7). We therefore only need to further argue that these properties survive the necessary limit process. For most terms in the integral equality (2.6) concerning u, this is fairly straightforward to show using the convergence properties established in Lemma 3.6, but we nonetheless give the full argument for at least the newly introduced term (compared to the Stokes case) as an example. This is especially pertinent as we needed to establish stronger convergence properties for u to handle this term compared to [24], namely strong L 2 as opposed to weak L 2 convergence.