On the Field-Induced Transport of Magnetic Nanoparticles in Incompressible Flow: Existence of Global Solutions

We prove global-in-time existence of weak solutions to a pde-model for the motion of dilute superparamagnetic nanoparticles in fluids influenced by quasi-stationary magnetic fields. This model has recently been derived in Grün and Weiß(On the field-induced transport of magnetic nanoparticles in incompressible flow: modeling and numerics, Mathematical Models and Methods in the Applied Sciences, in press). It couples evolution equations for particle density and magnetization to the hydrodynamic and magnetostatic equations. Suggested by physical arguments, we consider no-flux-type boundary conditions for the magnetization equation which entails H(div,curl)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H({\text {div}},{\text {curl}})$$\end{document}-regularity for magnetization and magnetic field. By a subtle approximation procedure, we nevertheless succeed to give a meaning to the Kelvin force (m·∇)h\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\mathbf {m}\cdot \nabla )\mathbf {h}$$\end{document} and to establish existence of solutions in the sense of distributions in two space dimensions. For the three-dimensional case, we suggest two regularizations of the system which each guarantee existence of solutions, too.

Here, (u, p) denote the hydrodynamic variables of the carrier fluid, c stands for the number density of the magnetic nanoparticles and m or h = ∇R describe magnetization or magnetic field, respectively. The vector field V part denotes the particle's velocity relative to the flow of the carrier fluid. It takes diffusive and magnetic effects into account. The system is driven by the external magnetic field h a which satisfies Maxwell's equations in the absence of matter, i.e. For further explanation of the parameters in the model, we refer the reader to Sect. 2. The model includes a two-domain-approach, i.e. the magnetic field is defined on a larger domain Ω ⊃ Ω, which has the advantage that one can account for stray field effects at the boundary of the fluid domain.
As already pointed out in [17], in the literature so far two pathways have been pursued to study problems of ferrohydrodynamics. The first one is concerned with phenomena for which the particle distribution can be assumed to be homogeneous in space (and consequently constant in time). For those ferrofluids, pde-models have been derived by Shliomis [25] and Rosensweig [24]. Both models couple evolution equations for momentum and magnetization to Maxwell's equations or to their simplifications from magnetostatics. Rosensweig takes in addition an evolution equation for the angular momentum of the fluid into account. In a series of papers, Amirat and Hamdache [1][2][3][4][5][6] developed a mathematical existence theory in the framework of the Shliomis model. Just recently, Nochetto, Salgado and Tomas [20] proposed numerical schemes for the Rosensweig model. In a second paper, they [19] considered two-phase flow with one ferrofluid involved to model the famous Rosensweig instability [23], and they provided a convergence proof in a simplified setting. All these publications have in common that there are no pathways suggested how to deal with non-homogeneous, non-steady particle densities.
In a second line of research, mathematical models have been suggested and investigated for the transport of magnetic nanoparticles with particle densities varying in space and time. These models have in common that evolution equations for the magnetization are not considered separately. Instead, authors assume the magnetization to be given explicitly as a function of particle density and magnetic field. Polevikov and Tobiska [22] were interested in a steady-state diffusion problem for particles in a ferrofluid. Most recently, Himmelsbach, Neuss-Radu and Neuß [18] proposed a new model featuring an evolution equation for the particle density coupled to the magnetostatic equations and assuming the macroscopic flow field to be given. In the radial symmetric case they show existence and uniqueness of solutions and provide numerical simulations, too.
In this spirit, model (1.1) is the first model which takes both non-constant particle densities and magnetization fields into account. Note that the boundary conditions (1.2c), (1.2d) are of no-flux-type, combining the no-flux boundary condition (1.2b) of the number density c with the addditional flux originating from the diffusive term −σΔm = −σ (∇ div m − curl curl m). These boundary conditions may be favorable from a physical point of view. They allow for energy estimates in a natural way and they do not prohibit tangential or normal traces of m to be different from zero on ∂Ω. This is in contrast to the work of [1] where the normal component of the magnetization has been prescribed to vanish and therefore H 1 -regularity holds for m and h globally.
In our framework, however, the problem arises that magnetization and magnetic field have only H(div, curl)-regularity, cf. [8]. This is the more an issue, as the Kelvin force (m · ∇)h requires control of gradients of h and as it enters all the evolution equations (1.1a), (1.1c), and (1.1f). In this situation, it seems natural to derive H 1 loc -regularity for m and h and to consider solutions in the sense of distributions. Due to the intricate coupling of the evolution equations for m and c, estimates on gradients of (appropriate powers of) c depend on the integrability of (m · ∇)h and of (m · ∇)m. As a consequence, we expect results only in L p loc -spaces. As such a localization seems to become indispensable already on the level of Galerkin approximations, the appropriate choice of approximation spaces is a central topic of this paper. This includes in particular a strategy how to deal with the two-domain modeling approach.
The outline of the paper is as follows. In Sect. 2, we explain further features of the model (1.1), we state an energy estimate, and we formulate our hypotheses on the domains Ω and Ω and on the data.
To cope with the nonlinearities (m·∇)h and (∇m) T m in the convective term of the evolution equation for the particle density c, we refrain to local arguments -which are needed already in the first limit passage Discrete to Regularized Continuous. This requires a subtle choice of Galerkin approximation spaces as in general the projections Π Xn ϕ, n ∈ N, of C ∞ 0 -functions onto ansatz spaces X n do not have compact support. For this reason, some techniques which we use for the passage to the limit Regularized Continuous to Continuous, cannot be applied at this earlier stage. The remedy is to choose the approximation spaces in such a way that L ∞ -convergence of the gradients of Π Xn ϕ to ∇ϕ is guaranteed. 1 . A sufficient condition is to require H 3 -convergence of Π Xn ϕ to ϕ. Inspired by the discussion of boundary conditions for the magnetization fields in [20], we prefer to take natural boundary conditions div m ∂Ω ≡ 0 and curl m ∂Ω ≡ 0 for our ansatz spaces. It is worth mentioning that the boundary conditions (1.2c), (1.2d) reduce to these conditions in those points (x, t) ∈ ∂Ω × (0, T ) where c(x, t) = 0. In addition, some effort is devoted to guarantee that ansatz functions for R are defined on Ω , having in addition gradients, the restriction of which to Ω is contained in the approximation space for m. We devote Sect. 3 to construction and decomposition of such ansatz spaces for magnetization and magnetic potential. Then, Sect. 4 collects the ansatz spaces for velocity and particle density.
In Sect. 5, we introduce a (T)ransport and (M)obility (R)egularized model -replacing the usual c log centropy by a strictly convex approximation with quadratic growth, and using a further density cut-off in the transport velocity V part , see (5.2) and (5.3). For this TMR-model, global existence of discrete solutions is established, and Sections 6 and 7 provide compactness results as well as the limit passage Discrete to Regularized Continuous.
In Sect. 8, we show that solutions to the TMR-model converge in the limit of vanishing regularization parameters to weak solutions of model (1.1). For this, it will be essential to choose the nonlinear mobility in the particle evolution in such a way that the flux V part has L 2 -regularity. It turns out that this can be achieved by choosing the mobility quadratic in c. Due to a certain regularity gap, we have to confine ourselves to the case of two space dimensions unless additional regularizing terms are considered. For this, see Remark 8.11. In this paper, we cannot avoid a rather involved notation. For the reader's convenience, the Appendix A.5 explains the notation and provides references on definitions and further properties.
For the reader's convenience, we include a formal derivation of the estimate in the appendix, see Sect. A.1, together with some general reflections on the physical background of the model. Let us formulate our assumptions on the spatial domains and on the data.
For the reader's convenience, we give a definition for two-dimensional curl-operator and cross product for vector fields or vectors, respectively. Definition 2.1. Let u : Ω → R 2 be such that the weak curl of (u x , u y , 0) T exists, then-without introducing new notation-the curl-operator of two dimensional vector fields is defined by which is a scalar function. Analogously, the vector product × will be defined for vectors a, b ∈ R 2 and scalar g ∈ R as follows.
We have In this section, we introduce the function spaces which will serve for the construction of approximation spaces for magnetization and magnetic field in the Faedo-Galerkin approach of Sect. 5. Our choice is guided by the following criteria. (C1) Formal energy estimates, compare (2.2), indicate that the magnetization is contained in L 2 ((0, T ); H(div, curl)(Ω)). (C2) The magnetic field h is a gradient field on Ω , satisfying h| Ω \Ω ∈ H(div 0 )(Ω \Ω) due to the magnetostatic equations curl h = 0, div(h + m) = 0, and h| Ω ∈ H(div)(Ω) due to the formal energy estimate (2.2). (C3) Approximation functions for the magnetization should satisfy the boundary conditions which allow for stability estimates. (C4) We require the sequence of approximation spaces for the magnetization to be dense in a closed subspace of H 3 (Ω) d . This leads us to consider the space as the function space related to the magnetization and the space as the function space for the scalar potentials of the magnetic field h. Note that • M is a dense subset of L 2 (Ω) d as it contains all C ∞ 0 (Ω) d -functions. • M is a closed subspace of H 3 (Ω) d as it is the preimage of closed sets under continuous mappings.
• R is complete with respect to the norm The latter is true as the norm · R dominates · H 1 (Ω ) (by using Poincaré's inequality for mean value free functions) and Δ(·) L 2 (Ω) . Hence, if {r k } k∈N is a Cauchy sequence in R, there are functions r ∈ H 1 mean (Ω ) and s ∈ L 2 (Ω) such that r k → r in H 1 mean (Ω ) and Δr k → s in L 2 (Ω). The identification s = Δr is done in a standard way via integration by parts. The fact that Δr vanishes on Ω \Ω follows analogously.

Construction of a Basis of M
The goal of this subsection is to construct a basis of M which is orthonormal with respect to the L 2scalar product. In particular, we aim at substructures which can be exploited to be the starting point to construct a basis of the set R as well. For consistency reasons with respect to the Galerkin procedure to be studied later on, we require that approximation spaces R n ⊂ R and M n ⊂ M satisfy ∇R n | Ω ⊂ M n .
Let us first introduce two subspaces of M.
Proof. By the completeness of H 1 0 (Ω) and Poincaré's inequality, the L 2 -closure of S is a subset of ∇[H 1 0 (Ω)]. Hence, its intersection with H is a subset of S by definition, see (3.1.2). The converse inclusion is obvious. By a similar reasoning, the L 2 -closure of H is a subset of ∇[H 1 mean (Ω)] and its intersection with M guarantees H 3 -regularity as well as vanishing of the divergence on the boundary. Again, the converse inclusion is obvious. Now, we formulate a general decomposition lemma, which might be of independent interest, and which yields the desired basis with a gradient substructure. We postpone its proof to the appendix. Lemma 3.1.3. Let U, X, Y be separable Hilbert spaces that satisfy the following assumptions.
• U ⊂ X is a closed subspace, and we have Then, the following holds true: and {u i } i∈N1 , {v i } i∈N2 form orthonormal sets in Y and orthogonal sets in X, respectively. (iii) We have the decomposition and for each i ∈ N 1 and j ∈ N 2 , we have In particular, the family Proof. By Lemma 3.1.3-applied to S ⊂ H-we have bases that are L 2 -orthonormal and orthogonal with respect to the inner product of For later use, we rename and relabel the basis functions of M.
By construction we have the following properties, with respect to the L 2 and H 3 scalar products, Ψ m 4i−2 ⊥Ψ m 4j−2 for i, j ∈ N, i = j, with respect to the L 2 and H 3 scalar products, Ψ m 4i ⊥Ψ m 4j for i, j ∈ N, i = j, with respect to the L 2 and H 3 scalar products, Ψ m i ⊥Ψ m j for i, j ∈ N, i = j, with respect to the L 2 scalar product. (3.1.14)

Construction of a Basis of R
Recalling criterion (C2) as well as our requirement that approximation spaces R n and M n should satisfy ∇R n | Ω ⊂ M n , it is natural to extend in a first step functions in H ⊂ A 0 to the whole of Ω , consistent with the definition of R. For this, we choose the uniquely determined mean-value-free potentials To explain our extension procedure, we begin with some considerations valid on a general bounded C 1,1domain V . Note that these ideas later on shall be applied both to Ω \Ω and to Ω.
be defined for f ∈ H 1/2 (∂V ) to be the unique solution of the inhomogeneous Dirichlet-Laplace problem Also, denote by {u V i } i∈N the eigenfunctions to positive eigenvalues (λ V i ) i∈N of the homogeneous Neumann-Laplace problem, (i) Using Definition 3.2.1, we can extend and augment the basis {h i } i∈N . Definẽ and normalize them by Note that by this choice we have R 4i | Ω \Ω ≡ const. for all i ∈ N. The reason for this is the fact that h 2i = s i , cf. (3.1.11), admits a trace-free potential, hence its mean-value-free potential is constant on ∂Ω. This constant, see c 2i above, has been subtracted from φ Ω 2i and therefore on Ω \Ω the potentialR 4i is equal to L −1 Ω \Ω 0, which is constant zero. As a first step to find suitable approximation spaces for the magnetic potential we choose (3.2.7) Let us show S to be constant. For this, we consider three types of testfunctions. We start with the following. Let {ψ dir i } i∈N be the eigenfunctions associated with positive eigenvalues (μ i ) i∈N of the homogenous Dirichlet Laplace problem, Due to (H1), those functions are H 4 -regular. Then define, for all i ∈ N, , which implies ∇S ∈ H(div 0 , curl 0 )(Ω). And consequently, we get for all i ∈ N 0 , using the functions (3.2.11), (∂Ω), this implies that ∇S is in H(div)(Ω ) globally and therefore ∇S ∈ H(div 0 , curl 0 )(Ω ). By the last class of functions we similarly get for all i ∈ N 0 . We finally conclude (on simply connected Ω with C 1,1 -boundary [8]) that ∇S ∈ H n0 (div 0 , curl 0 )(Ω ) = {0} and therefore S ∈ R is constant with zero mean, i.e. S = 0.
In Sect. 5 we will use Galerkin approximate solutions to obtain a solution to our model. The usual approach is to define a space given as the linear hull of only finitely many elements of our complete sets for M and R. For this, it is crucial that those sets are linearly independent. The generating set of M already is a basis. In the case of R linear independency is not evident. Moreover, for later purposes, we want to orthogonalize parts of the already established generating set {R i } i∈N . Those issues are addressed in the following lemma.

Lemma 3.2.4. There exists a basis {ψ
which satisfies the following.
The basis is a linearly independent selection of functions from {R i } i∈N . Therefore, by construction we will have i) and ii), see (3.2.1) and Lemma 3.1.4.
Consider the set {R 2i } i=1,...,n . This set is linearly independent as it is linearly independent on Ω due to ∇R 2i | Ω = ∇φ Ω i = h i . From the set {R 2i−1 } i∈N we can pick n elements by induction. We start with R i1 := R 1 .
It remains to find functions ψ R 2i−1 , i ∈ N, such that iii) is satisfied. By applying Gram-Schmidt-orthogonalization to the functions {ψ R 2i−1 } i∈N , we get new functions (3.2.14) Due to Poincaré's inequality the bilinear form above is an inner product and our procedure is well-posed.

Lemma 3.2.5. The basis {ψ
Proof. The application of Gram-Schmidt-orthogonalization in the proof of Lemma 3.2.4 only orthogonalizes the basis functions with odd indices, therefore the properties of ∇R 2i | Ω = h i do not change. Also, In conclusion, the claim follows.

Construction of Discrete Spaces for Velocity Field and Particle Density
In this section, we briefly specify our choice of Galerkin functions to approximate particle density and flow field. Observing that [15,Theorem 3.4]) and a closed subset of For the particle density, we follow the usual approach and take the complete set {ψ c i } i∈N of eigenfunctions of the Laplacian on Ω subjected to homogeneous Neumann boundary conditions. By standard results, they are H 2 -regular, H 1 -orthogonal and form a basis of . (4.7)

Existence of Discrete Solutions
In this section, we prove the existence of global solutions to a discretization of an appropriate transport and mobility regularized version of model (1.1). We call this regularization the TMR-model. It differs from model (1.1) by a new viscosity term in (1.1c) and a cut-off near zero applied to c-terms in the denominator of V part . In addition, we use a regularized entropy g L s as well. It reads Obviously, Here, The boundary condition (1.2b) is adapted accordingly, i.e. cV part ·ν| ∂Ω = σ c ∇c·ν| ∂Ω . In our weak solution concept the pressure vanishes by the use of solenoidal testfunctions in the Navier-Stokes equations. We will make use of the space Also recall the definitions of M and U from (3.3) and (4.1). We definẽ M to be identical to M but equipped with the norm of the sum which is the sum of all H 3 -norms of the individual summands.
hold and we adapt the corresponding result for the closure in H 1 (Ω) d which can be found in [12, p. 149]. The only difference is that (in the notation of [12]) we need an H 3 -estimate and an compactness result for solutions . This is a consequence of Theorem 3.2 in [12].

Galerkin approximation
In the following we fix σ c > 0 and 0 < s < e < L < ∞, where e is Euler's number.

Definition 5.3. Let
Remark 5.4. (i) The projection Π Hn , defined in (5.8), is well-defined, as can be seen by the following facts. A function ψ ∈ R can be written in terms of the basis (3.2.12) with convergence of infinite sums in the norm of R, which dominates the H 1 (Ω )-norm. Hence, ∇ and the infinite sum commute. Also restriction to Ω and the infinite sum commute. The projection is orthogonal due to L 2 -orthogonality, cf. (3.1.14), of the basis in (3.1.13). (ii) The projections Π Cn , Π Un , Π Mn are (at least L 2 -)orthogonal projections, see (4.5), (4.6) and (4.2), (4.3) and (3.1.14).
Obviously, analogous statements as in iii) hold for the spaces related to C n and M n , too. Concerning the space R n the argument is more involved. One needs to distinguish between two cases. On Ω, L 2 -orthogonality of the gradients of the basis functions proves the analogous claim for coefficients associated to ψ R 2i , i ∈ N, see (3.2.15) and L 2 -orthogonality of (3.1.10) combined with (3.1.11). For basis functions with odd indices one can use the orthogonality given in (3.2.14). Assume (H2). We know the following convergence behaviour of the projectors defined in (5.9), [16], Moreover, for orthogonal projections one has stability estimates in a standard way. Therefore independently of n ∈ N, For the last stability estimate one has to use Minkowski's inequality first, splitting the various components of the basis functions into the three groups that belong either to S, S o or V, respectively. One additionally has the stability estimate where C is independent of n ∈ N. Moreover, one has (5.14) 10 Page 16 of 54 G. Grün, P. Weiß JMFM We now introduce our Galerkin scheme for approximate solutions. For this, we make the ansatz for c n , R n , m n similarly. As an example, we write down the Galerkin scheme of the Navier-Stokes equations in detail. We look for α u := (α u i , . . . , α u 2n ) such that for all j = 1, . . . , 2n, the equations Therefore, (5.15) can be written as a system of n ordinary differential equations in explicit form. We prefer a short notation for the full system using u n , c n , R n , m n . This way, we are looking for (u n , c n , R n , m n ) ∈ X n such that for all j = 1, . . . , 2n.
Remark 5.5. In contrast to (1.1a) an additional term −D´Ω(c n (t)) s ∇Π Cn (g L s ) (c n (t)) · Ψ u j dx appears on the right-hand side of (5.16a). This term is required to obtain stability estimates. In the limit n → ∞, it becomes a gradient. Hence it vanishes as the test functions are assumed to be solenoidal.
Note that the stiffness matrix in (5.16c) is invertible as well due to Poincaré's inequality and (5.16c) is not part of the ordinary differential equation. Instead its solution is just a function of the other unknowns and therefore weak differentiability of h a and, consequently, h n := ∇R n is sufficient. Moreover, the susceptibility (2.1), and the operator (·) s (5.4) are Lipschitz-continuous. All other nonlinear terms are obviously locally Lipschitz-continuous.
As our space dimension is at most d ≤ 3, by Sobolev's embedding all terms are well-defined. In detail, the regularity of our basis functions implies that the unknowns ∇u n , c n , ∇m n , ∇∇R n | Ω are L ∞functions. Therefore, integrability of the terms in (5.16) is evident. Naturally, due to the density results Lemma 5.6. For any initial data u 0 n ∈ U n , c 0 n ∈ C n , m 0 n ∈ M n , system (5.16) has a global solution. Proof. By the Picard-Lindelöf Theorem, the system above has a unique local solution that attains any prescribed initial data u 0 n ∈ U n , c 0 n ∈ C n , m 0 n ∈ M n . The local solution is indeed a global solution on [0, T ] due to a priori estimates which we will prove next.
For fixed t ∈ [0, T ], we multiply (5.16) by the coefficient functions α (·) j (t) at time t and sum up over all j = 1, . . . , 2n. Note that we have (∇R n )| Ω ⊂ M n , which makes it possible to take the magnetic field h = ∇R| Ω , R ∈ R n as testfunction for the magnetization m ∈ M n .
As a first step we test • and the weak time derivative-note that the coefficients α (·) (·) are weakly differentiable and h a is weakly differentiable in time-of (5.16c) with μ 0 α 1 R n . For the ease of presentation we use the abbreviations A rather involved computation is related to the nonlinear testfunction g s (c n ) which has to be projected onto C n . By H 1 -regularity of c n and linear growth of (g L s ) (c n ), the term Π Cn (g L s ) (c n ) is well defined. One gets, The other computations are straightforward and we easily get and and Step 2: We verify the identity −ˆΩ div m n div ∇R n dx = div ∇R n 2 L 2 (Ω \∂Ω) .
First, we need more information about the term (div ∇R n − (div ∇R n ) Ω ) ∈ R, where (div ∇R n ) Ω is the mean value of div ∇R n on Ω . We quickly check that div ∇R n | Ω \Ω ≡ 0 and div ∇R n | Ω ∈ H 1 0 (Ω) which is a consequence of the choice of the basis functions for R, (3.2.12) 2 . Their gradients on Ω are associated with the basis of H, see (3.2.1), and the divergence of those functions vanishes at the boundary ∂Ω. Hence, we can easily deduce weak differentiability on the whole domain Ω . In preparation, we consider the part of the basis (3.2.12) whose gradients generate the space S.
(3.2.15)-just as div ∇R n − (div ∇R n ) Ω does. With the aforementioned basis functions ψ R 4i , we represent div ∇R n − (div ∇R n ) Ω with convergence of the infinite sum in the norm of R. We now know that there is a sequence (ψ k ) k∈N , such that Due to the boundary conditions of (3.2.12), the L 2 -orthogonality of (3.1.13) and the considerations above, we get In order to proceed, we need more information about ψ n . We easily obtain ψ n | Ω \Ω ≡ const. =: C as the basis functions ψ R 4i are constant on Ω \Ω. Therefore, and Combining the computations we get (5.22).
Step 3: Estimates on terms containing h a . Next, we will integrate in time and apply Young's inequality on terms of the type h a · ∇R n , ∂ t h a · ∇R n , m n · Π Hn h a .
In detail, Note that integration in time is possible due to continuity of the solutions t → α (·) (·) (t). Actually, we can only integrate in time until some t * ∈ (0, T ), but we can later deduce that we were able to integrate until T as the solutions will be bounded in time and therefore still exist for even larger times. We arrive at In the following, we first describe how to proceed, then give the intermediate steps.
Recall, the solution of (5.16c) is just a function of the other unknowns, i.e. R n (0) ∈ R n is defined as mean value free solution ofˆΩ Therefore, for any ε > 0, Due to the convergence m 0 n → m 0 in L 2 (Ω) d one can easily bound m 0 n in L 2 (Ω) d . Analogously, u 0 n is bounded in L 2 (Ω) d . The regularised entropy can be bounded by a quadratic function, therefore analogously g L s (c 0 n ) is bounded in L 2 (Ω). The projection Π Hn h a is bounded in L 2 (I × Ω) d , as a consequence of the stability (5.13) of Π Hn and regularity of the external magnetic field h a , (H2).
It remains to deal with the first four terms of the right-hand side. Therefore, we compute The h a -term is bounded due to convergence (5.14) and h a | Ω ∈ L 2 (I; H 3 (Ω)), cf. (H2). Hence, absorption is possible. The termˆT 0ˆΩ χ(c n , ∇R n )∇R n · Π Hn h a dx dt can be estimated analogously where we exploit the boundedness of χ, (H4). For the term ∂ t m n · Π Hn h a we rewriteˆΩ On the right-hand side the newly added term is bounded due to convergence Π Hn h a (T ) → h a (T )| Ω in L 2 (Ω) d . On the right-hand side the term − μ0β 2τ rel´Ω m 0 n · Π Hn h a (0) dx appears but can be estimated by Young's inequality and the same arguments as before. Note that the terms h a (0), h a (T ) are well-defined, as the space of the time variable is one-dimensional and h a is weakly differentiable. The second term on the right-hand side of (5.24) can be dealt with by absorption and boundedness of ∂ t Π Hn h a L 2 (I×Ω) d . Indeed, for any function ψ ∈ C ∞ 0 ((0, T )) and the basis representations Hence, ∂ t α 2i = β 2i and therefore ∂ t Π Hn h a = Π Hn ∂ t h a . The boundedness of that term follows analogously as before. Young's inequality applied to the term of the type χ(c n , ∇R n )∇R n · m n combined with boundedness (H4) of χ makes it possible to achieve an estimate using Gronwall's inequality later on. We end up with JMFM Transport of Magnetic Nanoparticles in Flow Page 23 of 54 10 Note that from the considerations before it follows that the h a -terms and the initial data are bounded. By Gronwall's inequality, μ 0 α 3 4 m n (T ) 2 L 2 (Ω) d ≤ C ha,initial (1 + e CT ). Note that the final time T is arbitrary and in a standard way L ∞ -in-time-estimates for u n , c n , m n can be achieved. We note, that easily one estimates all pure h a -terms (without projector) by the H 1 (I; Therefore, the solution exists on the whole time interval [0, T ].

Compactness Results
In this section, we establish compactness in time and in space necessary for the limit procedures in the discrete model and in the TMR-model. The starting point is the estimate and the constant C > 0 does not depend on σ c , L, s but only on h a , initial data and T .

Compactness in time
In this paragraph we establish estimates for ∂ t u n , ∂ t c n and ∂ t m n .
The claim follows asM → H 3 (Ω) d → W 1,∞ (Ω) d and Proof. We use the stability of the projection operator Π Un and standard estimates to get the result. Let Ψ ∈ L 2 (I; U). Also note that T 0ˆΩ ∂ t u n · Ψ dx dt =ˆT 0ˆΩ ∂ t u n · Π Un Ψ dx dt =: J Then, we compute . Due to H 3 -stability of Π Un , cf. (5.11), the claim follows.
Next, we are concerned with time-compactness of the magnetic field. Denote Proof. Equation (5.16c) can be differentiated with respect to time, so we get Let ∇ψ ∈ L 2 (I; H special ), then ∇ψ(t)| Ω ∈ S for almost all t ∈ I and therefore from the three types of basis functions within (3.2.12) one only needs those with even index, and from those only every second element-those whose gradients generate S. This implies that ∇Π Rn ψ| Ω \Ω ≡ 0, cf. (3.2.15).
Recall the L 2 -orthogonality of (3.1.13), which has been used above and also recall the relation between the bases of M, H, S and R from Section 3. In particular note that on Ω the potentials of the basis functions of H are used to define the basis functions of R, see (3.2.1). We use (5.16c) and estimate = C( Π Hn ∇ψ L 2 (I×Ω) d + Π Hn ∇ψ L 2 (I;H 3 (Ω) d ) ).
The first term on the right-hand side is bounded by the other one and the proof is finished by using stability of the projections due to H 3 -orthogonality of the corresponding basis functions of S.
Let us specify an assumption on f p (c n ) such that the conditions in Lemmas 6.1, 6.3 and 6.4 are satisfied. For this, we first state (without proof) an immediate consequence of our assumption of g L s (c). we infer the condition m ∈ [0, 2]. Note that if we do not use the L ∞ (I; L 2 (Ω))-bound of c but the L ∞ (I; L 1 (Ω))-bound which is independent of L, then we arrive at the condition m ∈ [1, 2] and a slightly less regular time derivative ∂ t c n .
As the gradients of m n and h n are only locally bounded in Ω, cf. [9], for an application of the Aubin-Lions lemma we will consider the functions m n φ and h n |Ω, where φ ∈ C ∞ 0 (Ω) is a cut-off function and Ω ⊂⊂ Ω. We have global estimates for the time derivatives of m n and h n . Therefore, we expect to obtain estimates of ∂ t (m n φ) and ∂ t (h n |Ω) as well. This is done in the subsequent corollaries.
Therefore, we try to recycle as much as possible from the original computations in Lemma 6.3 where we put φ besides the testfunction. LetΨ ∈M. Then we obviously have (Ψφ) ∈M. Therefore, we can just plug in the new testfunction (Ψφ), where Ψ ∈ L 2 (I;M) and proceed like in the proof of Lemma 6.3. Corollary 6.9. Under the assumptions of Lemma 6.5, for any V ⊂ Ω there is a constant C > 0 such that and the result follows.

Existence for the TMR-Model
We proceed taking the limit and identifying the terms of (5.16) with the terms of (5.7). Our strategy is to use sufficiently regular testfunctions, which will be projected onto the finite dimensional Galerkin approximation spaces as testfunctions in (5.16). The kind of testfunctions we are going to consider are Here, we used a notation like C ∞ 0 (Ω) d ∩ M to emphasize that such a space is considered under the norm of M.
We label all the terms of (5.16a) from left to right with L u 1 , . . . , L u 4 for the terms of the left-hand side and R u 1 , . . . , R u 3 for the terms of the right-hand side. The terms of the other equations will be labeled analogously, i.e. they will be labeled by L c does not exist in (5.7a) and is supposed to vanish. All other terms from (5.16) are in the same order as in (5.7). Moreover, we assume f 2 to be continuous with growth 0 < a 0 ≤ f 2 (c) ≤ a 1 |c| m + a 2 , m ∈ [0, 2] for some a 0 , a 1 , a 2 > 0 and |χ(·, ·)| ≤ χ max < ∞. We abbreviate h n := ∇R n and h := ∇R for the limit R of R n . From (6.1) we get the following convergence results for subsequences, which will not be relabeled for the ease of notation, and where s, L, σ c are fixed, in a standard way, see explanations below.
For the weak- * convergence in L ∞ (I; L 2 (Ω) d ) of c n , one can use Lemma 6.6. For the determination of the exponent q d , one uses Hölder's inequality as in the following Lemma 7.1 in the case d = 3 or a slightly better result from [11] in the case d = 2 in combination with the Aubin-Lions lemma [26] and Vitali's convergence theorem.
For application of the Aubin-Lions lemma estimates of the time derivatives need to be used, see Lemma 6.1, Lemma 6.3 and Corollary 6.8, Lemma 6.4, Lemma 6.5 and Corollary 6.9. We have a look at ∂ t c n , ∂ t u n first. From the time-compactness estimates we get converging subsequences in some L a (I; L b (loc) (Ω) (d) )-spaces, e.g. a = 2 in case of ∂ t u n or a = 5 4 in the case of ∂ t c n , b = 2, that converge pointwise almost everywhere. By uniform estimates and Vitali's convergence theorem the convergence in L q d − (I; L q d − (Ω)) or L q d − (I; L q d − (Ω)) d , respectively, can be achieved. For the convergence of m n we have to use local estimates. The local weak H 1 -convergence of the magnetic variables comes from the well-known result see [14,Lemma 2.5]. Using the formulas div(φm) = ∇φ · m + φ div m, we find (φm) ∈ L 2 (I; H 1 0 (Ω)) for any scalar φ ∈ C ∞ 0 (Ω) and m ∈ H(div, curl)(Ω), according to (7.5). With Corollary 6.8 we obtain a strongly converging subsequence of (m n φ) that converges pointwise almost everywhere. By uniform (local) estimates we get the convergence of m n in L q d − (I; L q d − loc (Ω) d ). For ∂ t h n we use the fact that on any V ⊂⊂ Ω we have ∂ t h n | V and h n | V bounded in L 2 (I; (∇[H 4 0 (V )]) ) or L 2 (I; H 1 (V ) d ), respectively, and obtain by the same methods a strongly converging subsequence in L q d − (I; L q d − (V ) d ) and therefore the local strong convergence.
With those convergence results-(7.3) and (5.14)-at hand, we can easily identify the limits in the linear terms L u 2 , L c 3 , L m 3 , L m 4 and the nonlinear terms L u 3 , L u 4 . For the terms L R 1 , R R 1 , R R 2 we proceed as follows. Let S ∈ R, then from (5.16c) we infer Ω ∇R n (t) · ∇Π Rn S dx =ˆΩ h a (t) · ∇Π Rn S dx −ˆΩ m n (t) · ∇Π Rn S dx. We multiply (7.7) with a function ϕ ∈ C ∞ 0 ((0, T )) and integrate in time. As ϕ does not depend on the spatial variable, we can writê We take the limit by exploiting weak convergence of ∇R n , m n in L 2 (I × Ω) d and strong convergence of the gradients of the projections in L 2 (I × Ω) d . Indeed, the projections ∇Π Rn (Sϕ) = ϕ∇Π Rn S converge pointwise in time and ϕ ∇Π Rn S − ∇S L 2 (Ω ) d is dominated (in time) by a constant, which is integrable on [0, T ]. By Lebesgue's dominated convergence theorem the projections converge in L 2 (I × Ω) d and we end up withˆT By the fundamental lemma of calculus of variations, for almost all t ∈ I the equationŝ From (6.1) we also get a subsequence (which will not be relabeled for the ease of notation) such that As c n converges strongly in L q− (I × Ω), there is a subsequence (not relabeled) that converges pointwise almost everywhere, and so does because (·) s and f 2 are continuous. The same argument holds for the other terms occurring below. By uniform bounds-based on (7.2) and (6.1)-and Vitali's convergence theorem one can deduce the strong convergences are bounded uniformly in L 2 (I ×Ω) d . We are going to identify W now. Testing with smooth and compactly supported testfunctions We consider the first term of the right-hand side. ˆT The term J 1 tends to zero because of the strong convergence of f 2 (c n ) and the L 2 -boundedness of ∇Π Cn (g L s ) (c n ). The latter follows easily from the H 1 -stability of the projector Π Cn and the computation ∇(g L s ) (c n ) = (g L s ) (c n )∇c n , where (g L s ) is a bounded function and ∇c n is bounded in For the term J 2 we will prove the weak convergence ∇Π Cn (g L s ) (c n ) ∇(g L s ) (c) (for a subsequence without relabeling) in L 2 (I × Ω) d . From the boundedness of the term ∇Π Cn (g L s ) (c n ) we also get ∇Π Cn (g L s ) (c n ) w in L 2 (I × Ω) d . We identify the limit by testing with where we used the convergence of (g L s ) (c n ) in L 2 (I × Ω) and combined with the L 2 -stability of Π Cn , cf. (5.11). Hence, J 2 → 0. For the second term on the righthand side of (7.8) we make use of the compact support of Φ and the thereby applicable higher regularity/convergence results for the magnetic variables. All three magnetic contributions are of the same kind, so we will only look at (∇h n ) T m n as an example. The factor √ f2(cn) (cn)s converges strongly in any L r (I × Ω), r ∈ [1, ∞), the magnetic field part ∇h n converges weakly in L 2 loc (I × Ω) d×d and the magnetization part m n converges strongly in L q d − (I; L q d − loc (Ω) d ), note that q d > 2. Then, the identification of the limit is straight forward. Hence, Exploiting the local regularity of the magnetic variables we prove convergence of the convective term in (5.16d) It suffices to only consider the less regular part ˆT The first term of the right-hand side converges to zero due to weak convergence of (V part ) n in L 2 (I × Ω) d and the sufficient integrability of the smooth and compactly supported testfunction θ and m ∈ L q d − (I; L q d − loc (Ω)). In the second term higher regularity of m n will be used, hence by boundedness of ((V part ) n ) n∈N in L 2 (I × Ω) d and strong (local) convergence of m n → m in L q d − (I; L q d − loc (Ω) d ) this term converges to zero. The local strong convergence of m n was applicable due to the compact support of θ. In the last term we use boundedness of ((V part ) n ) n∈N in L 2 (I × Ω) d and (m n ) n∈N in L ∞ (I; L 2 (Ω) d ). For convergence it is sufficient that ∇Π Mn θ → ∇θ in L 2 (I; L ∞ (Ω) d ). But as Π Mn θ(t) → θ(t) in H 3 (Ω) d → W 1,∞ (Ω) d , this is true, see (5.14). The terms R m 1 and R u 2 can be dealt with analogously (note that in the setting of the Navier-Stokes equations Π Un v(t) → v(t) in H 3 (Ω) d as well) and for R m 2 only the partˆT 0ˆΩ χ(c n , h n )h n · Π Mn θ dx dt needs to be considered. For this, we extract a pointwise almost everywhere in I ×Ω converging subsequence of h n . The susceptibility is a continuous function (it has a continuous extension when the second argument is zero), hence for a pointwise almost everywhere in I × Ω converging subsequence (not relabeled) of c n and h n one easily gets pointwise convergence almost everywhere in I × Ω for χ(c n , h n ). By assumption (7.2) the susceptibility χ is bounded, hence χ(c n , h n ) converges in any L r (I × Ω), r ∈ [1, ∞). Then the convergence of this term is an easy consequence. The sum R u 1 + R u 3 is linked to the convective velocity (V part ) n . We have, see e.g. along the lines of the proof of Lemma 6.4, Note that we also used´Ω(∇m n ) T m n · Π Un v dx = 0 for this, which is true as we have div Π Un v = 0. The first factor converges strongly in L q d − (I × Ω), the second factor-with brackets around-converges weakly in L 2 (I × Ω) d and the projection of the testfunction converges strongly in L 5 (I × Ω) d . Note that 1 2 + 1 q d ≥ 1 5 . Hence, we have convergence towards the term where we used div v = 0 and (∇m) T m = 1 2 ∇|m| 2 and Note that due to the compact support of v in Ω the terms (∇m) T m, (∇h) T m and c s ∇(g L s ) (c) are sufficiently regular to be separated from each other. Rewriting term L c 4 with the definition of (V part ) n , (6.2), one gets a term similar to L c 2 but with less regularity. Hence it suffices to consider The first factor converges strongly in L q d − (I × Ω), the second converges weakly in L 2 (I × Ω) d , so the last factor has to converge strongly in L 5 (I × Ω) d . Concerning the convergence with respect to the spatial variable we use the embedding H 2 (Ω) d → L 6 (Ω) d (for d ≤ 3) and L 5 (I; H 2 (Ω))-convergence of Π Cn ψ, see (5.14). Hence, the term L c 4 converges. The convergence of the terms L u 1 , L c 1 and L m 1 are trivial consequences of Lemma 6.4, Lemma 6.1 and Lemma 6.3. Also, the convergence of solutions at time t = 0 towards initial data in the weak sense is obvious.
As the limit functions obey an energy estimate corresponding to weak lower semi-continuity of norms and sufficient weak convergence of all terms on the left-hand side of (6.1), one could prove timecompactness again without the necessity of projectors. Therefore, the stability results related to the H 3 -norm are not needed and a better time-regularity could be achieved. However, to keep it simple, we do not pursue such an approach. However, for easier accessibility, instead of ∂ t h ∈ L 2 (I; (H special ) ), see as well as χ ∞ < ∞ (see (7.2)) and d ∈ {2, 3} there exists a weak solution (u, c, R, m) as specified in Definition 5.1. Moreover, the weak solution satisfies the energy estimate + m L 2 (I;H(div,curl)(Ω)) + h L 2 (I;H(div,curl)(Ω \∂Ω)) Proof. This follows from the considerations made so far. The convergence of the terms in the Galerkin scheme to the terms of the weak formulation has been discussed before this theorem. The energy estimate follows from the weak lower semi-continuity of norms and the convergences of initial data and projections of various h a -terms. For further details, recall the right-hand side (5.25), from which the lim inf needs to be considered. Also, for non-negative initial data c 0 one easily obtains the estimate´Ω g L s (c 0 ) dx ≤ K c 0 2 L 2 (Ω) for some K > 0 and for terms of the kind h a (s) 2 L 2 (Ω) , s ∈ {0, T }, one can use Sobolev's embedding (with respect to time variable) in order to estimate those terms by the H 1 (I; L 2 (Ω ))-norm. The L 2 (I; H(div)(Ω))-norm of h a can be estimated by the H 1 (I; L 2 (Ω ))-norm as well due to div h a = 0. Hence, we obtain the estimate (7.11) of the constant on the right-hand side.

The Non-regularized Case
In this section, we study the limit problem (s, L −1 , σ c ) = (0, 0, 0). This requires a different approach to obtain regularity of particle densities c n . To fix further notation, we choose sequences σ c = σ c (n) := 1 n → 0, s = s(n) := 1 n → 0 , L = L(n) := 3n → ∞. We write {u n , c n , R n , m n } for the solutions of the regularized system that exist according to the Theorem 7.2. In this section we confine ourselves to the special case d = 2 and f 2 (c) ∼ c 2 , in detail we approximatively choose for any n ∈ N, where g s := g ∞ s and the bounds of c n follow from Lemma 6.6. Moreover, the particle velocity field is redefined as which is a consequence of (7.10), we deduced the . Terms like 1 3n 2 ∇c n L 2 (I×Ω) on the left-hand side have been omitted as they are not controlled uniformly. The remaining generic constant C on the right-hand side of (8.3) does not depend on s, σ c nor L. It only depends on h a and other given data.
As H(div, curl)(Ω) is not compactly embedded into L 2 (Ω) d (cf. [7]) and both m and h enter the system in a nonlinear way, we work with localized spaces-this way guaranteeing applicability of Aubin-Lion-type arguments to deduce strong convergence results. Analytically, we formulate our results both for d = 2 and d = 3 space dimensions. For a passage to the limit in three space-dimensions, however, apparently an additional regularization is needed. We discuss some options in Remark 8.11. Let us now improve regularity. . Proof. By boundedness of (V part ) n in L 2 (I × Ω) d and the identity (8.4) we infer that ∇c n has at least the regularity the three terms (∇m n ) T m n , (∇h n ) T m n , (∇h a ) T m n come along with (or regularity of (V part ) n if the latter is worse). Comparing regularity of m n , h n , h a it suffices to consider the term (∇m n ) T m n , only. According to (7.5), (7.6) we have m n ∈ L 2 (I; H 1 loc (Ω) d ). Together with m n ∈ L ∞ (I; L 2 (Ω) d ) we obtain 10 Page 36 of 54 G. Grün, P. Weiß JMFM for γ ∈ (0, 2) we infer-setting a := |∇m n |, b := |m n | and choosing γ = d+2 d+1 -that ∇c n ∈ L γ (I; L γ loc (Ω) d ). By Sobolev's embedding, Using c n ∈ L ∞ (I; L 1 (Ω)), see Lemma 6.6, the claim follows by Lemma 7.1. Note that we used the estimate in Lemma 6.6 which is independent of L = L(n).
In order to identify the limit of the fluxes (c s ) n (V part ) n via strong convergence of c n and weak convergence of (V part ) n in L 2 (I × Ω) d higher regularity for c n is needed.
c > e. (8.5) Proof. We split the integration into two parts.
Note that time-compactness of the Galerkin solutions m n established before in the case d = 3 carries over to this setting. The reason is the fact that σ c (n), L(n), s(n) occur in the magnetization equation (5.7d) only as part of (V part ) n . But (V part ) n is bounded in L 2 (I × Ω) d , see (8.3), as it was in the Galerkin setting, too. Also, it is possible to carry over the estimates for ∂ t h n . The only reason for the restriction to the case d = 2 is the regularity of the particle density c n . Sufficient time-compactness of u n can be proven in dimension d = 3 as well. Note that we cannot use the same time-compactness estimates as in the Galerkin-setting for the velocity field right away, because the weak formulation and the Galerkin scheme are not analogous to each other. We are going to estimate globally in Ω, and therefore do not gain any improvements in the case d = 2.
Lemma 8.6. Let d = 3 and let (u n , m n , h n ) n∈N be bounded in and assume h a ∈ L ∞ (I; L 2 (Ω) d ). Then there is a constant C > 0 such that for all n ∈ N ∂ t u n L 2 (I;U ) ≤ C.
Proof. Take v ∈ L 2 (I; U), then ˆT So, we can conclude ˆT If we combine the Lemmas and Corollaries 8.1-8.6 of this paragraph, we get the following statement. In the case d = 2 there exists a subsequence n k → ∞ such that for some limit function c ∈ L 2 (I; L 2 loc (Ω)) c n k → c in L 2 (I; L 2 loc (Ω)), c n k → c pointwise almost everywhere in Ω.
Note that the time-compactness estimate of ∂ t h n | V , for V ⊂⊂ Ω, carries over as well. The convergence behavior of the Galerkin solutions carries over to our sequences u n , m n and h n due to analogous energy estimates. For the functions c n new results have been obtained in this paragraph. Therefore as a starting point we use the convergences Based on the stability estimate (8.3), we obtain a non-negativity result for the limit c of the particle density functions. Therefore c ≥ 0.
Additionally, we have the following straight forward pointwise convergence results. The regularised functions (c n ) s satisfy the same estimates from Lemma 8.2, as the function G that was used in the lemma cannot distinguish between c n and (c n ) s(n) . Moreover the L ∞ (I; L 1 (Ω))-estimate and the L 4 3 (I; L 4 loc (Ω))-estimate of c n can be carried over to (c n ) s trivially. Hence, (c n ) s(n) → c in L 2 (I; L 2 loc (Ω)), (8.8) too.
We will now pass to the limit. We will plug in C ∞ 0 ([0, T ); C ∞ 0 (Ω)) 2 -testfunctions v for the Navier-Stokes equations and easily identify the left-hand side of the equations. Exploiting the compact spatial support of the testfunction, we have no problems with the right-hand side, either. As an example let us consider the termˆT 0ˆΩ (m n × h n ) · curl v dx dt.
As we can use h n → h in L q (I; L q loc (Ω) d ) and m n → m in L q (I; L q loc (Ω) d ), where q > 2 in both cases d = 2 or d = 3, there is no obstacle in taking the limit. For (5.7b) we analogously take smooth and compactly supported testfunctions.
Many considerations are identical to the case when the Galerkin solutions converged to the weak solution of the regularized system. Note that the convergence behavior in two dimensions is at least as good as in three dimensions and therefore we only need to concentrate on the terms with c n in it. First, we will restore the weak convergence result of (V part ) n in L 2 (I × Ω) d . We easily get (V part ) n W JMFM Transport of Magnetic Nanoparticles in Flow Page 41 of 54 10 in L 2 (I × Ω) d for some W ∈ L 2 (I × Ω) d , see (8.3). We identify W in the same way as done in the Galerkin-setting. Let Φ ∈ C ∞ 0 (I × Ω) d , then

→1
(α 1 ∇h n + β 2 ∇h a − α 2 ∇m n ) T m n · Φ dx dt. (8.9) First we note that the two terms that had been combined into one term before can be separated into individual terms due to the regularity of ∇c n . The second term of the right-hand side converges as the gradients converge weakly and locally with respect to L 2 (I × Ω) d and m n converges strongly and locally with respect to L 2+ (I × Ω) d (in any case d ∈ {2, 3}) and the bounded quotient 0 < (cn) s(n) ≤ 1 converges strongly in any L p (I × Ω)-norm, p ∈ [1, ∞), due to pointwise convergence (for a non-relabeled subsequence), see Lemma 8.9, and Vitali's convergence theorem. Those convergences are sufficient to identify the limit. The first term converges, obviously, and we arrive at W = −KD∇c + μ 0 K(α 1 ∇h + β 2 h a − α 2 ∇m) T m. With this at hand, we can proceed as in the Galerkin case for any term except for those terms of the particle density equation (5.7b). We can use our bound on ∇c n ∈ L 5/4 (I; L 5/4 (Ω) d ) (in both cases d ∈ {2, 3}) and the smoothness of the testfunction in order to prove that the third term of (5.7b)-the regularizing term-vanishes. Now, consider the first term. Integration by parts giveŝ T 0 ∂ t c n , ψ (H 2 * (Ω))×H 2 * (Ω) dt = −ˆT 0ˆΩ c n ∂ t ψ dx dt +ˆΩ c 0 ψ(0) dx with testfunctions ψ ∈ C 1 0 ([0, T ); C 2 0 (Ω)). This can be justified already on the level of the Galerkin approximation. Taking the limit is straightforward. The fourth term in (5.7b) simplifies tô T 0ˆΩ (c n ) s (V part ) n · ∇ψ dx dt, which converges as (c n ) s(n) → c in L 2 (I; L 2 loc (Ω)) and (V part ) n V part in L 2 (I × Ω) d . This step was the only one where we needed to restrict ourselves to the case d = 2. The second term is easier and analogous to the fourth term (no restriction to d = 2 needed). Hence, the limit has been taken and we obtain the following result.  Remark 8.11. It would be desirable to have existence of weak solutions in the three-dimensional setting, too. Recall that the bottleneck is the space-time integrability of c. Due to the intricate coupling between the evolution equations for c and m, this integrability is in a subtle way related to the regularity of (∇m) T m or (∇h) T m and hence as a consequence related to the regularity of m. Let us discuss two methods of regularization to overcome this issue.