Existence of weak solutions to a diffuse interface model involving magnetic fluids with unmatched densities

Abstract

In this article we prove the global existence of weak solutions for a diffuse interface model in a bounded domain (both in 2D and 3D) involving incompressible magnetic fluids with unmatched densities. The model couples the incompressible Navier–Stokes equations, gradient flow of the magnetization vector and the Cahn–Hilliard dynamics describing the partial mixing of two fluids. The density of the mixture depends on an order parameter and the modelling (specifically the density dependence) is inspired from Abels et al. (Models Methods Appl Sci 22(3):1150013, 2011).

1 Introduction

This article is devoted to the mathematical analysis of a system modeling the flow of two viscous incompressible fluids with magnetic properties and unmatched densities undergoing partial mixing in a bounded domain $\Omega \subset {\mathbb {R} }^d$, $d=2,3.$ with the boundary $\partial \Omega$ of class $C^{2}$. Let $T>0$ and define the space time cylinder as $Q_T = \Omega \times (0,T)$. Further let $\Sigma _{T}$ denote $\partial \Omega \times (0,T)$.

The mixing of the fluids is described by an order parameter $\phi :Q_T \rightarrow {\mathbb {R} }$, which is the difference of the volume fractions of the fluids involved and undergoes a smooth and rapid transition in an interfacial region between the two fluids. We denote by $v:Q_T\rightarrow {\mathbb {R} }^d$ the mean fluid velocity, by $\rho =\rho (\phi ):Q_{T}\rightarrow {\mathbb {R} }$ the mean mass density, $p: Q_T\rightarrow {\mathbb {R} }$ the pressure, $M:Q_{T}\rightarrow \mathbb {R}^{3}$ the magnetization and $\mu :Q_T \rightarrow {\mathbb {R} }$ the chemical potential. The system we consider reads as follows

$$\begin{aligned} \begin{array}{llll} &{}\displaystyle \partial _{t}(\rho v)+\text{ div }(\rho v\otimes v)-\text{ div }\,(2\nu (\phi )\mathbb {D}(v))+\text{ div }(v\otimes J)+\nabla p&{}\\ &{}\quad \displaystyle =\mu \nabla \phi +\frac{\xi (\phi )}{\alpha ^{2}}((|M|^{2}-1)M)\nabla M-{{\,\textrm{div}\,}}(\xi (\phi )\nabla M)\nabla M&{} \text{ in } Q_{T},\\ &{}\displaystyle {{\,\textrm{div}\,}}v=0&{} \text{ in } Q_{T},\\ &{} \displaystyle \partial _{t}M+(v\cdot \nabla )M={{\,\textrm{div}\,}}(\xi (\phi )\nabla M)-\frac{\xi (\phi )}{\alpha ^{2}}(|M|^{2}-1)M&{} \text{ in } Q_{T},\\ &{}\displaystyle \partial _{t}\phi +(v\cdot \nabla )\phi =\Delta \mu &{} \text{ in } Q_{T},\\ &{}\displaystyle \mu =-\eta \Delta \phi +\Psi '(\phi )+\xi '(\phi )\frac{|\nabla M|^{2}}{2}+\frac{\xi '(\phi )}{4\alpha ^{2}}(|M|^{2}-1)^{2}&{} \text{ in } Q_{T},\\ &{}\displaystyle v=0,\ \partial _{n}M=0,\ \partial _{n}\phi =\partial _{n}\mu =0&{} \text{ on } \Sigma _{T},\\ &{}\displaystyle (v,M,\phi )(\cdot ,0)=(v_{0},M_{0},\phi _{0})&{} \text{ in } \Omega . \end{array} \end{aligned}$$

(1.1)

where J is a relative flux related to the diffusion of the components and is given by the following expression:

$$\begin{aligned} \begin{array}{l} \displaystyle J=-\frac{{\widetilde{\rho }}_{2}-{\widetilde{\rho }}_{1}}{2}\nabla \mu , \end{array} \end{aligned}$$

(1.2)

with ${\widetilde{\rho }}_{i}$ ($i=1,2$) denoting the specific density of the $i-$th fluid. In system (1.1), $\alpha$ and $\eta$ are positive constants, $\nu (\phi )$ is the concentration dependent viscosity coefficient which is assumed to be non degenerate and $\mathbb {D}(v)=\frac{1}{2}\left( \nabla v+(\nabla v)^\top \right)$ denotes the symmetric part of the velocity gradient. The factor $\alpha$ penalizes the saturation condition of the length of the magnetization vector |M| from 1 and $\eta$ measures the thickness of the region where the two fluids mix. The function $\xi (\phi )$ denotes the exchange parameter, which reflects the tendency of the magnetization to align in one direction. We assume that $\xi (\cdot )$ is non degenerate, i.e., that it has a positive lower bound and further both $\xi$ and $\xi '$ are bounded from above, cf. (1.6). The homogeneous free energy density of the fluid mixture is denoted by $\Psi (\phi ).$ We will consider a class of singular free energies and our consideration (cf. (1.8)) will include the homogeneous free energy of the form

$$\begin{aligned} \begin{array}{l} \displaystyle \Psi (s)=\frac{a}{2}\left( (1+s)\ln (1+s)+(1-s)\ln (1-s)\right) -\frac{b}{2}s^{2}, \end{array} \end{aligned}$$

(1.3)

where $s\in [-1,1]$ and $a,b>0,$ introduced in [17]. The logarithmic terms in (1.3) relate to the entropy of the system. We note that the function $\Psi$ in (1.3) is convex iff $a\geqslant b$ and $\Psi '$ shows singular behavior at $s=\pm 1$.

The present article is devoted to prove the existence of a global weak solution (i.e. without any restriction on time and the size of the initial data) of the model (1.1). In Sect. 1.1 we will first introduce some notations corresponding to the functional spaces, the Leray projector which will be essential to deal with our incompressible bi-fluid model and then present our result on the global well posedness of the system (1.1). In Sect. 1.2 we will comment on the strategies of the proof and related approaches. After a discussion on the physical background of the system in Sect. 1.3, we devote Sect. 1.4 to bibliographical notes.

1.1 Functional framework and main result

1.1.1 Functional settings

The Lebesgue and Sobolev spaces are denoted by the notations $L^{p}(\Omega )$ and $W^{s,p}(\Omega )$ respectively. The standard notations $C^k({\overline{\Omega }})$ and $C^{k,\gamma }({\overline{\Omega }})$ ($\gamma \in (0,1]$) are used to denote respectively the spaces of k–times continuously differentiable functions and Hölder continuous functions. The functional spaces with elements having compact support are denoted by using a subscript c. The hooked arrow notations $X\hookrightarrow Y$ ($X{\mathop {\hookrightarrow }\limits ^{C}} Y$) are used to write the continuous (compact) embedding of a Banach space X to a Banach space Y. The duality pairing between a Banach space X and its dual $X'$ is written as $\left\langle \cdot ,\cdot \right\rangle .$ The functional space $C_w([0,T];X)$ denote a subspace of $L^\infty (0,T;X)$ containing functions f for which the mapping $t\mapsto \left\langle \phi , f(t)\right\rangle$ is continuous on [0, T] for each $\phi \in X'$. We set

$$\begin{aligned}&L^{2}_{{{\,\textrm{div}\,}}}(\Omega )=\overline{\{v\in C^\infty _c(\Omega ):\ {{\,\textrm{div}\,}}v=0\text { in }\Omega \}}^{\Vert \cdot \Vert _{L^2}},\\&W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )=\overline{\{v\in C^\infty _c(\Omega ):\ {{\,\textrm{div}\,}}v=0\text { in }\Omega \}}^{\Vert \cdot \Vert _{W^{1,2}}},\\&W^{2,2}_n(\Omega )=\{u\in W^{2,2}(\Omega ): \partial _n u=0\text { on }\partial \Omega \},\\&V(\Omega )=W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )\cap W^{2,2}(\Omega ). \end{aligned}$$

From now onward we will use C to denote a generic positive constant which might vary from line to line. For simplicity of notations we will always use $\cdot$ to denote both the scalar product of vectors and tensor products.

Let us now introduce the Leray projector $\mathbb {P}_{{{\,\textrm{div}\,}}}:L^{2}(\Omega )\rightarrow L^{2}_{{{\,\textrm{div}\,}}}(\Omega )$ defined as

$$\begin{aligned} \begin{array}{l} \displaystyle \mathbb {P}_{{{\,\textrm{div}\,}}}(f)=f-\nabla p\,\, \text{ for }\,\, f\in L^{2}(\Omega )\, \text{ where }\, p\in W^{1,2}(\Omega )\,\, \text{ with } \,\,\int _{\Omega } p=0 \end{array} \end{aligned}$$

(1.4)

solves the weak Neumann problem

$$\begin{aligned} \begin{array}{l} \displaystyle \left( \nabla p,\nabla \varphi \right) _{\Omega }=\left( f,\nabla \varphi \right) \,\,\text{ for } \text{ all }\,\,\varphi \in C^{\infty }(\overline{\Omega }). \end{array} \end{aligned}$$

(1.5)

1.1.2 Existence of weak solutions

In order to state the precise definition of a weak solution to (1.1) we summarize the assumptions on the functions $\xi ,$ $\nu ,$ $\Psi$ and also make the dependence of the mean density $\rho (\phi )$ on $\phi$ precise.

Assumption 1.1

The function $\xi \in C^1({\mathbb {R} })$ satisfies

$$\begin{aligned} \begin{aligned}&0<c_1\leqslant \xi \leqslant c_2\text { on }{\mathbb {R} }, \text { for some }c_1,c_2>0,\\&\xi '\leqslant c_3\text { on }{\mathbb {R} }, \text { for some }c_3>0.\ \end{aligned} \end{aligned}$$

(1.6)

The viscosity coefficient $\nu \in C^1({\mathbb {R} })$ satisfies

$$\begin{aligned} 0<\nu _1\leqslant \nu \leqslant \nu _2\text { on }{\mathbb {R} }, \text { for some }\nu _1,\nu _2>0. \end{aligned}$$

(1.7)

The homogeneous free energy density $\Psi \in C([-1,1])\cap C^2((-1,1))$ solves

$$\begin{aligned} \begin{aligned}&\lim _{s\rightarrow -1}\Psi '(s)=-\infty ,\ \lim \limits _{s\rightarrow 1}\Psi '(s)=\infty ,\\ {}&\Psi ''(s)\geqslant -\kappa , \text { for some }\kappa \in \mathbb {R}. \end{aligned} \end{aligned}$$

(1.8)

The mean mass density $\rho$ and the phase field $\phi$ are related via

$$\begin{aligned} \rho (\phi )=\frac{1}{2}({\widetilde{\rho }}_{1}+{\widetilde{\rho }}_{2})+\frac{1}{2}({\widetilde{\rho }}_{2}-{\widetilde{\rho }}_{1})\phi \,\,\text{ in }\,\, {\overline{Q}}_{T}, \end{aligned}$$

(1.9)

where ${\widetilde{\rho }}_{i}>0$, $i=1,2$ are specific constant mass densities of the unmixed fluids.

Remark 1.2

The expression for $\rho$ in (1.9) implies that

$$\begin{aligned} 0<\min \{{{\widetilde{\rho }}}_1,{{\widetilde{\rho }}}_2\}\leqslant \rho (\phi )\leqslant \max \{{{\widetilde{\rho }}}_1,{{\widetilde{\rho }}}_2\} \end{aligned}$$

(1.10)

provided that $|\phi |\leqslant 1$.

Remark 1.3

Let $\xi _{1},\xi _{2}>0$ be the exchange constants for the magnetic fluids undergoing partial mixing. The function

$$\begin{aligned} \xi (\phi )=(1-{\mathcal {H}}_{\eta }(\phi ))\xi _{1}+{\mathcal {H}}_{\eta }(\phi )\xi _{2}, \end{aligned}$$

where $\eta >0$ correspond to the thickness of the interface and ${\mathcal {H}}_{\eta }(x)=\frac{1}{1+e^{-\frac{x}{\eta }}}$ is a regularized approximation of the Heaviside step function, provides an example of a non degenerate function $\xi$ satisfying the assumptions (1.6), cf. [36, 42, 50], in which such a regularized Heaviside function is used in a similar context.

Next we define the notion of weak solution to the system (1.1).

Definition 1.4

(Definition of weak solutions) Let $T>0.$ For a given triplet

$$\begin{aligned} (v_0,M_0,\phi _0)\in L^2_{{{\,\textrm{div}\,}}}(\Omega )\times W^{1,2}(\Omega )\times W^{1,2}(\Omega )\text { with }|\phi _0|\leqslant 1, \end{aligned}$$

(1.11)

we call the quadruple $(v,M,\phi ,\mu )$ possessing the regularity

$$\begin{aligned} \begin{aligned} v&\in C_w([0,T];L^2_{{{\,\textrm{div}\,}}}(\Omega ))\cap L^2(0,T;W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )),\\ M&\in C_w([0,T];W^{1,2}(\Omega ))\cap C([0,T];L^{2}(\Omega ))\cap W^{1,2}(0,T;L^{\frac{3}{2}}(\Omega )),\\ \phi&\in \displaystyle C_w([0,T];W^{1,2}(\Omega ))\cap C([0,T];L^{2}(\Omega ))\cap L^2(0,T;W^{2,1}(\Omega ))\\&\text { with }\Psi '(\phi )\in L^1(Q_T),\\ \mu&\in L^2(0,T;W^{1,2}(\Omega )), \end{aligned} \end{aligned}$$

(1.12)

a weak solution to (1.1) if it satisfies

$$\begin{aligned} \begin{aligned} \int _{\Omega }&\rho (t)v(t)\cdot \psi _1(t)-\int _{\Omega }\rho _0v_0\cdot \psi _{1}(0)=\int _0^t\int _{\Omega }\biggl (\rho v\cdot \partial _t\psi _{1} +(\rho v\otimes v\\&+v\otimes J)\cdot \nabla \psi _{1}-2\nu (\phi )\mathbb {D}v\cdot \mathbb {D}\psi _{1}-\nabla \mu \phi \cdot \psi _{1} +\Bigl (\frac{\xi (\phi )}{\alpha ^{2}}\bigl ((|M|^{2}-1)M\bigr )\nabla M\\&-{{\,\textrm{div}\,}}\bigl (\xi (\phi )\nabla M\bigr )\nabla M\Bigr )\cdot \psi _{1}\biggr ),\\ \int _\Omega&M(t)\cdot \psi _2(t)-\int _\Omega M_0\cdot \psi _{2}(0)=\int _0^t\int _\Omega \biggl (M\cdot \partial _t\psi _{2}-(v\cdot \nabla ) M\cdot \psi _{2}\\&-\xi (\phi )\nabla M\cdot \nabla \psi _{2}-\frac{1}{\alpha ^2}\bigl (\xi (\phi )(|M|^2-1)M\bigr )\cdot \psi _{2}\biggr ),\\ \int _\Omega&\phi (t)\psi _3(t)-\int _\Omega \phi _0\psi _3(0)=\int _0^t\int _\Omega \left( \phi \partial _t\psi _{3}-(v\cdot \nabla )\phi \psi _{3}-\nabla \mu \cdot \nabla \psi _{3}\right) ,\\&\mu -\xi '(\phi )\frac{|\nabla M|^{2}}{2}-\frac{\xi '(\phi )}{4\alpha ^{2}}(|M|^{2}-1)^{2}=-\eta \Delta \phi +\Psi '(\phi )\text { a.e. in }Q_T \end{aligned} \end{aligned}$$

(1.13)

for all $t\in (0,T)$, for all $\psi _{1}\in C^{1}_{c}([0,T);V(\Omega ))$, $\psi _{2}\in C^{1}_{c}([0,T);W^{1,2}(\Omega ))$ and all $\psi _{3}\in C^{1}_{c}\left( [0,T);W^{1,2}(\Omega )\right)$. The initial data are attained in the form

$$\begin{aligned} \lim _{t\rightarrow 0_+}\left( \Vert v(t)-v_0\Vert _{L^2(\Omega )}+\Vert M(t)-M_0\Vert _{W^{1,2}(\Omega )}+\Vert \phi (t)-\phi _0\Vert _{W^{1,2}(\Omega )}\right) =0.\nonumber \\ \end{aligned}$$

(1.14)

Now we present the central result of our article that concerns the global existence of a weak solutions of the model (1.1).

Theorem 1.5

Let $T>0$, $\Omega \subset {\mathbb {R} }^d$ be a bounded domain of class $C^{2}$ and let the initial datum $(v_0,M_0,\phi _0)$ satisfy (1.11). Then under Assumption 1.1 there exists a quadruple $(v,M,\phi ,\mu )$ which solves (1.1) in the sense of Definition 1.4. Moreover there exists a $p>2$ such that the triplet $(M,\phi ,\Psi '(\phi ))$ enjoys the following improved regularity

$$\begin{aligned} \begin{aligned} M&\in L^2(0,T;W^{1,p}(\Omega )),\\ \phi&\in L^2(0,T;W^{2,\frac{2p}{p+2}}(\Omega )),\\ \Psi '(\phi )&\in L^2(0,T;L^\frac{2p}{p+2}(\Omega )). \end{aligned} \end{aligned}$$

(1.15)

Further the following items hold:

• (i) The obtained weak solution $(v,M,\phi ,\mu )$ of (1.1) satisfies the following energy inequality:

  $$\begin{aligned} \begin{aligned}&E_{tot}(v(t),M(t),\phi (t))+\int _0^t\Biggl (\Vert \sqrt{2\nu (\phi )}\mathbb {D} v\Vert ^2_{L^2(\Omega )}+\Vert \nabla \mu \Vert ^2_{L^2(\Omega )}\\&\quad +\left\| {{\,\textrm{div}\,}}(\xi (\phi )\nabla M)-\frac{\xi (\phi )}{\alpha ^2}M(|M|^2-1)\right\| ^2_{L^2(\Omega )}\Biggr )\leqslant E_{tot}(v_0,M_0,\phi _0) \end{aligned} \end{aligned}$$

  (1.16)

  for all $t\in (0,T),$ where

  $$\begin{aligned} \begin{aligned} E_{tot}(v,M,\phi )=&\frac{1}{2}\int _{\Omega }\rho (\phi )|v|^{2}+\frac{1}{2}\int _{\Omega }\xi (\phi )|\nabla M|^{2}+\frac{1}{4\alpha ^{2}}\int _{\Omega }\xi (\phi )(|M|^{2}-1)^{2}\\&+\frac{\eta }{2}\int _{\Omega }|\nabla \phi |^{2}+\int _{\Omega }{\Psi }(\phi ) \end{aligned} \end{aligned}$$

  (1.17)

  with $\rho (\phi )$ being defined as in (1.9).

• (ii) The magnetization M attains the homogeneous Neumann boundary condition in a weak sense, i.e. for a.e. $t\in (0,T)$ the following holds

  $$\begin{aligned} \begin{array}{l} \gamma _{n}(\xi (\phi )\nabla M)=0\quad \text{ in }\quad (W^{\frac{1}{2},2}(\partial \Omega ))', \end{array} \end{aligned}$$

  (1.18)

  where $\gamma _{n}$ is the normal trace operator defined on $\{g\in L^{2}(\Omega ):{{\,\textrm{div}\,}}g\in L^2(\Omega )\}$ such that the following Stokes formula holds for a.e. $t\in (0,T):$

  $$\begin{aligned} \begin{aligned} \int _{\Omega }{{\,\textrm{div}\,}}(\xi (\phi )\nabla M)\cdot \psi _{2}=&-\int _{\Omega }\xi (\phi )\nabla M\cdot \nabla \psi _{2}\\&+\langle \gamma _{n}(\xi (\phi )\nabla M),\gamma _{0}\psi _{2}\rangle _{(W^{\frac{1}{2},2}(\partial \Omega ))',W^{\frac{1}{2},2}(\partial \Omega )} \end{aligned} \end{aligned}$$

  (1.19)

  for all $\psi _{2}\in W^{1,2}(\Omega )$, where $\gamma _{0}: W^{1,2}(\Omega )\rightarrow W^{\frac{1}{2},2}(\partial \Omega )$ is the Dirichlet trace operator.

• (iii) Moreover, if $M_{0}\in W^{1,2}(\Omega )\cap L^{r}(\Omega ),$ $r>6,$ then $M\in L^{\infty }(0,T;L^{r}(\Omega ))$ and

  $$\begin{aligned} \begin{array}{l} \displaystyle \Vert M(t)\Vert _{L^{r}(\Omega )}\leqslant \Vert M_{0}\Vert _{L^{r}(\Omega )}e^{\delta t}\quad \text{ for } \text{ all }\quad t\in [0,T] \end{array} \end{aligned}$$

  (1.20)

  for some positive constant $\delta >0.$ Additionally, if $M_{0}\in L^{\infty }(\Omega ),$ then

  $$\begin{aligned} \begin{array}{l} \displaystyle \Vert M(t)\Vert _{L^{\infty }(\Omega )}\leqslant \Vert M_{0}\Vert _{L^{\infty }(\Omega )}e^{\delta t}\quad \text{ for } \text{ all }\quad t\in [0,T]. \end{array} \end{aligned}$$

  (1.21)

We stress on the fact that the improved regularity (1.15), more precisely the uniform bounds on the suitable approximates of $(M,\phi ,\Psi '(\phi ))$ in the spaces mentioned in (1.15) play a key role in the passage of limit and recovering the weak formulations (1.13) (more precisely (1.13)$_{4}$).

To the best of our knowledge, [36, 42, 50] are the only articles in the literature studying diffuse interface models for magnetic fluids. The article [42] develops a simplified model describing the behavior of two-phase ferrofluid flows using phase field-techniques and present an energy-stable numerical scheme for the same. The authors of [42] further analyse the stability and the convergence of the numerical scheme developed and as a by-product they prove the existence of weak solutions of their model. In the article [50] the authors propose a diffuse interface model and finite element approximation for two-phase magnetohydrodynamic (MHD) flows with different viscosities and electric conductivities. Their model involves the incompressible Navier-Stokes equations, the Maxwell equations of electromagnetism and the Cahn–Hilliard equations. Unlike [42] and [50], in one of our previous articles [36], we studied a diffuse interface magnetic fluid model where the magnetization vector M is modeled by a gradient flow dynamics.

As far as we know, our current article is the first mathematical study of a diffuse interface model for a magnetic fluid with unmatched densities. We consider the model (1.1) where the mean mass density of the mixture depends on the order parameter $\phi$ via the formula (1.9). We show that the mass density $\rho (\phi )$ is always strictly positive and bounded. In order to do so we prove that the order parameter $\phi$ solves the physically reasonable bound $\phi (x,t)\in [-1,1]$ for a.e. $(x,t)\in Q_{T}.$

To this end it is important to use a singular potential $\Psi (\cdot )$, cf. (1.8), as a homogeneous free energy density of the mixture. Often in the literature this singular free energy is approximated by a suitable smooth free energy density. For instance, in [36], we considered a Ginzburg–Landau double-well potential $\tfrac{1}{4\eta }(|\phi |^{2}-1)^{2}$ instead of the singular potential $\Psi (\cdot ).$ But using such a polynomial potential one can not ensure that the order parameter $\phi$ stays in the physical reasonable interval $[-1,1]$ due to the lack of a comparison principle for fourth order diffusion equation and hence in particular can not deal with the model (1.1) comprising of fluids with unmatched densities.

Unlike the model considered in [50], in the present case (and also in the one considered in [36]) the magnetization M enters into the Cahn–Hilliard dynamics. Due to the presence of $|\nabla M|^{2}$ in the Cahn–Hilliard part, cf. (1.1)$_{5}$ we can not obtain $L^{2}(0,T;W^{2,2}(\Omega ))$ regularity of the order parameter $\phi$ and we only recover $\phi \in L^{2}(0,T;W^{2,q}(\Omega ))$ for some $q>1$ by a bootstrap argument using $L^{2}(0,T;W^{1,p}(\Omega ))$ ($p>2$) elliptic regularity for M.

1.2 Ideas of proof and related discussion

The proof of Theorem 1.5 is given in Sect. 3. It relies on an unconditionally stable time discretization scheme designed in Sect. 2.2. Before going into the analysis of a time discrete problem we first write the singular potential $\Psi$ as a perturbation of a convex function. This helps in reformulating the Cahn–Hilliard equation (1.1)$_{5}$ as the subdifferential of a convex potential and to use the monotone operator theory and regularity results for Cahn–Hilliard equation developed in [1, 11]. The reformulation is done in Sect. 2.1. Roughly speaking the content of Sect. 2.2 is to consider a sequence $0=t_0<t_1<\ldots<t_k<t_{k+1}<\ldots$, $k\in \mathbb {N}_0$ of equidistant nodes and next to construct a solution $(v_{k+1},M_{k+1},\phi _{k+1},\mu _{k+1})$ to a stationary problem (cf. (2.12)) at the point $t_{k+1}$ using $(v_k,M_k,\phi _k)$ which corresponds to the solution at the time $t_k$.

There is no common rule to write a time discretization of a nonlinear PDE model. It is generally done in a way so that the discrete system admits an energy inequality which is in close proximity with the formal energy balance of the original unsteady model. Here we follow a strategy devised in our previous article [36] to suitably discretize the term $\displaystyle \frac{\xi (\phi )}{\alpha ^{2}}(|M|^{2}-1)M$ appearing in (1.1)$_{3}.$ In order to obtain an energy type inequality for (2.12) one in particular tests $\displaystyle \frac{M_{k+1}-M_{k}}{h}$ (which is the discretization of the time derivative $\displaystyle \partial _{t}{M}$, with $h=t_{k+1}-t_{k}$ in the discrete magnetization equation (2.12)$_{3}$) with an approximation of $(|M|^{2}-1)M.$ Since the map $\displaystyle M\mapsto (|M|^{2}-1)M$ is not monotone one can check that

$$\begin{aligned} (M_{k+1}-M_{k})\left( |M_{k+1}|^{2}M_{k+1}-M_{k+1}\right) \ngeqslant \frac{1}{4}(|M_{k+1}|^{2}-1)^{2}-\frac{1}{4}(|M_{k}|^{2}-1)^{2} \end{aligned}$$

and hence the discretization $\displaystyle \frac{\xi (\phi )}{\alpha ^{2}}(|M|^{2}-1)M\approx \frac{{\xi (\phi _{k})}}{\alpha ^{2}}(|M_{k+1}|^{2}M_{k+1}-M_{k+1})$ does not lead to an unconditionally stable scheme. Following the convex splitting scheme used in our previous article [36] for vector valued functions, which is itself inspired from the convex splitting used in [50] for scalar functions, we use the approximation $\displaystyle \frac{\xi (\phi )}{\alpha ^{2}}(|M|^{2}-1)M\approx \frac{{\xi (\phi _{k})}}{\alpha ^{2}}(|M_{k+1}|^{2}M_{k+1}-M_{k})$, which along with Lemma 2.5 provides the desired estimate

$$\begin{aligned} (M_{k+1}-M_{k})\left( |M_{k+1}|^{2}M_{k+1}-M_{k}\right) \geqslant \frac{1}{4}(|M_{k+1}|^{2}-1)^{2}-\frac{1}{4}(|M_{k}|^{2}-1)^{2}. \end{aligned}$$

We then deal with the time discrete system (2.12) by considering it as a perturbation of a certain nonlinear operator and solving the operator equation by employing a fixed point argument.

After the proof of an existence result for the discrete problem (2.12) in Sect. 2.2, we define piecewise constant interpolants in Sect. 3 which approximate $(v,M,\phi ,\mu ),$ the solution to (1.13). The weak compactness of the interpolants are obtained from an energy type inequality and the strong compactness properties result by using the classical Aubin-Lions lemma and some suitable interpolation estimates. At this point a crucial observation is the obtainment of the strong convergence of $\{\nabla M^{N}\}_{N}$, where $M^{N}$ approximates M. This convergence plays a key role to pass to the limit in the terms approximating $\displaystyle \text{ div }(\xi (\phi )\nabla M)\nabla M$ (cf. (1.1)$_{1}$) and $\displaystyle \xi '(\phi )\frac{|\nabla M|^{2}}{2}$ (cf. (1.1)$_{5}$), as was presented in our previous article [36]. These arguments can be directly adapted to the current setting. We comment on this at the end of Sect. 3.1.1.

To pass to the limit in the approximate of the nonlinear term ${\widetilde{\Psi }}'_{0}(\phi )$, where ${\widetilde{\Psi }}_{0}(\cdot )$ corresponds to the convex part of $\Psi (\cdot )$ and is defined on entire $\mathbb {R},$ cf. (2.2), one first needs to show an apriori estimate of the same in $L^1(Q_{T})$ and next identify the weak limit for a non relabeled subsequence with ${\widetilde{\Psi }}'_{0}(\phi ).$ For the first part the idea roughly is to write the Cahn–Hilliard equation as

$$\begin{aligned} \begin{aligned} \displaystyle -\eta \Delta \phi ^{N}+\widetilde{{\Psi }}'_{0}(\phi ^{N})=&\text {lower order terms}+\frac{|\nabla M^{N}|^{2}}{2}{} & {} \text { in }\Omega \\ \displaystyle \partial _{n}\phi ^{N}=&0{} & {} \text { on }\partial \Omega \end{aligned} \end{aligned}$$

(1.22)

and to use the elliptic structure to obtain a suitable uniform bound for $\phi ^{N}$ and $\widetilde{{\Psi }}'_{0}(\phi ^{N}).$ In that direction one needs to estimate the right hand side of (1.22)$_{1}$ in $L^{q}(\Omega )$ with $q>1.$ The boundedness of $E_{tot}$ (defined in (1.17)) alone does not provide this information and hence we exploit the dissipative part of the energy and the equation solved by $M^{N}$ to obtain that

$$\begin{aligned} {{\,\textrm{div}\,}}\,({{\xi }(\phi ^{N})\nabla M^{N}})\in L^{2}(0,T;L^{\frac{3}{2}}(\Omega )). \end{aligned}$$

(1.23)

Next one would expect to recover an improved bound on $\{M^{N}\}$ from (1.23) but in view of the uniform bound on $\{\phi ^{N}\}$ in $L^{\infty }(0,T;W^{1,2}(\Omega ))$) one can only use the fact that $\{\xi (\phi ^{N})\}$ is bounded in $L^{\infty }(Q_{T})$ and nondegenerate. With this setup we can use [32], which deals with the regularity of weak solutions to elliptic problems with nondegenerate, bounded and measurable coefficients, to obtain a uniform estimate of $M^{N}$ in $L^{2}(0,T;W^{1,p}(\Omega ))$ for some $p>2,$ in fact p is slightly greater than two and tends to two as the operator in (1.23) degenerates. At this moment one can recover a uniform estimate of $|\nabla M^{N}|^{2}$ in $L^{2}(0,T;L^{\frac{2p}{p+2}}(\Omega ))$ since $\displaystyle \nabla M^{N}\in L^{\infty }(0,T;L^{2}(\Omega ))\cap L^{2}(0,T;L^{p}(\Omega ))$. This in turn allows us to use (1.22) and to obtain a uniform bound of ${\widetilde{\Psi }}'_{0}(\phi ^{N})$ in $L^{2}(0,T;L^{\frac{2p}{p+2}}(\Omega )).$ The details of obtaining these uniform bounds of $M^{N}$ and ${\widetilde{\Psi }}'_{0}(\phi ^{N})$ can be found in Sect. 3.1.2. To identify the limit of a subsequence of $\{{\widetilde{\Psi }}'_{0}(\phi ^{N})\}$ with ${\widetilde{\Psi }}'_{0}(\phi ),$ we adapt the ideas related to the estimate of the measure of the set $\{|\phi |=1\}$ developed in [20, 25].

After the recovery of the weak formulations (1.13) solved by $(v,M,\phi ,\mu )$ we prove that the obtained weak solution attains the initial data in a strong sense (cf. (1.14)).

Items (ii) and (iii) of Theorem 1.5, which correspond to the obtainment of boundary condition for M in a weak sense and some further regularity results of M in Lebesgue spaces, are proved in Sect. 3.4.

1.3 Physical background and comments on the derivation of the model

In this section we comment on the physical background of the model (1.1). In [36] we already have derived a model for diffuse interface magnetic fluids but with matched densities and a smooth double well potential for the mixing energy. In the present article the density consideration of the mixture is inspired from [10]. To deal with the density dependence it is important to consider a singular logarithmic potential for the mixing energy which is also more physical than a smooth double well potential considered in [36].

We assume that the mean velocity satisfies the homogeneous Dirichlet boundary condition on $\partial \Omega .$ For the derivation of a modified momentum balance equation solved by a solenoidal mean velocity field v and mean density $\rho (\phi )$ (cf. (1.9)) we refer to [10, Section 2]. The obtained momentum balance equation with a general stress tensor $\mathbb {S}$ is of the form:

$$\begin{aligned} \begin{array}{ll} \rho \partial _{t}v+\left( (\rho v+J)\cdot \nabla \right) v={{\,\textrm{div}\,}}\mathbb {S}\quad \text{ in }\,\, Q_{T}, \displaystyle \end{array} \end{aligned}$$

(1.24)

along with the incompressibility ${{\,\textrm{div}\,}}\,v=0$ where J is the relative diffusion flux defined in (1.2). We assume that the stress tensor $\mathbb {S}$ is the sum of the standard viscous Newtonian stress tensor $2\nu (\phi )\mathbb {D}(v)-\pi \mathbb {I}$ (where $\mathbb {I}$ is the identity matrix and $\pi$ being the mean pressure) and an extra contribution from the mixing energy and a simplified micromagnetic energy.

In order to derive the dependence of $\mathbb {S}$ on $\phi ,$ $\nabla \phi ,$ M and $\nabla M$ we start with the expression of ${\mathcal {E}}_{mix}$ (the mixing energy) and ${\mathcal {E}}_{mag}$ (the micromagnetic energy). The mixing energy reads as follows

$$\begin{aligned} \begin{aligned} {\mathcal {E}}_{mix}= \frac{\eta }{2}\int _{\Omega }|\nabla \phi |^{2}+\int _{\Omega }\Psi (\phi ), \end{aligned} \end{aligned}$$

(1.25)

where $\phi$ is the order parameter, $\eta >0$ denotes the thickness of the interface where the two fluids mix and $\Psi (\cdot )$ is the singular potential satisfying (1.8) (or one can in particular consider the expression (1.3)).

In the present article the magnetic energy contribution ${\mathcal {E}}_{mag}$ is inspired from micromagnetics, for details we refer to [21] and the references therein. We only consider the exchange energy contribution, which reflects the tendency of the magnetization to orient in one direction. We further consider the dependence of the micromagnetic energy on the order parameter $\phi$ which in turn allows us to study a system involving fluids with different magnetic behavior. The simplified micromagnetic energy reads as follows

$$\begin{aligned} \begin{aligned}&{\mathcal {E}}_{mag}=\int _{\Omega }\xi (\phi )\frac{|\nabla M|^{2}}{2}+\frac{1}{4\alpha ^{2}}\int _{\Omega }\xi (\phi )(|M|^{2}-1)^{2}, \end{aligned} \end{aligned}$$

where $\alpha >0$ is a parameter. The first term in the expression of ${\mathcal {E}}_{mag}$ is the exchange energy contribution and the second one is a penalization term punishing the derivation of |M| from one (which is a more physical constraint). The consideration of such a penalization term is standard in the literature, cf. [37, Section 1.2], [19] or [44]. The magnetic energy in our case is coupled with the order parameter via a regular, bounded and non degenerate function $\xi$ (we refer to (1.6) for the assumptions on $\xi$). In a little different situation, for the modeling and numerical analysis of liquid crystals, one can find specific expressions of degenerate functions $\xi (\cdot )$ in the articles [46] and [51]. For us it seems important to choose a non degenerate function $\xi$ which in turn plays a crucial role to obtain the strong compactness of $\nabla M.$

Now exactly as in [51, Section 2.3] one can use the principle of virtual work to compute that the contribution of ${\mathcal {E}}_{mix}$ and ${\mathcal {E}}_{mag}$ to the stress tensor $\mathbb {S}$ is given by

$$\begin{aligned} - \frac{\partial {\mathcal {E}}_{mag}}{\partial \nabla M} \odot \nabla M -\frac{\partial {\mathcal {E}}_{mix}}{\partial \nabla \phi } \otimes \nabla \phi = -\xi (\phi )(\nabla M\odot \nabla M)-\eta (\nabla \phi \otimes \nabla \phi ), \end{aligned}$$

where $(\nabla M\odot \nabla M)_{ij}=\sum _{k=1}^{3}(\nabla _{i}M_{k})(\nabla _{j}M_{k})$ and $(\nabla \phi \otimes \nabla \phi )_{ij}=\nabla _{i}\phi \nabla _{j}\phi$.

Hence altogether $\mathbb {S}$ has the following expression

$$\begin{aligned} \begin{array}{l} \mathbb {S}=2\nu (\phi )\mathbb {D}(v)-\xi (\phi )(\nabla M\odot \nabla M)-\eta (\nabla \phi \otimes \nabla \phi ). \end{array} \end{aligned}$$

(1.26)

The momentum balance (1.24) along with (1.26) reads as follows:

$$\begin{aligned} \begin{aligned} \rho \partial _{t}v&+\left( (\rho v+J)\cdot \nabla \right) v-{{\,\textrm{div}\,}}(2\nu (\phi )\mathbb {D}(v))+\nabla \pi \\ =&-{{\,\textrm{div}\,}}\bigl (\xi (\phi )(\nabla M\odot \nabla M)+\eta (\nabla \phi \otimes \nabla \phi )\bigr )\text { in } Q_{T}, \end{aligned} \end{aligned}$$

(1.27)

with ${{\,\textrm{div}\,}}\,v=0$ and $v=0$ on $\Sigma _{T}.$

Next we assume physically reasonable boundary conditions $\partial _{n}M=\partial _{n}\phi =\partial _{n}\mu =0$ on $\Sigma _{T}.$ The derivation of the magnetization equation (1.1)$_{3}$ is based on gradient flow dynamics. The obtainment of the Cahn–Hilliard equations (1.1)$_{4,5}$ relies on the generalized Fick’s law, i.e., the mass flux be proportional to the gradient of the chemical potential (we refer to [17, 45] for details). The detailed derivation of (1.1)$_{3,4,5}$ can be done by following the arguments presented in [36, p. 8], with modifications since here we use a singular potential in the mixing energy.

In view of (1.1)$_{4}$ one at once derives the following mass conservation

$$\begin{aligned} \begin{array}{l} \partial _{t}\rho +{{\,\textrm{div}\,}}\left( \rho v+J\right) =0\quad \text{ in }\,\,Q_{T}. \end{array} \end{aligned}$$

(1.28)

In the spirit of [10], we explain here the dynamics behind (1.28). The equation (1.28) implies that the flux of the density consists of two parts: $\rho v$, describing the transport by the mean velocity, and a relative flux J (cf. (1.2)) related to diffusion of the components. Hence for the unmatched density case, diffusion of the components leads to the diffusion of the mass density.

The modification of the momentum balance (cf. (1.24) and (1.27)) by adding the relative diffusion flux J was proposed in [10] to obtain a local dissipation inequality and global energy estimate for their model. It serves the same purpose for our case and we recover the following formal energy balance for the system (1.27)–(1.1)$_{2,3,4,5,6,7}:$

$$\begin{aligned} \begin{aligned}&\frac{d}{dt}E_{tot}+\Biggl (\Vert \sqrt{2\nu (\phi )}\mathbb {D} v\Vert ^2_{L^2(\Omega )}+\Vert \nabla \mu \Vert ^2_{L^2(\Omega )} \\&\quad +\left\| {{\,\textrm{div}\,}}(\xi (\phi )\nabla M)-\frac{\xi (\phi )}{\alpha ^2}M(|M|^2-1)\right\| ^2_{L^2(\Omega )}\Biggr )=0, \end{aligned} \end{aligned}$$

where $E_{tot}$ is given by (1.17).

Finally in view of the mass balance (1.28), the incompressibility of v and the identities

$$\begin{aligned} \begin{aligned}&\eta {{\,\textrm{div}\,}}(\nabla \phi \otimes \nabla \phi )=\eta \Delta \phi \nabla \phi +\eta \nabla \left( \frac{|\nabla \phi |^{2}}{2}\right) ,\\&{{\,\textrm{div}\,}}(\xi (\phi )(\nabla M\odot \nabla M))={{\,\textrm{div}\,}}(\xi (\phi )\nabla M)\nabla M+\xi (\phi )\nabla \left( \frac{|\nabla M|^{2}}{2}\right) , \end{aligned} \end{aligned}$$

(1.29)

one can rewrite (1.27) (we refer to [36, p. 8] for the use of the identities (1.29)) in the form (1.1)$_{1},$ of course with a modified pressure p. We emphasize that the term $\mu \nabla \phi$ has better compactness properties compared to ${{\,\textrm{div}\,}}(\nabla \phi \otimes \nabla \phi )$ and this will be exploited in the following analysis.

1.4 Bibliographical remarks

Diffuse interface models without magnetization and comprising of fluids with matched densities date back to the works [12, 33, 34, 47]. The article [33] provides a unified framework for coupled Navier–Stokes and Cahn–Hilliard equations using the balance law for microforces with a mechanical version of the second law of thermodynamics. The mathematical analysis of such a model first appeared in [47], where the author deals with strong solutions and stability of stationary solutions (as $t\rightarrow \infty$) in the setting $\Omega =\mathbb {R}^{2}$ assuming a smooth double well potential for the mixing energy. For a detailed review of the subject we also refer to [12].

A detailed study is then performed in the article [15] where the author proves the existence of global weak solution of the model both in dimension two and three in a channel with a smooth double well potential for the mixing energy. The article [15] further proves that the model under consideration (with non degenerate mobility) admits a strong solution which is global (in time) in dimension two and local (in time) in dimension three. The case of degeneracy in the Cahn–Hilliard equation is also considered in [15]. A more complete mathematical description (existence, uniqueness, regularity of solutions and asymptotic behavior) of a similar model with the physically relevant logarithmic potential is discussed in [2]. The article [2] proves the existence of strong solutions for an initial velocity $v_{0}$ in interpolation spaces between $W^{1,2}_{0}(\Omega )$ and $W^{2,2}(\Omega )\cap W^{1,2}_{0}(\Omega )$ and satisfying $\text{ div }\,v_{0}=0$ on $\Omega .$ Further the author shows that any weak solution of the system becomes regular for large times and the order parameter converges to a solution of the stationary Cahn–Hilliard system/a critical point of the mixing energy and the mean velocity tends to zero. We would also like to quote the article [31] where the authors study the uniqueness and regularity of weak and strong solutions of the model with the logarithmic potential. The article proves the existence and uniqueness of a global strong solution in dimension two under the assumption that the initial velocity $v_{0}\in W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )$ and the local in time existence and uniqueness of strong solution in dimension three. Unlike [2, 31] we study the case where the fluids involved have different density and they show magnetic behavior. The coupling of the magnetization with the Cahn–Hilliard dynamics do not allow us to obtain sufficient regularity of the unknowns to prove the uniqueness of weak solutions in dimension two.

In the literature there has been several approaches modeling diffuse interface systems (without magnetization) where the density is not constant. The article [39] derived a quasi-incompressible diffuse interface model where the velocity field is not divergence free. Analytical results for the model considered in [39] first appeared in [3, 4]. We also refer to [16, 22] for other diffuse interface models involving fluids with unmatched densities. For a slightly non homogeneous case, i.e., with the assumption that the densities of the fluids undergoing partial mixing tend to be equal, [16] proved the global existence of a weak solution and of a unique local strong solution (which is global in two dimension) for their model. The non homogeneous magnetic fluid dynamics we consider in the present article is inspired from the model introduced in [10]. The global existence of weak solutions for the model introduced in [10] can be found in [6] (the case of non degenerate mobility) and [7] (the case of degenerate mobility). In the recent article [29] the author proves the local in time existence of a strong solution in a bounded domain in dimension two of the model introduced in [10]. The author also shows the global existence in the space periodic set-up. Other variants of the Abels, Garcke and Grün model [10] can be found in [23, 24, 26]. The article [26] considers a general diffuse interface model for incompressible two-phase flows with unmatched densities describing the evolution of free interfaces in contact with the solid boundary whereas [23, 24] deal with a non-local version of the model derived in [10]. We would further like to quote some very recent articles viz. [9, 27, 30] on the analysis of the celebrated Abels, Garcke and Grün model. The article [30] deals with the analysis of strong solutions to the ’AGG’ model in three dimensions whereas global regularity and asymptotic stabilization is investigated in [9]. The authors of [27] analyze the global well-posedness and convergence to equilibrium for the ’AGG’ model with nonlocal free energy.

For other analytical results on varying density diffuse interface models we also refer the readers to [8, 18] (compressible fluids), [5] (non Newtonian fluids), [35] (coupling between Allen–Cahn and Navier–Stokes) and [38] (non-isothermal diffuse interface model).

Concerning the analysis of single phase magnetic fluids the readers can consult [13] (magnetization is modeled by LLG dynamics), [44] (gradient flow of magnetization) and for various models involving diffuse interface magnetic fluids with matched densities and smooth double well potential we refer to [36, 42, 50].

A crucial idea in our article is to recover the hidden $L^{2}(0,T;W^{1,p})$ ($p>2$) regularity of $M^{N},$ uniform with respect to N (where $M^{N}$ is the approximate of M defined in Sect. 3) by using the fact that ${{\,\textrm{div}\,}}(\xi (\phi ^{N})\nabla M^{N})$ is uniformly bounded in $L^{2}(0,T;L^{\frac{3}{2}}(\Omega ))$ and $\partial _{n}M^{N}\mid _{\partial \Omega }=0.$ In view of the regularity of $\phi ^{N}$ we can only use that $\xi (\phi ^{N})$ is bounded, measurable and non degenerate uniformly in N. To the best of our knowledge for such an elliptic boundary value problem there are two classes of result. One deals with the Hölder regularity ($C^{0,\gamma }(\Omega )$) of solutions and the second one proving $W^{1,p}(\Omega )$ ($p>2$) regularity. None of these results imply one another. For the first type of result we refer to [28, 41] whereas for the second we quote [32, 40]. In the present article we have used in particular the result proved in [32].

2 Existence of weak solutions to a time discrete model

In this section we prove the existence of weak solution to a time discrete problem corresponding to system (1.1). In that direction we need some regularity results for the Cahn–Hilliard equations proved in [1]. In order to access the regularity results from [1] one needs to reformulate the equation (1.1)$_{5}$ by using the subdifferential of a convex potential. This will be done in the next section and the arguments are inspired from [6].

2.1 Reformulation of the problem using subdifferential of a convex potential

First we define a potential ${\widetilde{\Psi }}$ as follows:

$$\begin{aligned} {\widetilde{\Psi }}:\mathbb {R}\longrightarrow \mathbb {R},\qquad {\widetilde{\Psi }}(s)={\left\{ \begin{array}{ll} \Psi (s)&{}\text { if }s\in [-1,1],\\ +\infty &{}\text { else}\end{array}\right. } \end{aligned}$$

(2.1)

where $\Psi$ is introduced in (1.8). Since $\Psi ''(s)\geqslant -\kappa ,\,\,\text{ for } \text{ some }\,\,\kappa \in \mathbb {R},$ ${\widetilde{\Psi }}$ is not necessarily convex. In order to use the theory of subdifferentials we introduce a convex function

$$\begin{aligned} {\widetilde{\Psi }}_{0}(r)={\widetilde{\Psi }}(r)+\frac{\kappa }{2}r^{2}\,\,\text{ for }\,\,r\in [-1,1],\,\,{\widetilde{\Psi }}_{0}\in C([-1,1])\cap C^{2}((-1,1)). \end{aligned}$$

(2.2)

With this new function ${\widetilde{\Psi }}_{0},$ the equation (1.1)$_{5}$ can be equivalently written as

$$\begin{aligned} \begin{array}{l} \displaystyle \mu +{\kappa }\phi ={\widetilde{\Psi }}'_{0}(\phi )-\eta \Delta \phi +\xi '(\phi )\frac{|\nabla M|^{2}}{2}+\frac{\xi '(\phi )}{4\alpha ^{2}}(|M|^{2}-1)^{2} \text{ in } Q_{T}. \end{array} \end{aligned}$$

(2.3)

Inspired by [6], we define the energy ${\widetilde{E}}:L^{2}(\Omega )\longrightarrow \mathbb {R}\cup \{+\infty \}$ with the domain

$$\begin{aligned} \begin{array}{l} \text{ dom }\,{\widetilde{E}}=\{\phi \in W^{1,2}(\Omega )\;|\;-1\leqslant \phi \leqslant 1\,\,\text{ a.e. }\} \end{array} \end{aligned}$$

(2.4)

as

$$\begin{aligned} {\widetilde{E}}(\phi )={\left\{ \begin{array}{ll} \displaystyle \frac{\eta }{2}\int \nolimits _{\Omega }|\nabla \phi |^{2}+\int \nolimits _{\Omega }{\widetilde{\Psi }}_{0}(\phi )&{}\text { for }\phi \in \text{ dom }\,{\widetilde{E}},\\ +\infty &{}\text { else.} \end{array}\right. } \end{aligned}$$

(2.5)

Proposition 2.1

[1, Theorem 3.12.8] Let ${\widetilde{E}}:L^{2}(\Omega )\longrightarrow \mathbb {R}\cup \{+\infty \}$ be as defined in (2.4)–(2.5). Then $\partial {\widetilde{E}}(\phi )=-\eta \Delta \phi +{\widetilde{\Psi }}'_{0}(\phi )$ and

$$\begin{aligned} \begin{array}{l} \displaystyle {\mathcal {D}}(\partial {\widetilde{E}})=\{\phi \in W^{2,2}_{n}(\Omega )\;|\;{\widetilde{\Psi }}'_{0}(\phi )\in L^{2}(\Omega ),\,\,{\widetilde{\Psi }}''_{0}(\phi )|\nabla \phi |^{2}\in L^{1}(\Omega )\} \end{array} \end{aligned}$$

(2.6)

is the domain of definition of the subgradient $\partial {\widetilde{E}}.$ Moreover, there exists a positive constant C such that

$$\begin{aligned} \begin{aligned}&\Vert \phi \Vert ^{2}_{W^{2,2}(\Omega )}+\Vert {\widetilde{\Psi }}'_{0}(\phi )\Vert ^{2}_{L^{2}(\Omega )}+\int _{\Omega }{\widetilde{\Psi }}''_{0}(\phi (\cdot ))|\nabla \phi (\cdot )|^{2}\\&\quad \leqslant C\left( \Vert \partial {\widetilde{E}}(\phi )\Vert ^{2}_{L^{2}(\Omega )}+\Vert \phi \Vert ^{2}_{L^{2}(\Omega )}+1\right) . \end{aligned} \end{aligned}$$

(2.7)

Further for every $1<p\leqslant 2$, there exists a constant $C_{p}>0$ such that

$$\begin{aligned} \begin{array}{l} \Vert \phi \Vert _{W^{2,p}(\Omega )}+\Vert {\widetilde{\Psi }}'_{0}(\phi )\Vert _{L^{p}(\Omega )}\leqslant C_{p}\left( \Vert \partial {\widetilde{E}}(\phi )\Vert _{L^{p}(\Omega )}+\Vert \phi \Vert _{L^{2}(\Omega )}+1\right) . \end{array} \end{aligned}$$

(2.8)

Remark 2.2

It follows from (1.8) and the definitions (2.1) and (2.2) that if $\phi \in {\mathcal {D}}(\partial {\widetilde{E}}),$ then $|\phi |\leqslant 1$ a.e.

With the help of the subgradient $\partial {\widetilde{E}}$ the equation (2.3) can be written as

$$\begin{aligned} \begin{array}{l} \displaystyle \mu +{\kappa }\phi =\partial {\widetilde{E}}(\phi )+\xi '(\phi )\frac{|\nabla M|^{2}}{2}+\frac{\xi '(\phi )}{4\alpha ^{2}}(|M|^{2}-1)^{2} \text{ in } Q_{T}. \end{array} \end{aligned}$$

(2.9)

It is interesting to note that for $\phi \in \text{ dom }\,{\widetilde{E}}$ the free energy corresponding to the system (1.1) is related to ${\widetilde{E}}(\phi )$ via the following relation

$$\begin{aligned} \begin{array}{ll} \displaystyle {\mathcal {E}}_{free}&\displaystyle =\frac{1}{2}\int _{\Omega }\xi (\phi )|\nabla M|^{2}+\frac{1}{4\alpha ^{2}}\int _{\Omega }\xi (\phi )(|M|^{2}-1)^{2}+{\widetilde{E}}(\phi )-\frac{\kappa }{2}\Vert \phi \Vert ^{2}_{L^{2}(\Omega )}. \end{array} \end{aligned}$$

2.2 Analysis of a time discrete model

To begin with we define a suitable time discretization of the model (1.1) keeping in mind the reformulation (2.3) (or (2.9)) of (1.1)$_{5}.$

Let $h>0$ be a constant,

$$\begin{aligned} v_{k}\in L^{2}_{{{\,\textrm{div}\,}}}(\Omega ),\,\,M_{k}\in W^{2,2}_{n}(\Omega ),\,\,\phi _{k}\in {\mathcal {D}}(\partial {\widetilde{E}}) \end{aligned}$$

(2.10)

with ${\mathcal {D}}(\partial {\widetilde{E}})$ as in (2.6) and

$$\begin{aligned} \begin{array}{l} \rho _{k}=\frac{1}{2}({\widetilde{\rho }}_{1}+{\widetilde{\rho }}_{2})+\frac{1}{2}({\widetilde{\rho }}_{2}-{\widetilde{\rho }}_{1})\phi _{k} \end{array} \end{aligned}$$

(2.11)

be the information at time step $t_{k},$ $k\in \mathbb {N}_{0}$. The quadruple $(v_{k+1},M_{k+1},\phi _{k+1}$, $\mu _{k+1}),$ solution at the time step $t_{k+1},$ is determined as a weak solution to the following system

$$\begin{aligned} \begin{aligned} \frac{\rho _{k+1}v_{k+1}-\rho _{k}v_{k}}{h}+{{\,\textrm{div}\,}}(\rho _{k}v_{k+1}\otimes v_{k+1}&)+\nabla p_{k+1}-\mu _{k+1}\nabla \phi _{k}\\ +{{\,\textrm{div}\,}}(v_{k+1}\otimes J_{k+1})-{{\,\textrm{div}\,}}(2\nu (\phi _{k})\mathbb {D}v_{k+1}&)\\ -\frac{{\xi (\phi _{k})}}{\alpha ^{2}}(|M_{k+1}|^{2}M_{k+1}\!-\!M_{k})\nabla M_{k{+}1} {=}&\!-\!{{\,\textrm{div}\,}}(\xi (\phi _{k})\nabla M_{k{+}1})\nabla M_{k{+}1}{} & {} \text { in }\Omega \\ {{\,\textrm{div}\,}}v_{k{+}1}{=}&0{} & {} \text { in }\Omega \\ \frac{M_{k+1}-M_{k}}{h}+(v_{k+1}\cdot \nabla )M_{k+1}=&{{\,\textrm{div}\,}}(\xi (\phi _{k})\nabla M_{k+1})\\ -&\frac{{\xi (\phi _{k})}}{\alpha ^{2}}(|M_{k+1}|^{2}M_{k+1}\!-\!M_{k}){} & {} \text { in }\Omega \\ \frac{\phi _{k+1}-\phi _{k}}{h}+(v_{k+1}\cdot \nabla )\phi _{k}=&\Delta \mu _{k+1}{} & {} \text { in }\Omega \\ \mu _{k+1}+\kappa \frac{\phi _{k+1}+\phi _{k}}{2}+\eta \Delta \phi _{k+1}-{\widetilde{\Psi }}'_{0}&(\phi _{k+1})\\ - \frac{1}{4\alpha ^{2}}{H_{0}(\phi _{k+1},\phi _{k})}(|M_{k+1}|^{2}-1)^{2} =&H_{0}(\phi _{k+1},\phi _{k})\frac{|\nabla M_{k+1}|^{2}}{2}{} & {} \text { in }\Omega \\ v_{k+1}=0,\ \partial _{n}M_{k+1}=0,\ \partial _{n}\phi _{k+1}=&\partial _{n}\mu _{k+1}=0{} & {} \text { on }\partial \Omega , \end{aligned} \end{aligned}$$

(2.12)

where

$$\begin{aligned} \begin{array}{l} \displaystyle J=J_{k+1}=-\frac{{\widetilde{\rho }}_{2}-{\widetilde{\rho }}_{1}}{2}\nabla \mu _{k+1},\,\,\,\, \rho _{k+1}=\frac{1}{2}({\widetilde{\rho }}_{1}+{\widetilde{\rho }}_{2})+\frac{1}{2}({\widetilde{\rho }}_{2}-{\widetilde{\rho }}_{1})\phi _{k+1} \end{array}\nonumber \\ \end{aligned}$$

(2.13)

and $H_{0}:\mathbb {R}\times \mathbb {R}\rightarrow \mathbb {R}$ is defined as

$$\begin{aligned} \ H_{0}(a,b)={\left\{ \begin{array}{ll} \frac{\xi (a)-\xi (b)}{a-b} &{}\text{ if }\qquad a\ne b,\\ \xi '(b) &{} \text{ if }\qquad a= b. \end{array}\right. } \end{aligned}$$

(2.14)

Now let us introduce the notion of weak solution to the time discrete system (2.12). In the following definition of weak solution, the term $\displaystyle \int _{\Omega }(\textrm{div}(v_{k+1}\otimes J_{k+1})){\widetilde{\psi }}_{1}$ (originated from the fifth term of (2.12)$_{1}$) is replaced using the identity

$$\begin{aligned} \begin{aligned}&\int _{\Omega }({{\,\textrm{div}\,}}(v_{k+1}\otimes J_{k+1})){\widetilde{\psi }}_{1}\\&\quad =\int _{\Omega }\left( \textrm{div}J_{k+1} -\frac{\rho _{k+1}\!-\!\rho _{k}}{h}\!-\!v_{k+1}\cdot \nabla \rho _{k}\right) \frac{v_{k+1}}{2}\cdot {\widetilde{\psi }}_{1}\\&\qquad +\int _{\Omega }(J_{k+1}\cdot \nabla )v_{k+1}\cdot {\widetilde{\psi }}_{1} \end{aligned} \end{aligned}$$

(2.15)

(we refer to (2.17)) and the derivation can be found in [6, Remark 4.1 (i)]. This reformulation helps mainly to obtain later the discrete energy estimate (2.29).

Definition 2.3

(Weak solution to the problem (2.12)) Let (2.10)–(2.11) hold. The quadruple

$$\begin{aligned} \begin{array}{l} (v_{k+1},M_{k+1},\phi _{k+1},\mu _{k+1})\in W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )\times W^{2,2}_{n}(\Omega )\times {\mathcal {D}}(\partial {{\widetilde{E}}})\times W^{2,2}_{n}(\Omega ), \end{array}\nonumber \\ \end{aligned}$$

(2.16)

is a weak solution to system (2.12) if the following identities are true

$$\begin{aligned} \begin{aligned}&\int _{\Omega }\frac{\rho _{k+1}v_{k+1}-\rho _kv_{k}}{h}\cdot {{\widetilde{\psi }}}_{1}+\int _{\Omega }{{\,\textrm{div}\,}}(\rho _{k}v_{k+1}\otimes v_{k+1})\cdot {\widetilde{\psi }}_{1}\\&\qquad +\int _{\Omega }\left( {{\,\textrm{div}\,}}J_{k+1} -\frac{\rho _{k+1}-\rho _{k}}{h}-v_{k+1}\cdot \nabla \rho _{k}\right) \frac{v_{k+1}}{2}\cdot {\widetilde{\psi }}_{1}\\&\qquad {+}\int _{\Omega }(J_{k+1}\cdot \nabla )v_{k+1}\cdot {\widetilde{\psi }}_{1}-\int _{\Omega }\left( \frac{{\xi (\phi _{k})}}{\alpha ^{2}}(|M_{k+1}|^{2}M_{k+1}-M_{k})\nabla M_{k+1}\right) \cdot {{\widetilde{\psi }}}_{1}\\&\qquad +\int _{\Omega }\left( \textrm{div}(\xi (\phi _{k})\nabla M_{k+1})\nabla M_{k+1}\right) \cdot {{\widetilde{\psi }}}_{1}\\&\quad =-2\int _{\Omega }\nu (\phi _{k})\mathbb {D} v_{k+1}\cdot \mathbb {D}{{\widetilde{\psi }}}_{1}-\int _{\Omega }\nabla \mu _{k+1}\phi _{k}\cdot {{\widetilde{\psi }}}_{1} \end{aligned}\nonumber \\ \end{aligned}$$

(2.17)

for all ${{\widetilde{\psi }}}_{1}\in W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ),$

$$\begin{aligned} \begin{aligned}&\int _{\Omega }\frac{M_{k+1}-M_{k}}{h}\cdot {{\widetilde{\psi }}}_{2}+\int _{\Omega }(v_{k+1}\cdot \nabla )M_{k+1}\cdot {{\widetilde{\psi }}}_{2}\\&\quad =\int _{\Omega }\left( {{\,\textrm{div}\,}}\left( \xi (\phi _{k})\nabla M_{k+1}\right) -\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M_{k+1}|^{2}M_{k+1}-M_{k})\right) \cdot {{\widetilde{\psi }}}_{2} \end{aligned} \end{aligned}$$

(2.18)

for all ${{\widetilde{\psi }}}_{2}\in L^2(\Omega ),$

$$\begin{aligned} \begin{aligned}&\frac{\phi _{k+1}-\phi _{k}}{h}+(v_{k+1}\cdot \nabla )\phi _{k}=\Delta \mu _{k+1} \end{aligned} \end{aligned}$$

(2.19)

a.e. in $\Omega$ and

$$\begin{aligned} \begin{aligned}&\mu _{k+1}+\kappa \frac{\phi _{k+1}+\phi _{k}}{2}\!-\!H_{0}(\phi _{k+1},\phi _{k})\frac{|\nabla M_{k+1}|^{2}}{2}\!-\!\frac{H_{0}(\phi _{k+1},\phi _{k})}{4\alpha ^{2}}(|M_{k+1}|^{2}\!-\!1)^{2}\\&\quad =-\eta \Delta \phi _{k+1}+{\widetilde{\Psi }}'_{0}(\phi _{k+1}) \end{aligned}\nonumber \\ \end{aligned}$$

(2.20)

a.e. in $\Omega ,$ where $J_{k+1}$ and $\rho _{k+1}$ are as defined in (2.13) and $H_{0}$ is defines in (2.14).

In the next lemma we first prove an estimate of the $L^{1}$ norm of ${\widetilde{\Psi }}^{'}_{0}(\phi _{k+1})$ assuming the existence of a triplet $(M_{k+1},\phi _{k+1},\mu _{k+1})$ solving (2.20). Then using the obtained estimate of $\Vert {\widetilde{\Psi }}^{'}_{0}(\phi _{k+1})\Vert _{L^{1}(\Omega )}$ we further prove an estimate of $\displaystyle \left| \int _{\Omega }\mu _{k+1}\right| .$ This estimate will be specifically used in Sect. 3 to show (3.13). We will also use a similar estimate while showing (2.66).

Lemma 2.4

Let $\phi =\phi _{k+1}\in {\mathcal {D}}(\partial {\widetilde{E}})$ and $\mu =\mu _{k+1}\in W^{1,2}(\Omega )$ solve (2.20) with $\phi _{k}\in W^{2,2}_{n}(\Omega ),$ $|\phi _{k}|\leqslant 1$ in $\Omega$ such that

$$\begin{aligned} \frac{1}{|\Omega |}\int _{\Omega }\phi _{k}=\frac{1}{|\Omega |}\int _{\Omega }\phi \in (-1,1), \end{aligned}$$

and $M=M_{k+1}\in W^{2,2}_{n}(\Omega ).$ Then there exists a constant $\displaystyle C=C\left( \int _{\Omega }\phi _{k}\right) >0,$ such that

$$\begin{aligned} \begin{aligned} \Vert {\widetilde{\Psi }}'_{0}(\phi )\Vert _{L^{1}(\Omega )}+\left| \int _{\Omega }\mu \right| \leqslant&C\biggl (\Vert \nabla \mu \Vert _{L^{2}(\Omega )}+\Vert \nabla \phi \Vert ^{2}_{L^{2}(\Omega )} \\ {}&+\Vert M\Vert ^{2}_{W^{1,2}(\Omega )}+\Vert M\Vert ^{4}_{W^{1,2}(\Omega )}+1\biggr ). \end{aligned} \end{aligned}$$

(2.21)

Proof

Without the magnetization vector M a similar result was shown in [6, Lemma 4.2]. We will suitably adapt the line of arguments used in proving [6, Lemma 4.2] in our case.

We test (2.20) by $(\phi -\overline{\phi }),$ where $\displaystyle \overline{\phi }=\frac{1}{|\Omega |}\int _{\Omega }\phi$ and obtain

$$\begin{aligned} \begin{aligned}&\int _{\Omega }\mu (\phi -\overline{\phi })+\kappa \int _{\Omega }\frac{\phi +\phi _{k}}{2}(\phi -\overline{\phi }) -\int _{\Omega } H_{0}(\phi ,\phi _{k})\frac{|\nabla M|^{2}}{2}(\phi -\overline{\phi })\\&\qquad -\int _{\Omega }\frac{H_{0}(\phi ,\phi _{k})}{4\alpha ^{2}}(|M|^{2}-1)^{2}(\phi -\overline{\phi }) -\eta \int _{\Omega }\nabla \phi \cdot \nabla (\phi -\overline{\phi })\\&\quad =\int _{\Omega }{\widetilde{\Psi }}'_{0}(\phi )(\phi -\overline{\phi }). \end{aligned} \end{aligned}$$

(2.22)

One observes that $\displaystyle \int _{\Omega }\mu (\phi -\overline{\phi })=\int _{\Omega }(\mu -\overline{\mu })\phi ,$ where $\displaystyle \overline{\mu }=\frac{1}{|\Omega |}\int _{\Omega }\mu .$ Since $\overline{\phi }\in (-1+\epsilon ,1-\epsilon )$ for some $\epsilon >0$ (note that $\epsilon$ is independent of $\phi$) and $\lim \limits _{\phi \rightarrow \pm 1}{\widetilde{\Psi }}'_{0}(\phi )=\pm \infty ,$ one has the inequality

$$\begin{aligned} \begin{array}{l} {\widetilde{\Psi }}_{0}^{'}(\phi )(\phi -\overline{\phi })\geqslant C_{1}|{\widetilde{\Psi }}'_{0}(\phi )|-C_{2}, \end{array} \end{aligned}$$

(2.23)

for constants $C_{1}>0$ and $C_{2}.$ The inequality (2.23) can be proved by dividing $[-1,1]$ into three intervals $[-1,-1+\frac{\epsilon }{2}],$ $[-1+\frac{\epsilon }{2},1-\frac{\epsilon }{2}],$ $[1-\frac{\epsilon }{2},1],$ arguing by the blow up behavior of ${\widetilde{\Psi }}'_{0}$ at the endpoints $\{-1,1\}$ and the fact that ${\widetilde{\Psi }}'_{0}\in C([-1+\frac{\epsilon }{2},1-\frac{\epsilon }{2}]).$ Hence integrating (2.23) in $\Omega ,$ we have the estimate

$$\begin{aligned} \begin{array}{l} \displaystyle \int _{\Omega } {\widetilde{\Psi }}'_{0}(\phi )(\phi -\overline{\phi })\geqslant C_{1}\int _{\Omega }|\Psi '_{0}(\phi )|-C_{3} \end{array} \end{aligned}$$

(2.24)

for constants $C_{1}>0$ and $C_{3}$. Further using (1.6)$_{2}$ and the definition (2.14) of $H_{0}(\cdot ,\cdot )$ one has $|H_{0}(\phi _{k+1},\phi _{k})|\leqslant c_{3}.$ Hence in view of (2.22), (2.24) and the fact that $|\phi |,|\phi _{k}|\leqslant 1,$ we deduce

$$\begin{aligned} \begin{array}{ll} \displaystyle \int _{\Omega }|{\widetilde{\Psi }}'_{0}|&{}\displaystyle \leqslant C\left( \Vert \mu -\overline{\mu }\Vert _{L^{2}(\Omega )}+\Vert \nabla M\Vert _{L^{2}(\Omega )}^{2}+\Vert M\Vert ^{4}_{L^{4}(\Omega )}+\Vert \nabla \phi \Vert _{L^{2}(\Omega )}^{2}+1\right) \\ &{}\displaystyle \leqslant C\left( \Vert \nabla \mu \Vert _{L^{2}(\Omega )}+\Vert M\Vert ^{2}_{W^{1,2}(\Omega )}+\Vert M\Vert ^{4}_{W^{1,2}(\Omega )}+\Vert \nabla \phi \Vert _{L^{2}(\Omega )}^{2}+1\right) , \end{array}\nonumber \\ \end{aligned}$$

(2.25)

where we have used Poincaré’s inequality to obtain the final step.

Now we want to use the inequality (2.25) to obtain an estimate of $\displaystyle \left| \int _{\Omega }\mu \right| .$ In that direction we integrate (2.20) to obtain

$$\begin{aligned} \begin{array}{ll} \displaystyle \left| \int \nolimits _{\Omega } \mu \right|&\leqslant \displaystyle C\left( \int _{\Omega }|{\widetilde{\Psi }}'_{0}|+\Vert M\Vert ^{2}_{W^{1,2}(\Omega )}+\Vert M\Vert ^{4}_{W^{1,2}(\Omega )}+\Vert \nabla \phi \Vert _{L^{2}(\Omega )}^{2}+1\right) . \end{array}\nonumber \\ \end{aligned}$$

(2.26)

Next using (2.25) in (2.26) we furnish

$$\begin{aligned} \begin{array}{l} \displaystyle \left| \int \nolimits _{\Omega } \mu \right| \leqslant C\left( \Vert \nabla \mu \Vert _{L^{2}(\Omega )}+\Vert \nabla \phi \Vert ^{2}_{L^{2}(\Omega )}\displaystyle +\Vert M\Vert ^{2}_{W^{1,2}(\Omega )}+\Vert M\Vert ^{4}_{W^{1,2}(\Omega )}+1\right) . \end{array}\nonumber \\ \end{aligned}$$

(2.27)

Combining (2.25) and (2.27) we conclude the proof of Lemma 2.4. $\square$

Next we recall the following result from [36, Lemma 4.1], which will be used to obtain (2.29) (a discrete analogue of energy dissipation) in Theorem 2.6.

Lemma 2.5

[36, Lemma 4.1] Let $A,B\in \mathbb {R}^{3}.$ The following relation is true

$$\begin{aligned} \begin{aligned}&\frac{1}{4}\bigl (|A|^{2}-1\bigr )^{2}-\frac{1}{4}\bigl (|B|^{2}-1\bigr )^{2}+\frac{1}{4}\bigl (|A|^{2}-|B|^{2}\bigr )^{2}+\frac{1}{2}|A\cdot (A-B)|^{2} \\&\quad +\frac{1}{2}|A-B|^{2}\leqslant (A-B)\cdot \bigl (|A|^{2}A-B\bigr ). \end{aligned} \end{aligned}$$

(2.28)

Now we state and prove the central result of this section which corresponds to the existence of weak solution to the time discrete system (2.12).

Theorem 2.6

(Existence of weak solution to the problem (2.12)) Let Assumption 1.1, (2.10) and (2.11) hold. Then there exists a quadruple $(v_{k+1}$, $M_{k+1}, \phi _{k+1},\mu _{k+1})$ which satisfies (2.16) and solves the identities (2.17)–(2.20). Moreover, the following discrete version of the energy estimate holds

$$\begin{aligned} \begin{aligned}&E_{tot}(v_{k+1},M_{k+1},\phi _{k+1})+2h\int _{\Omega }\nu (\phi _{k})|\mathbb {D}v_{k+1}|^{2} +h\int _{\Omega }|\nabla \mu _{k+1}|^{2}\\&\qquad \!+\!h\int _{\Omega }\left| \textrm{div}(\xi (\phi _{k})\nabla M_{k+1})\!-\!\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M_{k+1}|^{2}M_{k+1}\!-\!M_{k})\right| ^{2}\\&\quad \leqslant E_{tot}(v_{k},M_{k},\phi _{k}), \end{aligned} \end{aligned}$$

(2.29)

where $E_{tot}(v,M,\phi )$ is as defined in (1.17).

Proof

For simplicity in notations we will omit the subscript $k+1$ and we use the notation $(v,M,\phi ,\mu ,J,\rho )=(v_{k+1},M_{k+1},\phi _{k+1},\mu _{k+1},J_{k+1},\rho _{k+1})$ in the rest of the proof. We will perform the proof in two steps (c.f. Sect. 2.3 and 2.4).

2.3 Any weak solution $(v, M, \phi , \mu )$ of (2.12) in the sense of Definition 2.3 satisfies (2.29)–(1.17)

In the following computations we will need some identities in the spirit of [6]. We gather those identities in the following and refer to the proof of [6, Lemma 4.3] for details.

$$\begin{aligned} \begin{array}{ll} &{}(i)\,\,\displaystyle \int _{\Omega }\left( (\text{ div }\,J)\frac{v}{2}+\left( J\cdot \nabla \right) v\right) \cdot v=\int _{\Omega }\text{ div }\,\left( J\frac{|v|^{2}}{2}\right) =0,\\ &{}(ii)\,\,\displaystyle \int _{\Omega } \left( \text{ div }(\rho _{k}v\otimes v)-(\nabla \rho _{k}\cdot v)\frac{v}{2}\right) \cdot v=0,\\ &{}(iii)\,\,\displaystyle (\rho v-\rho _{k}v_{k})\cdot v=\left( \rho \frac{|v|^{2}}{2}-\rho _{k}\frac{|v_{k}|^{2}}{2}\right) +(\rho -\rho _{k})\frac{|v|^{2}}{2}+\rho _{k}\frac{|v-v_{k}|^{2}}{2}. \end{array} \end{aligned}$$

(2.30)

First we consider the test function ${\widetilde{\psi }}_{1}=v$ in (2.17) and use the identities (2.30) to render

$$\begin{aligned} \begin{aligned}&\int _{\Omega } \frac{\rho |v|^{2}-\rho _{k}|v_{k}|^{2}}{2h} +\int _{\Omega }\rho _{k}\frac{|v-v_{k}|^{2}}{2h}-\int _{\Omega }\left( \frac{{\xi (\phi _{k})}}{\alpha ^{2}}(|M|^{2}M-M_{k})\nabla M\right) \cdot v\\&\quad +\int _{\Omega }\left( \textrm{div}(\xi (\phi _{k})\nabla M)\nabla M\right) \cdot v =-2\int _{\Omega }\nu (\phi _{k})|\mathbb {D} v|^{2}-\int _{\Omega }(v\cdot \nabla )\mu \phi _{k}. \end{aligned} \end{aligned}$$

(2.31)

Next choosing ${{\widetilde{\psi }}}_2=-\text{ div }(\xi (\phi _{k})\nabla M)+\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})$ in (2.18) we infer

$$\begin{aligned} \begin{array}{ll} &{}\displaystyle \int _{\Omega }\frac{(M-M_{k})}{h}\cdot \Bigl (-\text{ div }(\xi (\phi _{k})\nabla M) +\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})\Bigr )\\ &{}\quad \displaystyle +\int _{\Omega }(v\cdot \nabla )M\cdot \Bigl (-\text{ div }(\xi (\phi _{k})\nabla M) +\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})\Bigr )\\ &{}\quad \displaystyle +\int _{\Omega }\Bigl |\text{ div }(\xi (\phi _{k})\nabla M)-\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})\Bigr |^{2}=0. \end{array} \end{aligned}$$

(2.32)

Multiplying (2.19) by $\mu$ and integrating over $\Omega$ one has

$$\begin{aligned} \begin{array}{l} \displaystyle \int _{\Omega }\frac{\phi -\phi _{k}}{h}\mu -\int _{\Omega }(v\cdot \nabla )\mu \phi _{k}=-\int _{\Omega }|\nabla \mu |^{2}. \end{array} \end{aligned}$$

(2.33)

Finally, multiplying (2.20) by $\displaystyle -\frac{\phi -\phi _{k}}{h}$ and integrating in $\Omega$ we have

$$\begin{aligned} \begin{aligned}&-\int _{\Omega }\mu \frac{\phi -\phi _{k}}{h}-\kappa \int _{\Omega }\frac{\phi ^{2}-\phi _{k}^{2}}{2h}+\int _{\Omega } H_{0}(\phi ,\phi _{k})\frac{|\nabla M|^{2}}{2}\frac{(\phi -\phi _{k})}{h}\\&\qquad +\int _{\Omega }\frac{H_{0}(\phi ,\phi _{k})}{4\alpha ^{2}}(|M|^{2}-1)^{2}\frac{(\phi -\phi _{k})}{h} +\eta \int _{\Omega }\nabla \phi \cdot \nabla \frac{(\phi -\phi _{k})}{h}\\&\quad =-\int _{\Omega }{\widetilde{\Psi }}'_{0}(\phi )\frac{(\phi -\phi _{k})}{h}. \end{aligned} \end{aligned}$$

(2.34)

Adding the expressions (2.31)–(2.34) and recalling (2.14), we have

$$\begin{aligned} \begin{aligned}&\frac{1}{2}\int _{\Omega }(\rho |v|^{2}{-}\rho _{k}|v_{k}|^{2}) +\int _{\Omega }\rho _{k}\frac{|v-v_{k}|^{2}}{2} +2h\int _{\Omega }\nu (\phi _{k})|\mathbb {D} v|^{2} \\&\quad +\int _{\Omega }(M-M_{k})\cdot \Bigl (-\text{ div }(\xi (\phi _{k})\nabla M)+\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})\Bigr )\\&\quad {+}h\int _{\Omega }\Bigl |\text{ div }(\xi (\phi _{k})\nabla M){-}\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M{-}M_{k})\Bigr |^{2}\!{+}\!\frac{1}{2}\int _{\Omega }\bigl (\xi (\phi )\!{-}\!\xi (\phi _{k})\bigr )|\nabla M|^{2}\\&\quad +\frac{1}{4\alpha ^{2}}\int _{\Omega }\bigl (\xi (\phi )-\xi (\phi _{k})\bigr )(|M|^{2}-1)^{2}-\frac{\kappa }{2}\int _{\Omega }({\phi ^{2}-\phi _{k}^{2}})\\&\quad +\eta \int _{\Omega }\nabla \phi \cdot (\nabla \phi -\nabla \phi _{k})+\int _{\Omega }{\widetilde{\Psi }}_{0}'(\phi )(\phi -\phi _{k})+h\int _{\Omega }|\nabla \mu |^{2}=0. \end{aligned}\nonumber \\ \end{aligned}$$

(2.35)

Integrating by parts the fourth term of (2.35), using

$$\begin{aligned} \begin{aligned} \ A\cdot (A-B)=\frac{|A|^{2}}{2}-\frac{|B|^{2}}{2}+\frac{|A-B|^{2}}{2}\,\,\text { for all }\,\,A,\,B\in \mathbb {R}^{m},\,m\in \mathbb {N}, \end{aligned}\nonumber \\ \end{aligned}$$

(2.36)

to expand $\xi (\phi _{k})\nabla M\cdot (\nabla M-\nabla M_{k})$ and $\nabla \phi \cdot (\nabla \phi -\nabla \phi _{k})$ respectively and using Lemma 2.5 on the term $(M-M_{k})\cdot \frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})$, we render

$$\begin{aligned} \begin{aligned}&\frac{1}{2}\int _{\Omega }(\rho |v|^{2}\!-\!\rho _{k}|v_{k}|^{2}) \!+\!\int _{\Omega }\rho _{k}\frac{|v-v_{k}|^{2}}{2}\!+\!2h\int _{\Omega }\nu (\phi _{k})|\mathbb {D} v|^{2}\!+\!\frac{1}{2}\int _{\Omega }\xi (\phi )|\nabla M|^{2}\\&\qquad {-}\frac{1}{2}\int _{\Omega }\xi (\phi _{k})|\nabla M_{k}|^{2}\!{+}\!\frac{1}{2}\int _{\Omega }\xi (\phi _{k})|\nabla M{-}\nabla M_{k}|^{2}\!{+}\!\frac{1}{4\alpha ^{2}}\int _{\Omega }\xi (\phi )\bigl (|M|^{2}\!{-}\!1\bigr )^{2}\\&\qquad -\frac{1}{4\alpha ^{2}}\int _{\Omega }\xi (\phi _{k})\bigl (|M_{k}|^{2}-1\bigr )^{2} \qquad +\frac{1}{4\alpha ^{2}}\int _{\Omega }\xi (\phi _{k})\bigl (|M|^{2}-|M_{k}|^{2}\bigr )^{2}\\&\qquad +\frac{1}{2\alpha ^{2}}\int _{\Omega }\xi (\phi _{k})|M\cdot (M-M_{k})|^{2} +\frac{1}{2\alpha ^{2}}\int _{\Omega }\xi (\phi _{k})|M-M_{k}|^{2}\\&\qquad {+}h\int _{\Omega }\Bigl |\text{ div }(\xi (\phi _{k})\nabla M)\!{-}\!\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M\!{-}\!M_{k})\Bigr |^{2}\!{+}\!\frac{\eta }{2}\int _{\Omega }|\nabla \phi |^{2}\!{-}\!\frac{\eta }{2}\int _{\Omega }|\nabla \phi _{k}|^{2}\\&\qquad +\!\frac{\eta }{2}\int _{\Omega }|\nabla \phi \!-\!\nabla \phi _{k}|^{2}\!+\!\int _{\Omega }\left( {\widetilde{\Psi }}_{0}(\phi )\!-\!{\widetilde{\Psi }}_{0}(\phi _{k})\right) \!-\!\frac{\kappa }{2}\int _{\Omega }({\phi ^{2}\!-\!\phi _{k}^{2}})\!+\! h\int _{\Omega }|\nabla \mu |^{2}\\&\quad \leqslant 0, \end{aligned}\nonumber \\ \end{aligned}$$

(2.37)

where we have used

$$\begin{aligned} \ \int _{\Omega }{\widetilde{\Psi }}_{0}'(\phi )(\phi -\phi _{k})\geqslant \int _{\Omega }\left( {\widetilde{\Psi }}_{0}(\phi )-{\widetilde{\Psi }}_{0}(\phi _{k})\right) , \end{aligned}$$

(2.38)

which follows from the convexity of ${\widetilde{\Psi }}_{0}$. Dropping some positive terms from the left hand side of the inequality (2.37) and recalling (2.2), we conclude the obtainment of the discrete energy estimate (2.29).

2.4 Proof of the existence of weak solutions to (2.12)

We will apply the Leray-Schauder fixed point principle to prove the existence of a weak solution to the discretized system (2.12). We start by considering the following spaces

$$\begin{aligned} \begin{aligned}&X=W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )\times W^{2,2}_{n}(\Omega )\times {\mathcal {D}}(\partial {\widetilde{E}})\times W^{2,2}_{n}(\Omega ),\\&Y= \bigl (W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )\bigr )'\times L^2(\Omega )\times L^{2}(\Omega )\times L^{2}(\Omega ), \end{aligned} \end{aligned}$$

(2.39)

with the norm defined as the sum of the individual components of the Cartesian products. We will write (2.12) in operator notation and for that we introduce ${\mathcal {N}}_{k},{\mathcal {F}}_{k}:X\rightarrow Y$. For $w=(v,M,\phi ,\mu )\in X,$ the operator ${\mathcal {N}}_{k}$ is defined as follows

$$\begin{aligned} \begin{aligned} {\mathcal {N}}_{k}(w)=\begin{pmatrix} {\mathcal {A}}v\\ \displaystyle -{{\,\textrm{div}\,}}(\xi (\phi _{k})\nabla M)+\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})+\int _\Omega M\\ \partial {\widetilde{E}}(\phi )+\phi \\ \displaystyle -\Delta \mu +\int _\Omega \mu \end{pmatrix}, \end{aligned} \end{aligned}$$

(2.40)

where ${\mathcal {A}}: W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )\rightarrow \bigl (W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )\bigr )'$ is given for all $v\in W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )$ by

$$\begin{aligned} \langle {\mathcal {A}}v,{{\widetilde{\psi }}}_1\rangle&=2\int _{\Omega }\nu (\phi _{k})\mathbb {D} v\cdot \mathbb {D} {{\widetilde{\psi }}}_1\text { for all }{{\widetilde{\psi }}}_1\in W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ). \end{aligned}$$

For $w=(v,M,\phi ,\mu )\in X,$ the operator ${\mathcal {F}}_{k}$ is defined as

$$\begin{aligned} {\mathcal {F}}_{k}(w)\!=\!\begin{pmatrix} \displaystyle -\frac{\rho v\!-\!\rho _{k}v_{k}}{h}\!-\!\text{ div }(\rho _{k}v\otimes v)\!-\!\nabla \mu \phi _{k}\!+\!\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M\!-\!M_{k})\nabla M\\ \displaystyle -\!{{\,\textrm{div}\,}}(\xi (\phi _{k})\nabla M)\nabla M \!-\!\frac{1}{2}\left( {{\,\textrm{div}\,}}J\!-\!\frac{\rho \!-\!\rho _{k}}{h}\!-\!v\cdot \nabla \rho _{k}\right) v\!-\!\left( J\cdot \nabla \right) v \\ \displaystyle -\frac{M-M_{k}}{h}-(v\cdot \nabla )M+\int _\Omega M\\ \displaystyle \mu +\kappa \frac{(\phi +\phi _{k})}{2}- H_{0}(\phi ,\phi _{k})\frac{|\nabla M|^{2}}{2}-\frac{H_{0}(\phi ,\phi _{k})}{4\alpha ^{2}}(|M|^{2}-1)^{2}+{\phi }\\ \displaystyle \ -\frac{\phi -\phi _{k}}{h}-(v\cdot \nabla )\phi _{k}+\int _\Omega \mu \end{pmatrix}.\nonumber \\ \end{aligned}$$

(2.41)

One observes that $w=w_{k+1}(=(v_{k+1},M_{k+1},\phi _{k+1},\mu _{k+1})\in X)$ is a weak solution to the time discrete problem (2.12) iff the following holds

$$\begin{aligned} {\mathcal {N}}_{k}(w)={\mathcal {F}}_k(w)\,\,\text{ in }\,\, Y. \end{aligned}$$

To prove the existence of a solution to this operator equation we next show some properties of the operators ${\mathcal {N}}_{k}$ and ${\mathcal {F}}_{k}.$

2.4.1 Invertibility and continuity of the inverse of ${\mathcal {N}}_{k}$ between suitable spaces

Let us consider the operator ${\mathcal {N}}_{k}$ component wise. It is well known that ${\mathcal {A}}: W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )\rightarrow \bigl (W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )\bigr )'$ is bounded and invertible, and ${\mathcal {A}}^{-1}:\bigl (W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )\bigr )'\longrightarrow W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )$ is bounded and continuous. For a proof one can adapt the arguments we are going to use to prove similar issues for the second component of ${\mathcal {N}}_{k}$. We choose to present a detailed argument for the second component since it is more involved than the first one.

We remark that the second component of ${\mathcal {N}}_{k}$ defines a bounded operator from $W^{2,2}_{n}(\Omega )\longrightarrow L^{2}(\Omega ).$ As a first step towards proving that the second component of ${\mathcal {N}}_{k}$ admits of a bounded and continuous inverse from $L^{2}(\Omega )$ to $W^{2,2}_{n}(\Omega ),$ we first claim that the operator

$$\begin{aligned} \begin{aligned} {\mathcal {B}}_{k}(M)=&-{{\,\textrm{div}\,}}_{N}(\xi (\phi _{k})\nabla M)+\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})+\int _\Omega M\\&:W^{1,2}(\Omega )\longrightarrow (W^{1,2}(\Omega ))', \end{aligned} \end{aligned}$$

(2.42)

where ${{\,\textrm{div}\,}}_{N}$ is interpreted in the following weak sense

$$\begin{aligned} \langle -{{\,\textrm{div}\,}}_N \Phi ,\psi _2\rangle =\int _\Omega \Phi \cdot \nabla \psi _2\text { for all }\psi _2\in W^{1,2}(\Omega )\,\,\text{ and }\,\, \Phi \in L^{2}(\Omega ), \end{aligned}$$

is invertible and ${\mathcal {B}}_{k}^{-1}:(W^{1,2}(\Omega ))'\longrightarrow W^{1,2}(\Omega )$ is bounded and continuous. In order to show this, we consider an arbitrary couple $M_{1},\,\,M_{2}\in W^{1,2}(\Omega ).$ Then we compute the following duality product

$$\begin{aligned} \begin{aligned} \langle {\mathcal {B}}_{k}&(M_{1})-{\mathcal {B}}_{k}(M_{2}),M_{1}-M_{2}\rangle \\ =&\int _\Omega \xi (\phi _k)|\nabla (M_1-M_2)|^2 +\int _\Omega \frac{\xi (\phi _k)}{\alpha ^2}\bigl (|M_1|^2M_1-|M_2|^2M_2\bigr )\cdot (M_1-M_2)\\&+\left( \int _\Omega M_1-M_2\right) ^2=I_{1}+I_{2}+I_{3}. \end{aligned}\nonumber \\ \end{aligned}$$

(2.43)

One observes $I_1\geqslant c_1\Vert \nabla (M_1-M_2)\Vert ^2_{L^2(\Omega )}$ since $\xi$ is non degenerate (c.f. (1.6)). The monotonicity of $\alpha \mapsto |\alpha |^2\alpha$ implies $I_2\geqslant 0.$ Employing the lower bound on $\xi$ again we infer $I_1+I_3\geqslant c\Vert M_1-M_2\Vert ^2_{W^{1,2}(\Omega )}.$

Hence

$$\begin{aligned} \displaystyle \langle {\mathcal {B}}_{k}(M_{1})-{\mathcal {B}}_{k}(M_{2}),M_{1}-M_{2}\rangle \geqslant c\Vert M_{1}-M_{2}\Vert ^{2}_{W^{1,2}(\Omega )}, \end{aligned}$$

implying ${\mathcal {B}}_{k}:W^{1,2}(\Omega )\longrightarrow (W^{1,2}(\Omega ))'$ is strongly monotone.

Since $W^{1,2}(\Omega )\hookrightarrow L^{6}(\Omega ),$ one justifies the boundedness of ${\mathcal {B}}_{k}:W^{1,2}(\Omega )\longrightarrow (W^{1,2}(\Omega ))'.$ Now using the Lebesgue dominated convergence theorem one checks that ${\mathcal {B}}_k$ is radially continuous on $W^{1,2}(\Omega )$, i.e., for each pair $M,{{\widetilde{M}}}\in W^{1,2}(\Omega )$ the function $t\in \mathbb {R}\mapsto \langle {\mathcal {B}}_k(M+t{\widetilde{M}}),{{\widetilde{M}}}\rangle$ is continuous. It is not hard to check that $\langle {\mathcal {B}}_k(M),M\rangle \geqslant c\Vert M\Vert ^2_{W^{1,2}(\Omega )}-c_k$, for any $M\in W^{1,2}(\Omega )$ with $c_k$ depending on $M_k$. The latter inequality implies that ${\mathcal {B}}_k$ is coercive on $W^{1,2}(\Omega )$, i.e.,

$$\begin{aligned} \lim _{\Vert M\Vert _{W^{1,2}(\Omega )}\rightarrow \infty }\frac{\langle {\mathcal {B}}_k(M),M\rangle }{\Vert M\Vert _{W^{1,2}(\Omega )}}=\infty . \end{aligned}$$

Using [43, Theorem 2.14] we obtain the existence of the inverse operator ${\mathcal {B}}_k^{-1}:(W^{1,2}(\Omega ))'\rightarrow W^{1,2}(\Omega )$ that is bounded and Lipschitz continuous.

We now claim that

$$\begin{aligned} \begin{array}{l} {\mathcal {B}}^{-1}_{k}: L^{2}(\Omega )\longrightarrow W^{2,2}_{n}(\Omega )\,\,\text{ is } \text{ bounded } \text{ and } \text{ continuous. } \end{array} \end{aligned}$$

(2.44)

The proof of this claim can be performed by using a boot-strap argument and using the regularity results for the following set of equations

$$\begin{aligned} \begin{aligned} \Delta M=&\frac{1}{\xi (\phi _k)}\Bigl (F-\xi '(\phi _{k})\nabla M\nabla \phi _k+\frac{\xi (\phi _{k})}{\alpha ^{2}}\bigl (|M|^{2}M-M_{k}\bigr )\Bigr )\quad{} & {} \text { in }\Omega ,\\ \partial _{n}M=&0{} & {} \text { on }\partial \Omega , \end{aligned} \end{aligned}$$

for any $F\in L^{2}(\Omega )$, where we apply the fact that ${\mathcal {B}}_k^{-1}:(W^{1,2}(\Omega ))'\rightarrow W^{1,2}(\Omega )$ is bounded and continuous. We refer the readers to [36, Section 4.2.1, p 14-15] for a concrete proof. This proves our claim that the second component of ${\mathcal {N}}_{k}$ has a bounded and continuous inverse from $W^{2,2}_{n}(\Omega )$ and $L^{2}(\Omega ).$

Next we consider the third component of ${\mathcal {N}}_{k}.$ Recalling the definition of ${\widetilde{E}}$ from (2.4)–(2.5), it can be justified in view of Proposition 2.1 that

$$\begin{aligned} \partial {\widetilde{E}}+I:{\mathcal {D}}(\partial {\widetilde{E}})\longrightarrow L^{2}(\Omega )\text { is invertible with bounded inverse} \end{aligned}$$

(2.45)

(${\mathcal {D}}(\partial {{\widetilde{E}}})$ is identified with a subspace of $W^{2,2}_{n}(\Omega )$ as in (2.6)), where $I:{\mathcal {D}}(\partial {\widetilde{E}})\longrightarrow L^{2}(\Omega )$ is the inclusion map. Moreover, we can follow arguments from [6, pp. 466–467] to show that the inverse operator

$$\begin{aligned} (\partial {{\widetilde{E}}}+I)^{-1}:L^{2}(\Omega )\longrightarrow W^{2-s,2}_n(\Omega )\text { is continuous for any }s\in (0,\tfrac{1}{4}). \end{aligned}$$

(2.46)

From now on we fix $s=\frac{1}{8}$ for definiteness.

Finally let us consider the last component of ${\mathcal {N}}_{k}.$ From standard elliptic theory, the operator

$$\begin{aligned} -\Delta (\cdot )+\int _{\Omega }\cdot : W^{2,2}_{n}(\Omega )\longrightarrow L^{2}(\Omega ) \end{aligned}$$

is bounded invertible with bounded and continuous inverse.

In summary we have shown that

$$\begin{aligned} \text {the map }{\mathcal {N}}_{k}:X\longrightarrow Y\text { is bounded, invertible, the inverse is bounded}\nonumber \\ \end{aligned}$$

(2.47)

and the inverse map is continuous from Y to ${\widetilde{X}},$ where

$$\begin{aligned} \begin{array}{l} {\widetilde{X}}=W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )\times W^{2,2}_{n}(\Omega )\times W^{\frac{15}{8},2}_{n}(\Omega )\times W^{2,2}_{n}(\Omega ). \end{array} \end{aligned}$$

(2.48)

Next we show that ${\mathcal {F}}_{k}:{\widetilde{X}}\longrightarrow Y$ is continuous and compact.

2.4.2 The operator ${\mathcal {F}}_{k}:{\widetilde{X}}\longrightarrow Y$ is continuous and compact

For the operator ${\mathcal {F}}_{k}=({\mathcal {F}}_{k}^{1},{\mathcal {F}}_{k}^{2},{\mathcal {F}}_{k}^{3},{\mathcal {F}}_{k}^{4}),$ we will verify the continuity and compactness component wise.

Let us start with ${\mathcal {F}}^{1}_k:{\widetilde{X}}\longrightarrow (W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))'.$ In this direction we will collect estimates of terms appearing in the expression of ${\mathcal {F}}_{k}^{1}.$ The following estimates are obtained by using Hölder’s inequality and standard Sobolev embeddings.

$$\begin{aligned} \begin{aligned} \Vert \rho v\Vert _{L^{\frac{3}{2}}(\Omega )}&\leqslant C\Vert v\Vert _{W^{1,2}(\Omega )}\left( \Vert \phi \Vert _{L^{2}(\Omega )}+1\right) ,\\ \Vert \text{ div }(\rho _{k}v\otimes v)\Vert _{L^{\frac{3}{2}}(\Omega )}&\leqslant C_{k}\Vert v\Vert ^{2}_{W^{1,2}(\Omega )},\\ \Vert \nabla \mu \phi _{k}\Vert _{L^{\frac{3}{2}}(\Omega )}&\leqslant C_{k}\Vert \mu \Vert _{W^{2,2}(\Omega )},\\ \Vert \xi (\phi _{k})(|M|^{2}M-M_{k})\nabla M\Vert _{L^{\frac{3}{2}}(\Omega )}&\leqslant C_{k}\left( \Vert M\Vert ^{4}_{W^{2,2}(\Omega )}+\Vert M\Vert _{W^{2,2}(\Omega )}\right) ,\\ \Vert \textrm{div}(\xi (\phi _{k})\nabla M)\nabla M\Vert _{L^{\frac{3}{2}}(\Omega )}&\leqslant C_{k}\Vert M\Vert ^{2}_{W^{2,2}(\Omega )},\\ \Vert \left( \text{ div }J\right) v\Vert _{L^{\frac{3}{2}}(\Omega )}&\leqslant C\Vert v\Vert _{W^{1,2}(\Omega )}\Vert \mu \Vert _{W^{2,2}(\Omega )},\\ \Vert \left( J\cdot \nabla \right) v\Vert _{L^{\frac{3}{2}}(\Omega )}&\leqslant C\Vert v\Vert _{W^{1,2}(\Omega )}\Vert \mu \Vert _{W^{2,2}(\Omega )},\\ \Vert (v\cdot \nabla \rho _{k})v\Vert _{L^{\frac{3}{2}}(\Omega )}&\leqslant C_{k}\Vert v\Vert ^{2}_{W^{1,2}(\Omega )}. \end{aligned} \end{aligned}$$

(2.49)

Hence the estimates above prove the boundedness of ${\mathcal {F}}_{k}^{1}$ from ${\widetilde{X}}$ to $L^{\frac{3}{2}}(\Omega ).$ One can use similar estimates to show that ${\mathcal {F}}_{k}^{1}$ is continuous from ${\widetilde{X}}$ to $L^{\frac{3}{2}}(\Omega ).$ Next the compact embedding $L^{\frac{3}{2}}(\Omega )\longrightarrow (W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))'$ furnishes the continuity and compactness of ${\mathcal {F}}^{1}_{k}:{\widetilde{X}}\longrightarrow (W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))'.$

Now we prove that ${\mathcal {F}}_{k}^{2}:{\widetilde{X}}\longrightarrow L^{2}(\Omega )$ is continuous and compact. One first observes

$$\begin{aligned} \displaystyle \Vert (v\cdot \nabla )M\Vert _{W^{1,\frac{3}{2}}(\Omega )}\leqslant C\Vert v\Vert _{W^{1,2}(\Omega )}\Vert M\Vert _{W^{2,2}(\Omega )} \end{aligned}$$

(2.50)

as in [36, p. 16]. The continuity of ${\mathcal {F}}_{k}^{2}:{\widetilde{X}}\longrightarrow W^{1,\frac{3}{2}}(\Omega )$ relies on a similar estimate and can be concluded without any difficulty. Further in view of the compactness of the embedding of $W^{1,\frac{3}{2}}(\Omega )$ to $L^2(\Omega )$ the asserted continuity and compactness follows.

Next we show that ${\mathcal {F}}_{k}^{3}: {\widetilde{X}}\longrightarrow L^{2}(\Omega )$ is continuous and compact. For the proof we take a different route than the ones used for ${\mathcal {F}}^{i}_{k},$ $i\in \{1,2\}.$ First we introduce

$$\begin{aligned} {\widetilde{Y}}=L^{2}(\Omega )\times W^{\frac{15}{8},2}(\Omega )\times W^{\frac{7}{4},2}(\Omega )\times W^{\frac{15}{8},2}(\Omega ). \end{aligned}$$

One observes that the embedding ${\widetilde{X}}\longrightarrow {\widetilde{Y}}$ is compact (since in dimension three the embedding $W^{m+k,2}\hookrightarrow W^{m,2}$ is compact for any $0<k<\frac{3}{2}$). We study the boundedness and continuity of the operator ${\mathcal {G}}:{{\widetilde{Y}}}\longrightarrow L^2(\Omega )$ that is defined as in the third component of (2.41). Once we verify its boundedness and continuity we immediately conclude that $\mathcal F^3_k:{\widetilde{X}}\longrightarrow L^{2}(\Omega )$ is bounded, continuous and compact as the composition of the embedding ${\widetilde{X}}\longrightarrow {\widetilde{Y}}$ and ${\mathcal {G}}$. In order to conclude the boundedness of ${\mathcal {G}}$, we collect the following estimates:

$$\begin{aligned} \begin{aligned} \Vert H_{0}(\phi ,\phi _{k})|\nabla M|^{2}\Vert _{L^{2}(\Omega )}&\leqslant C_{k}\Vert M\Vert ^{2}_{W^{\frac{15}{8},2}(\Omega )},\\ \Vert H_{0}(\phi ,\phi _{k})(|M|^{2}-1)^{2}\Vert _{L^{2}(\Omega )}&\leqslant C_{k}\left( \Vert M\Vert ^{4}_{W^{\frac{15}{8},2}(\Omega )}+1\right) . \end{aligned} \end{aligned}$$

(2.51)

Indeed, since $\xi (\cdot )\in C^{1}(\mathbb {R})$, one infers $\Vert H_{0}(\phi ,\phi _{k})\Vert _{L^{\infty }(\Omega )}\leqslant C_{k}$ by using the mean value theorem and the upper bound of $\xi '(\cdot )$ (cf. assumption (1.6)). Further the fact that $\nabla M$ is bounded in $W^{\frac{7}{8},2}(\Omega )$ and the continuous embedding $W^{\frac{7}{8},2}(\Omega )\hookrightarrow L^{4}(\Omega )$ proves (2.51)$_{1}.$ Whereas (2.51)$_{2}$ is a consequence of the continuous embedding $W^{\frac{15}{8},2}(\Omega )\hookrightarrow L^{\infty }(\Omega ).$ The boundedness and continuity of the linear terms in the expression of ${\mathcal {F}}_{k}^{3}$ are trivially concluded.

In the spirit of the boundedness estimates (2.51), the continuity of the nonlinear terms in ${\mathcal {F}}_{k}^{3},$ i.e. $\displaystyle \frac{1}{2}H_{0}(\phi ,\phi _{k})|\nabla M|^{2}$ and $\displaystyle \frac{1}{4\alpha ^2}H_{0}(\phi ,\phi _{k})(|M|^{2}-1)^{2},$ can be proved by the arguments as in [36, p. 16] adjusted to the current set-up. Hence we have proved that ${\mathcal {F}}_{k}^{3}: {\widetilde{X}}\longrightarrow L^{2}(\Omega )$ is continuous and compact.

Next we consider ${\mathcal {F}}_{k}^{4}.$ In the spirit of (2.50) we derive the following estimate

$$\begin{aligned} \displaystyle \Vert (v\cdot \nabla )\phi _{k}\Vert _{W^{1,\frac{3}{2}}(\Omega )} \leqslant C_{k}\Vert v\Vert _{W^{,2}(\Omega )}, \end{aligned}$$

(2.52)

which verifies the boundedness of ${\mathcal {F}}_{k}^{4}:{\widetilde{X}}\longrightarrow W^{1,\frac{3}{2}}(\Omega ).$ By (2.52) the continuity of ${\mathcal {F}}_{k}^{4}:{\widetilde{X}}\longrightarrow W^{1,\frac{3}{2}}(\Omega )$ can be concluded in a straight forward manner. Next using the compact embedding $W^{1,\frac{3}{2}}(\Omega )\hookrightarrow L^{2}(\Omega ),$ one at once renders that ${\mathcal {F}}_{k}^{4}:{\widetilde{X}}\longrightarrow L^{2}(\Omega )$ is continuous and compact.

In view of the arguments above we finally have proved that ${\mathcal {F}}_{k}:{\widetilde{X}}\longrightarrow Y$ is continuous and compact.

2.4.3 The fixed point argument

We now show the existence of a $w\in X$ (we recall the definition of X from (2.39)$_{1}$) satisfying

$$\begin{aligned} \begin{array}{l} {\mathcal {N}}_{k}(w)={\mathcal {F}}_{k}(w)\,\,\text{ in }\,\, Y. \end{array} \end{aligned}$$

(2.53)

For that purpose it is sufficient to prove the existence of a fixed point of the operator ${\mathcal {F}}_{k}\circ {\mathcal {N}}^{-1}_{k}$ on Y, i.e., the existence of $z\in Y$ satisfying

$$\begin{aligned} \begin{array}{l} z=({\mathcal {F}}_{k}\circ {\mathcal {N}}^{-1}_{k})z \,\,\text{ in }\,\,Y, \end{array} \end{aligned}$$

(2.54)

since the invertibility of the operator ${\mathcal {N}}_{k}:{X}\rightarrow Y$ implies the obtainment of $w\in {X}$ satisfying (2.53) by using $w={\mathcal {N}}^{-1}_{k}(z)$.

In order to show the existence of a fixed point of the operator equation (2.54) we apply the Leray-Schauder fixed point theorem [28, Theorem 10.3] to the continuous and compact operator ${\mathcal {F}}_{k}\circ {\mathcal {N}}^{-1}_{k}:Y\longrightarrow Y.$ To this end we verify that:

$$\begin{aligned} \begin{aligned}&\text {There exists }r\!>\!0 \text { such that if } z\in Y \text { solves } z\!=\!\lambda ({\mathcal {F}}_{k}\circ {\mathcal {N}}^{-1}_{k})z\text { with }\lambda \!\in \![0,1],\\&\text{ then } \text{ it } \text{ holds }\,\,\Vert z\Vert _{Y}\leqslant r. \end{aligned}\nonumber \\ \end{aligned}$$

(2.55)

Let $z\in Y$ satisfy $z=\lambda ({\mathcal {F}}_{k}\circ {\mathcal {N}}^{-1}_{k})z$ in Y with some $\lambda \in [0,1]$. Then

$$\begin{aligned} w=(v,M,\phi ,\mu )={\mathcal {N}}_{k}^{-1}z, \end{aligned}$$

solves

$$\begin{aligned} \begin{array}{l} {\mathcal {N}}_{k}(w)-\lambda {\mathcal {F}}_{k}(w)=0\,\,\text{ in }\,\, Y. \end{array} \end{aligned}$$

(2.56)

Let us first prove that

$$\begin{aligned} \Vert w\Vert _{{\widetilde{X}}}\leqslant C_{k}, \end{aligned}$$

(2.57)

with $C_k$ independent of $\lambda \in [0,1].$ Then we will bootstrap the regularity to have

$$\begin{aligned} \Vert w\Vert _{{X}}\leqslant C_{k}, \end{aligned}$$

(2.58)

with $C_k$ independent of $\lambda \in [0,1],$ from which (2.55) follows due to the boundedness of ${\mathcal {N}}_k: X\longrightarrow Y$, cf. (2.47). One recalls the definitions of ${\mathcal {N}}_{k}$ and ${\mathcal {F}}_{k}$ from (2.40) and (2.41), tests the first component of (2.56) by v, the second component by $-{{\,\textrm{div}\,}}(\xi (\phi _k)\nabla M)+\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})$, the third component by $\frac{\phi -\phi _{k}}{h}$ and the fourth component by $\mu$. The application of (2.36), Lemma 2.5 and identities (2.30) (similar arguments leading to (2.37) from (2.35)) yield the following after dropping some positive terms from the left hand side (similar to the obtainment of (2.29) from (2.37))

$$\begin{aligned} \begin{aligned}&\frac{\lambda }{h}\left( \frac{1}{2}\int _{\Omega }\rho |v|^{2}-\frac{1}{2}\int _{\Omega }\rho _{k}|v_{k}|^{2}\right) +2\int _{\Omega }\nu (\phi _{k})|\mathbb {D}v|^{2}\\&\quad +\frac{\lambda }{h}\left( \frac{1}{2}\int _{\Omega }\xi (\phi )|\nabla M|^{2}-\frac{1}{2}\int _{\Omega }\xi (\phi _{k})|\nabla M_{k}|^{2}\right) \\&\quad +\frac{\lambda }{h}\left( \frac{1}{4\alpha ^{2}}\int _{\Omega }\xi (\phi )\bigl (|M|^{2}-1\bigr )^{2}-\frac{1}{4\alpha ^{2}}\int _{\Omega }\xi (\phi _{k})\bigl (|M_{k}|^{2}-1\bigr )^{2}\right) \\&\quad +\int _{\Omega }\left| \text{ div }(\xi (\phi _{k})\nabla M)-\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})\right| ^{2} \\&\quad +\frac{1}{h}\left( \frac{\eta }{2}\int _{\Omega }|\nabla \phi |^{2}-\frac{\eta }{2}\int _{\Omega }|\nabla \phi _{k}|^{2}\right) +\frac{1}{h}\int _{\Omega }\left( {\widetilde{\Psi }}_{0}(\phi )-{\widetilde{\Psi }}_{0}(\phi _{k})\right) \\&\quad -\frac{\lambda }{h}\int _{\Omega }\kappa \frac{(\phi ^{2}-\phi _{k}^{2})}{2}\!+\!\int _{\Omega }|\nabla \mu |^{2}\!+\!(1-\lambda )\left( \int _\Omega \mu \right) ^{2}\!+\!\frac{(1-\lambda )}{h}\int _{\Omega }\left( \frac{\phi ^{2}}{2}-\frac{\phi ^{2}_{k}}{2}\right) \\&\quad +(1-\lambda )\int _\Omega M \int _\Omega \left( -\text{ div }(\xi (\phi _{k})\nabla M)+\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})\right) \leqslant 0. \end{aligned} \end{aligned}$$

(2.59)

Once again we recall that in obtaining the above inequality one expands $\xi (\phi _{k})\nabla M\cdot (\nabla M-\nabla M_{k})$ and $\nabla \phi \cdot (\nabla \phi -\nabla \phi _{k})$ by using (2.36) and uses Lemma 2.5 to expand the term $(M-M_{k})\cdot \frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k}).$ In order to obtain the inequality (2.59) we also have used (2.38).

When $0\leqslant \lambda <1,$ we use Young’s and Hölder’s inequality to estimate the term appearing in the last line of (2.59) to infer:

$$\begin{aligned} \begin{aligned}&\left| (1-\lambda )\int _\Omega M \int _{\Omega }\left( -\text{ div }(\xi (\phi _{k})\nabla M)+\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})\right) \right| \\&\quad \leqslant (1\!-\!\lambda )\epsilon |\Omega |\int _{\Omega }\left| \text{ div }(\xi (\phi _{k})\nabla M)\!-\!\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M\!-\!M_{k})\right| ^{2}\!+\!c\frac{1-\lambda }{\epsilon }\left( \int _\Omega M\right) ^{2} \end{aligned}\nonumber \\ \end{aligned}$$

(2.60)

for some positive parameter $\epsilon >0$.

Since $|1-\lambda |\leqslant 1,$ for small enough choice of the parameter $\epsilon >0,$ the first term on the right hand side (2.60) can be absorbed in the fifth summand appearing on the left hand side of (2.59). We will now estimate the second term on the right hand side of (2.60). We will now estimate the second term on the right hand side of (2.60). In that direction we follow the arguments used to show [36, p. 18, (4.32)] to obtain

$$\begin{aligned} \ \int _{\Omega }|\nabla M|^{2}+\int _{\Omega } |M|^{4}+(1-\lambda )\left( \int _\Omega M\right) ^{2}\leqslant C_{k}, \end{aligned}$$

(2.61)

where $C_{k}>0$ is independent of $\lambda >0.$

For a small enough choice of the parameter $\epsilon >0,$ using (2.61) and (2.60) in (2.59) we obtain in particular

$$\begin{aligned} \begin{aligned}&2h\int _{\Omega }\nu (\phi _{k})|\mathbb {D}v|^{2} +\frac{h}{2}\int _{\Omega }\left| \text{ div }(\xi (\phi _{k})\nabla M)-\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})\right| ^{2}\\&\qquad +\frac{\eta }{2}\int _{\Omega }|\nabla \phi |^{2}+\int _{\Omega }{\widetilde{\Psi }}_{0}(\phi )+h\int _{\Omega }|\nabla \mu |^{2}+h(1-\lambda )\left( \int _\Omega \mu \right) ^{2}-\lambda \int _{\Omega }\kappa \frac{\phi ^{2}}{2}\\&\quad \leqslant \int _{\Omega }\frac{\rho _{k}|v_{k}|^{2}}{2}+\frac{1}{2}\int _{\Omega }\phi _{k}^{2}+\frac{\eta }{2}\int _{\Omega }|\nabla \phi _{k}|^{2}+\int _{\Omega }|\kappa |\frac{\phi ^{2}_{k}}{2}+\int _{\Omega }{\widetilde{\Psi }}_{0}(\phi _{k})\\&\qquad +\frac{1}{2}\int _{\Omega }\xi (\phi _{k})|\nabla M_{k}|^{2}+\frac{1}{4\alpha ^{2}}\int _{\Omega }\xi (\phi _{k})\bigl (|M_{k}|^{2}-1\bigr )^{2}+ C_{k} \leqslant C_{k}. \end{aligned}\nonumber \\ \end{aligned}$$

(2.62)

Indeed, we obtain (2.62) from (2.59) by using the positive lower bound on $\rho$ from (1.10), which follows since $w=(v,M,\phi ,\mu )={\mathcal {N}}_{k}^{-1}z\in X$, cf. (2.47), implying $\phi \in {\mathcal {D}}(\partial {\widetilde{E}})$ and hence $\phi \in [-1,1]$ a.e. The fact that $\phi \in [-1,1]$ a.e. and $\lambda \in [0,1]$ alongside ${\widetilde{\Psi }}_0\in C([-1,1])$ allows us to obtain

$$\begin{aligned} \left| \int _{\Omega }{\widetilde{\Psi }}_0(\phi )\right| \leqslant C\quad \text{ and }\quad \left| \int _{\Omega }\kappa \frac{\phi ^{2}}{2}\right| \leqslant C, \end{aligned}$$

for some positive constant C and hence the term $\displaystyle \int _{\Omega }{\widetilde{\Psi }}_{0}(\phi )$ and $\displaystyle -\lambda \int _{\Omega }\kappa \frac{\phi ^{2}}{2}$ can be dropped from the left hand side of (2.62). Hence

$$\begin{aligned} \begin{aligned}&2h\int _{\Omega }\nu (\phi _{k})|\mathbb {D}v|^{2} +\frac{h}{2}\int _{\Omega }\left| \text{ div }(\xi (\phi _{k})\nabla M)-\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})\right| ^{2}\\&\quad +\frac{\eta }{2}\int _{\Omega }|\nabla \phi |^{2}+h\int _{\Omega }|\nabla \mu |^{2}+h(1-\lambda )\left( \int _\Omega \mu \right) ^{2}\leqslant C_{k}. \end{aligned} \end{aligned}$$

(2.63)

One uses (2.63), (1.7), Korn’s and Poincaré’s inequality and the fact that $\phi \in [-1,1]$ a.e. to render $\Vert v\Vert _{W^{1,2}(\Omega )}+\Vert \phi \Vert _{W^{1,2}(\Omega )}\leqslant C_{k}.$ We conclude $\Vert M\Vert _{W^{1,2}(\Omega )}\leqslant C_{k}$ independently of $\lambda$ from (2.61). The second term of (2.63)$_{1}$ provides

$$\begin{aligned} \left\| {{\,\textrm{div}\,}}(\xi (\phi _{k})\nabla M)-\frac{\xi (\phi _{k})}{\alpha ^{2}}(|M|^{2}M-M_{k})\right\| _{L^{2}(\Omega )}\leqslant C_{k}. \end{aligned}$$

(2.64)

One can now use elliptic regularity results and a bootstrap argument (exactly as the one used to prove the claim (2.44), cf. [36, Section 4.2.1, pp. 14–15]) to furnish that $\Vert M\Vert _{W^{2,2}(\Omega )}\leqslant C_k$ from (2.64).

Next one obtains $\Vert \nabla \mu \Vert _{L^{2}(\Omega )}\leqslant C_{k}$ from (2.63). In view of Poincaré’s inequality it is sufficient to show

$$\begin{aligned} \left| \int _\Omega \mu \right| \leqslant C_k, \end{aligned}$$

(2.65)

in order to prove $\Vert \mu \Vert _{W^{1,2}(\Omega )}\leqslant C_{k}.$ For $\lambda \in [0,\frac{1}{2}),$ (2.65) follows from the estimate of the last term which appears in the left hand side of (2.63). For $\lambda \in [\frac{1}{2},1]$ we follow similar arguments used to show (2.21). Hence (2.65) holds independently of the values of $\lambda \in [0,1]$ and as a consequence $\Vert \mu \Vert _{W^{1,2}(\Omega )}\leqslant C_{k}$ follows.

So far we have obtained the following

$$\begin{aligned} \begin{array}{l} \displaystyle \Vert v\Vert _{W^{1,2}(\Omega )}+\Vert M\Vert _{W^{2,2}(\Omega )}+\Vert \phi \Vert _{W^{1,2}(\Omega )}+\Vert \mu \Vert _{W^{1,2}(\Omega )}\leqslant C_{k}. \end{array} \end{aligned}$$

(2.66)

Further one recalls that $\mu$ solves the following equation

$$\begin{aligned} \begin{aligned} \displaystyle -\Delta \mu +\int _\Omega {\mu }&=-\lambda \frac{\phi -\phi _{k}}{h}-\lambda (v\cdot \nabla )\phi _{k}+\lambda \int _\Omega \mu{} & {} \quad \text{ in }\quad \Omega ,\\ \partial _{n}\mu&=0{} & {} \quad \text{ in }\quad \partial \Omega . \end{aligned} \end{aligned}$$

(2.67)

In view of the fact that $\lambda \in [0,1]$ and the estimate (2.66) we observe that the right hand side of (2.67)$_1$ can be estimated in $L^{2}(\Omega )$ and hence by standard elliptic regularity theory

$$\begin{aligned} \Vert \mu \Vert _{W^{2,2}(\Omega )}\leqslant C_{k}. \end{aligned}$$

Further, from the identity

$$\begin{aligned} \begin{aligned} \displaystyle \partial {\widetilde{E}}(\phi )+\phi =&\lambda \phi +\lambda \mu +\lambda \kappa \frac{\phi +\phi _{k}}{2}-\lambda H_{0}(\phi ,\phi _{k})\frac{|\nabla M|^{2}}{2}\\&-\lambda \frac{H_{0}(\phi ,\phi _{k})}{4\alpha ^{2}}\left( |M|^{2}-1\right) ^{2} \end{aligned} \end{aligned}$$

(2.68)

and (2.66) one has

$$\begin{aligned} \Vert \partial {\widetilde{E}}(\phi )+\phi \Vert _{L^{2}(\Omega )}\leqslant C_{k}. \end{aligned}$$

(2.69)

Inequality (2.69) along with the estimate of $\phi$ from (2.66) imply that $\Vert \partial {\widetilde{E}}(\phi )\Vert _{L^{2}(\Omega )}\leqslant C_{k}.$ Next in view of the inequality (2.7) one in particular concludes that

$$\begin{aligned} \Vert \phi \Vert _{W^{\frac{15}{8},2}(\Omega )}\leqslant C_{k} \end{aligned}$$

and hence altogether we have

$$\begin{aligned} \begin{array}{l} \Vert w\Vert _{{\widetilde{X}}}+\Vert \partial {\widetilde{E}}(\phi )\Vert _{L^{2}(\Omega )}=\Vert (v,M,\phi ,\mu )\Vert _{\widetilde{X}}+\Vert \partial {\widetilde{E}}(\phi )\Vert _{L^{2}(\Omega )}\leqslant C_{k}. \end{array} \end{aligned}$$

(2.70)

Once again using (2.7) and (2.70) one at once concludes (2.58) and consequently proves (2.55). Finally the fact that ${\mathcal {N}}_{k}:X\rightarrow Y$ has a bounded inverse yields the existence of a fixed point to the operator equation (2.53). This finishes the proof of Theorem 2.6. $\square$

3 Proof of Theorem 1.5

Let $T>0$ be fixed. Let $0=t_0<t_1<\ldots<t_k<\ldots$, $k\in \mathbb {N}_0$ be a strictly increasing sequence such that for each $k\in \mathbb {N}_0$ $t_{k+1}-t_k=h$ where $h=\frac{1}{N}$ for $N\in \mathbb {N}$ fixed. Applying Theorem 2.6 successively, we construct a sequence $\{(v_{k+1},M_{k+1},\phi _{k+1},\mu _{k+1})\}$, $k\in \mathbb {N}_0$ of solutions to problem (2.12) assuming $(v_k,M_k,\phi _k)\in L^2_{{{\,\textrm{div}\,}}}(\Omega )\times W^{2,2}_n(\Omega )\times W^{2,2}_n(\Omega )$ with $-1\leqslant \phi _{k}\leqslant 1.$ Obviously, the assumption $(M_0,\phi _0)\in W^{1,2}(\Omega )\times W^{1,2}(\Omega )$ excludes the possibility of application of Theorem 2.6 directly. Instead, we consider sequences $\{M_0^N\}\subset W^{2,2}_n(\Omega )$, $\{\phi ^N_0\}\subset W^{2,2}_n(\Omega )$ with $|\phi ^{N}_{0}|\leqslant 1$ such that

$$\begin{aligned} \begin{aligned} M^N_0&\rightarrow M_0{} & {} \text { in }W^{1,2}(\Omega ),\\ \phi ^N_0&\rightarrow \phi _0{} & {} \text { in }W^{1,2}(\Omega ) \end{aligned} \end{aligned}$$

(3.1)

as $N\rightarrow \infty$. Such an approximating sequence $\{\phi ^{N}_{0}\}$ can be constructed by solving a heat equation with initial data $\phi _{0},$ setting $\phi ^{N}_{0}$ as the solution to the heat equation at $t=\frac{1}{N}$ and using parabolic regularity. The details of this construction can be found in [6, Section 5.1, p. 471]. Similar arguments apply in constructing $\{M^{N}_{0}\}.$ Adapting the notation from [36] we introduce two types of interpolants related to the unknowns. At first we fix $N\in \mathbb {N}$. The piecewise constant interpolants of $(v,M,\phi )$ are defined on $[-h,\infty )$ and the one of $\mu$ on $[0,\infty )$ via

$$\begin{aligned} v^N(t)=v_0,\ M^N(t)=M^N_0,\ \phi ^N(t)=\phi ^N_0\text { for }t\in [-h,0) \end{aligned}$$

(3.2)

and

$$\begin{aligned} f^N(t)=f_k\text { for }t\in [(k-1)h,kh), \end{aligned}$$

(3.3)

where $f^N$ stands for the interpolants $v^N$, $M^N$, $\phi ^N$, $\mu ^N$ and $f_k$ represents the corresponding $v_k$, $M_k$, $\phi _k$ and $\mu _k$, $k\in \mathbb {N}$. We note that

$$\begin{aligned} \begin{array}{ll} \rho ^N=\frac{1}{2}({{\widetilde{\rho }}}_1+{{\widetilde{\rho }}}_2)+\frac{1}{2}({{\widetilde{\rho }}}_2-{{\widetilde{\rho }}}_1)\phi ^N. \end{array} \end{aligned}$$

(3.4)

Next, a piecewise affine interpolant ${{\widetilde{f}}}^N$ is defined as

$$\begin{aligned} \widetilde{f}^N(t)=\frac{(k+1)h-t}{h}f^N(t-h)+\frac{t-kh}{h}f^N(t)\text { for }t\in [kh,(k+1)h),\,\,k\in \mathbb {N}_{0}.\nonumber \\ \end{aligned}$$

(3.5)

For our purposes it is sufficient to consider only the piecewise affine interpolants $\{\widetilde{\rho v}^N\}$, $\{\widetilde{M}^N\}$ and $\{{{\widetilde{\phi }}}^N\}$, where the convention $(\rho v)^{N}=\rho ^{N}v^{N}$ is used. We also introduce the notation for the shift and the difference quotient in time of a function f as follows

$$\begin{aligned} f_h(t)=&f(t-h),\\ \partial ^{-}_{t,h}f(t)=&\frac{1}{h}(f-f_h)(t). \end{aligned}$$

As immediate consequences of the latter definition and (3.5) one gets

$$\begin{aligned} \begin{aligned} \partial _t{{\widetilde{f}}}^N(t)&=\partial ^-_{t,h}f^N(t)\,\,\text{ for } \text{ all }\,t\in [kh,(k+1)h),\,k\in \mathbb {N}_{0},\\ \Vert {{\widetilde{f}}}^N\Vert _{L^p(0,\tau ;X)}&\leqslant \Vert f^N\Vert _{L^p(0,\tau ;X)}\!+\!\Vert f^N_h\Vert _{L^p(0,\tau ;X)}\text { for any }p\in [1,\infty ]\,\,\text{ and }\,\,\tau >0. \end{aligned} \end{aligned}$$

(3.6)

We will next state identities that are satisfied by interpolants. Let $\tau \in (0,\infty )$ be chosen arbitrarily. We find $k_\tau \in \mathbb {N}_0$ such that $\tau \in [k_\tau h,(k_{\tau }+1)h)$. Further, we fix an arbitrary $\psi _1\in L^2(0,\infty ;V(\Omega ))$, set ${{\widetilde{\psi }}}_1=\int _{kh}^b\psi _1$ in (2.17), where

$$\begin{aligned} b= {\left\{ \begin{array}{ll} (k+1)h &{} k<k_\tau ,\\ \tau &{} k=k_\tau , \end{array}\right. } \end{aligned}$$

and sum the resulting expressions over $k\in \{0,1,\dots , k_\tau \}$ to obtain

$$\begin{aligned} \begin{aligned}&\int _0^\tau \left( \int _{\Omega } \partial ^-_{t,h}(\rho ^Nv^N)\cdot \psi _{1}-\int _{\Omega }(\rho ^N_hv^N\otimes v^N)\cdot \nabla {\psi }_{1}-\int _{\Omega }(v^N\otimes J^N)\cdot \nabla \psi _{1}\right. \\&\qquad \left. -\int _{\Omega }\left( \frac{{\xi (\phi ^N_h)}}{\alpha ^{2}}(|M^N|^{2}M^N-M^N_h)\nabla M^N\right) \cdot \psi _{1}\right. \\&\qquad \left. +\int _{\Omega }\left( {{\,\textrm{div}\,}}\left( \xi (\phi ^N_h)\nabla M^N\right) \nabla M^N\right) \cdot \psi _{1}\right) \\&\quad =\int _0^\tau \left( -2\int _{\Omega }\nu (\phi ^{N}_{h})\mathbb {D}v^{N}\cdot \mathbb {D}\psi _{1}-\int _{\Omega }\nabla \mu ^N\phi ^N_h\cdot \psi _{1}\right) \end{aligned} \end{aligned}$$

(3.7)

for all $\tau \in (0,\infty )$ and $\psi _1\in L^2(0,\infty ;V(\Omega )),$ where $\displaystyle J^{N}=-\frac{{\widetilde{\rho }}_{2}-{\widetilde{\rho }}_{1}}{2}\nabla \mu ^{N}$. In obtaining (3.7) from (2.17), we once again use identity (2.15). Similarly, we get

$$\begin{aligned} \begin{aligned}&\int _{0}^{\tau }\left( \int _{\Omega }\partial ^{-}_{t,h}M^{N}\cdot \psi _{2}+\int _{\Omega }(v^{N}\cdot \nabla )M^{N}\cdot \psi _{2}\right) \\&\quad =\int _{0}^{\tau }\int _{\Omega }\left( {{\,\textrm{div}\,}}({\xi (\phi ^{N}_{h})}\nabla M^{N})-\frac{\xi (\phi ^{N}_{h})}{\alpha ^{2}}(|M^{N}|^{2}M^{N}-M^{N}_{h})\right) \cdot \psi _{2}, \end{aligned} \end{aligned}$$

(3.8)

for all $\tau \in (0,\infty )$ and $\psi _{2}\in L^2(0,\infty ;W^{1,2}(\Omega ))$. Moreover, by obvious manipulations we deduce from (2.19) that

$$\begin{aligned} \int _{0}^{\tau }\left( \int _{\Omega } \partial ^{-}_{t,h}\phi ^{N}\psi _{3}+\int _{\Omega }(v^{N}\cdot \nabla )\phi ^{N}_{h}\psi _{3}\right) =-\int _{0}^{\tau }\int _{\Omega }\nabla \mu ^{N}\cdot \nabla \psi _{3} \end{aligned}$$

(3.9)

for all $\tau \in (0,\infty )$ and $\psi _3\in L^2(0,\infty ;W^{1,2}(\Omega ))$ and from (2.20) it follows that

$$\begin{aligned} \begin{aligned}&\int _0^\tau \int _\Omega \left( \mu ^N+\kappa \frac{\phi ^N+\phi ^N_h}{2}-H_{0}(\phi ^N,\phi ^N_h)\frac{|\nabla M^N|^{2}}{2}\right. \\&\quad \left. -\frac{H_{0}(\phi ^N,\phi ^N_h)}{4\alpha ^{2}}(|M^N|^{2}-1)^{2}\right) \psi _4=\int _0^\tau \int _\Omega \left( -\eta \Delta \phi ^N+{\widetilde{\Psi }}'_{0}(\phi ^N)\right) \psi _4 \end{aligned} \end{aligned}$$

(3.10)

for all $\tau \in (0,\infty )$ and $\psi _4\in L^\infty (0,\tau ;L^\infty (\Omega ))$.

3.1 Compactness of sequences of interpolants

The goal of this section is to collect all the necessary convergences of (sub)sequences of interpolants allowing for the passage $h\rightarrow 0$ (equivalently $N\rightarrow \infty$) in order to show the existence of a weak solution to the original problem.

3.1.1 Uniform bounds on sequences of interpolants and compactness

In order to obtain the uniform bounds we begin with the energy inequality for the interpolants $v^N$, $M^N$, $\phi ^N$ and $\mu ^N$. Summing in (2.29) over $k\in \mathbb {N}_0$ we obtain

$$\begin{aligned} \begin{aligned}&E_{tot}(v^N(t),M^N(t),\phi ^N(t))+2\int _0^t\int _\Omega \nu (\phi ^{N}_{h})|\mathbb {D} v^{N}|^{2}+\int _0^t\int _\Omega |\nabla \mu ^N|^2\\&\qquad +\int _0^t\int _\Omega \left| {{\,\textrm{div}\,}}(\xi (\phi ^N_h)\nabla M^N)-\frac{\xi (\phi ^N_h)}{\alpha ^2}(|M^N|^2M^N-M^N_h)\right| ^2\\&\quad \leqslant E_{tot}(v_0,M^N_0,\phi ^N_0)\leqslant C\left( E_{tot}(v_{0},M_{0},\phi _{0})+1\right) \end{aligned} \end{aligned}$$

(3.11)

first for each $t\in h\mathbb {N}_0,$ where the inequality in the final line of (3.11) follows from (3.1). As all interpolants involved in the latter inequality are constant on intervals of the form $[kh,(k+1)h),$ one concludes that (3.11) holds for all $t\in [0,\infty )$. Therefore recalling the definition of $E_{tot}$ in (1.17) and using (1.10) we conclude from (3.11) that in particular uniform bounds on

$$\begin{aligned} \begin{aligned}&\{v^N\}\text { in } L^\infty (0,T+1;L^2(\Omega ))\cap L^2(0,T+1;W^{1,2}(\Omega )),\\&\{v^N\}\text { in }L^\frac{10}{3}(Q_{T+1}),\\&\{M^N\}\text { in } L^\infty (0,T+1;W^{1,2}(\Omega )),\\&\{\phi ^N\}\text { in } L^\infty (0,T+1;W^{1,2}(\Omega )),\\&\{\nabla \mu ^N\}\text { in } L^2(Q_{T+1}),\\&\{{{\,\textrm{div}\,}}(\xi (\phi ^N_h)\nabla M^N)\!-\!\frac{\xi (\phi ^N_h)}{\alpha ^2}(|M^N|^2M^N\!-\!M^N_h)\}\text { in }L^2(0,T\!+\!1;L^2(\Omega )),\\&|\phi ^{N}|\leqslant 1\,\,\text {a.e. in}\,\, Q_{T}. \end{aligned} \end{aligned}$$

(3.12)

The bound (3.12)$_{1}$ is obtained by using (1.7) and Korn inequality. All the sequences in (3.12)$_{1}$-(3.12)$_{6}$ in the respective norms are bounded by $C(E_{tot}(v_{0}$, $M_{0},\phi _{0})+1)^{\frac{1}{2}}$ as a consequence of (3.11). Since $\phi _{k+1}\in {\mathcal {D}}(\partial {{\widetilde{E}}})$ for all $k\in \mathbb {N}_{0}$, we have $|\phi _{k+1}|\leqslant 1$; one uses the fact that $|\phi ^{N}_{0}|\leqslant 1$ and the definition of the interpolants (3.2)–(3.3) to conclude (3.12)$_{7}.$ Let us note that (3.12)$_2$ is a consequence of bounds (3.12)$_{1},$ the embedding $W^{1,2}(\Omega )\hookrightarrow L^{6}(\Omega )$ and the following interpolation

$$\begin{aligned} {[}L^{\infty }(0,T+1;L^{2}(\Omega )),L^{2}(0,T+1;L^{6}(\Omega ))]_{\theta =\frac{3}{5}}=L^{\frac{10}{3}}(Q_{T+1}). \end{aligned}$$

The boundedness (3.12)$_3$ follows by (3.11), assumption (1.6) and the bound of $\left\{ \frac{\xi (\phi ^N_h)}{\alpha ^2}(|M^N|^2-1)^2\right\}$ in $L^\infty (0,T+1;L^1(\Omega ))$. Applying (2.21) we have

$$\begin{aligned} \int _0^{T+1}\left| \int _\Omega \mu ^{N}\right| \leqslant G(T+1) \end{aligned}$$

for a monotone function $G:(0,\infty )\rightarrow (0,\infty )$. Combining the latter bound with (3.12)$_5$ we get

$$\begin{aligned} \{\mu ^N\}\text { is bounded in }L^2(0,T+1;W^{1,2}(\Omega )). \end{aligned}$$

(3.13)

Moreover, by the definition of a time shifted function we get

$$\begin{aligned} \begin{aligned} \{v^N_h\}&\text { is bounded in }&L^\infty (0,T+1;L^2(\Omega )),\\ \{M^N_h\}&\text { is bounded in }&L^\infty (0,T+1;W^{1,2}(\Omega )),\\ \{\phi ^N_h\}&\text { is bounded in}&L^\infty (0,T+1;W^{1,2}(\Omega )). \end{aligned} \end{aligned}$$

(3.14)

All the sequences in (3.14) in the respective norms are bounded by $C(E_{tot}(v_{0}$, $M_{0},\phi _{0})+1)^{\frac{1}{2}}.$ We conclude directly from the definition of the interpolants that

$$\begin{aligned} \{\phi ^N\},\{\phi ^N_h\},\{{{\widetilde{\phi }}}^N\}\subset [-1,1]. \end{aligned}$$

(3.15)

Taking into account the definition of $\rho ^N$ (as defined in (3.4)) and $\rho ^N_h$ it follows that

$$\begin{aligned} \{\rho ^N\},\{\rho ^N_h\}\text { are bounded in }L^\infty (0,T+1;W^{1,2}(\Omega ))\text { and }L^\infty (Q_{T+1}).\nonumber \\ \end{aligned}$$

(3.16)

As a consequence of bounds (3.12)$_{1,2,3,4}$ and (3.13) one has up to subsequences that are not explicitly relabeled

$$\begin{aligned} \begin{aligned} v^N&\rightharpoonup v{} & {} \text { in }L^2(0,T;W^{1,2}(\Omega )),\\ v^N&\rightharpoonup ^* v{} & {} \text { in }L^\infty (0,T;L^2(\Omega )),\\ v^N&\rightharpoonup v{} & {} \text { in }L^\frac{10}{3}(Q_T),\\ M^N&\rightharpoonup ^* M{} & {} \text { in }L^\infty (0,T;W^{1,2}(\Omega )),\\ \phi ^N&\rightharpoonup ^* \phi{} & {} \text { in }L^\infty (0,T;W^{1,2}(\Omega )),\\ \mu ^N&\rightharpoonup \mu{} & {} \text { in }L^2(0,T;W^{1,2}(\Omega )). \end{aligned} \end{aligned}$$

(3.17)

Next we will collect some strong convergence results of the interpolants. Some of these results are already proved in [36]. First following the arguments from [36, (5.31), (5.32) and (5.33)] we obtain

$$\begin{aligned} \phi ^N,\phi ^N_h,{{\widetilde{\phi }}}^N\rightarrow \phi \text { in }L^2(0,T;L^4(\Omega )). \end{aligned}$$

(3.18)

By [36, (5.39) and (5.40)] we have

$$\begin{aligned} M^N\rightarrow M\text { in }L^8(0,T;L^4(\Omega )),\ M^N_h\rightarrow M\text { in }L^2(Q_T). \end{aligned}$$

(3.19)

Moreover, the convergence

$$\begin{aligned} \begin{aligned}&{{\,\textrm{div}\,}}(\xi (\phi ^N_h)\nabla M^N)-\frac{\xi (\phi ^N_h)}{\alpha ^2}(|M^N|^2M^N-M^N_h)\\&\quad \rightharpoonup {{\,\textrm{div}\,}}(\xi (\phi )\nabla M)-\frac{\xi (\phi )}{\alpha ^2}(|M|^2M-M)\text { in }L^2(Q_T) \end{aligned} \end{aligned}$$

(3.20)

follows by [36, (5.58)]. Combining (3.15) with (3.18) we get up to a nonrelabeled subsequence

$$\begin{aligned} \phi ^N\rightarrow \phi \text { in }L^p(Q_T) \text { for all }p\in [1,\infty )\text { and a.e.\ in }Q_T. \end{aligned}$$

(3.21)

To prove this claim one uses the strong convergence of $\phi ^{N}$ from (3.18), boundedness from (3.15) and the following inequalities

$$\begin{aligned} \begin{array}{ll} \displaystyle \Vert \phi ^{N}-\phi \Vert ^{p}_{L^{p}(Q_{T})}&{}\displaystyle \leqslant C \int _{0}^{T}\Vert \phi ^{N}-\phi \Vert ^{(1-\frac{4}{p})p}_{L^{\infty }(\Omega )}\Vert \phi ^{N}-\phi \Vert ^{{4}}_{L^{4}(\Omega )}\\ &{}\displaystyle \leqslant C \int _{0}^{T}\Vert \phi ^{N}-\phi \Vert ^{(1-\frac{4}{p})p}_{L^{\infty }(\Omega )}\Vert \phi ^{N}-\phi \Vert ^{{2}}_{L^{\infty }(\Omega )}\Vert \phi ^{N}-\phi \Vert ^{{2}}_{L^{4}(\Omega )}\\ &{}\displaystyle \leqslant C\Vert \phi ^{N}-\phi \Vert ^{2+(1-\frac{4}{p})p}_{L^{\infty }(Q_{T})}\Vert \phi ^{N}-\phi \Vert ^{2}_{L^{2}(0,T;L^{4}(\Omega ))} \end{array} \end{aligned}$$

for $p\in [4,\infty )$ and $\Vert \cdot \Vert _{L^{p}(\Omega )}\leqslant C\Vert \cdot \Vert _{L^{4}(\Omega )}$ for any $p\in [1,4).$ Taking into account the definition of $\rho ^N$ (we refer to (3.4)) and $\rho ^N_h$ we obtain

$$\begin{aligned} \rho ^N,\rho ^N_h\rightarrow \rho \text { in }L^p(Q_T) \text { for all }p\in [1,\infty )\text { and a.e. in }Q_T. \end{aligned}$$

(3.22)

We now focus on the proof of the compactness of interpolants for the velocity with respect to the topology of a suitable function space. This proof is not straightforward as no uniform bound is available on a sequence of time derivatives of piecewise affine interpolants for the velocity. We investigate the convergence of $\{\widetilde{\rho v}^N\}$. From (3.7) it follows in particular that

$$\begin{aligned} \begin{aligned} \int _0^T&\langle \partial ^-_{t,h}(\rho ^Nv^N),\psi _1\rangle \\ =&\int _0^T\int _\Omega \left( \rho ^N_hv^N\otimes v^N+v^N\otimes J^N\right) \cdot \nabla \psi _1-2\nu (\phi ^{N}_{h})\mathbb {D}v^{N}\cdot \mathbb {D}\psi _1\\&\!+\!\left( \!\left( \frac{\xi (\phi ^N_h)}{\alpha ^2}(|M^N|^2M^N\!-\!M^N_h)\!-\!{{\,\textrm{div}\,}}(\xi (\phi ^N_h)\nabla M^N)\!\right) \nabla M^N\!-\!\nabla \mu ^N\phi ^N_h\right) \!\cdot \!\psi _1 \end{aligned} \end{aligned}$$

(3.23)

for all $\psi _1\in L^8(0,\infty ;V(\Omega ))$. We note that $\left\{ \rho ^N_h v^N\otimes v^N\right\}$ is bounded in $L^2(0,T;L^{\frac{3}{2}}(\Omega ))$ and $\{v^N\otimes J^N\}$ is bounded in $L^\frac{8}{7}(0,T;L^\frac{4}{3}(\Omega )).$ The explanation of achieving these bounds can be found in [6, p. 474]. Further one easily shows that $\{\nabla \mu ^N\phi ^N_h\}$ is bounded in $L^2(Q_{T})$ by using that $\{\phi ^{N}_{h}\}$ is bounded in $L^{\infty }(Q_{T})$ and $\{\mu ^{N}\}$ is bounded in $L^{2}(0,T;W^{1,2}(\Omega )).$ Using (1.7) and (3.12)$_{1},$ $\nu (\phi ^{N}_{h})\mathbb {D}v^{N}$ is bounded in $L^{2}(Q_{T}).$ Moreover, by (3.12)$_{3,6}$ we have the bound on $\left\{ \left( \frac{\xi (\phi ^N_h)}{\alpha ^2}(|M^N|^2\,M^N-M^N_h)-{{\,\textrm{div}\,}}(\xi (\phi ^N_h)\nabla M^N)\right) \nabla M^N\right\}$ in $L^2(0,T;L^1(\Omega ))$. Since $\partial _{t}\widetilde{\rho v}^{N}=\partial ^{-}_{t,h}(\rho ^{N}v^{N})$, one uses (3.23) and the fact that the Leray projector $\mathbb {P}_{{{\,\textrm{div}\,}}}$ commutes with the time derivative to infer that

$$\begin{aligned} \left\{ \partial _t\mathbb {P}_{{{\,\textrm{div}\,}}}\left( \widetilde{\rho v}^N\right) \right\} \text { is bounded in }L^\frac{8}{7}(0,T;(V(\Omega ))'). \end{aligned}$$

(3.24)

Taking into account the uniform bound on $\left\{ \mathbb {P}_{{{\,\textrm{div}\,}}}\left( \widetilde{\rho v}^N\right) \right\}$ in $L^\infty (0,T;L^2(\Omega ))$, which follows from the continuity of $\mathbb {P}_{{{\,\textrm{div}\,}}}$, (3.12)$_{1}$ and (3.16), the Aubin-Lions lemma gives the compactness of $\left\{ \mathbb {P}_{{{\,\textrm{div}\,}}}\left( \widetilde{\rho v}^N\right) \right\}$ with respect to the strong topology of $L^2(0,T;(W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))')$. Moreover, in view of (3.17)$_{1}$ and the almost everywhere convergence (3.22) we can follow line by line the arguments presented in [49, pp. 90–91, (3.95)] to conclude that

$$\begin{aligned} \mathbb {P}_{{{\,\textrm{div}\,}}}\left( \widetilde{\rho v}^N\right) \rightharpoonup \mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v)\text { in }L^2(0,T;L^2_{{{\,\textrm{div}\,}}}(\Omega )). \end{aligned}$$

(3.25)

Consequently, for a nonrelabeled subsequence we obtain

$$\begin{aligned} \mathbb {P}_{{{\,\textrm{div}\,}}}\left( \widetilde{\rho v}^N\right) \rightarrow \mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v)\text { in }L^2(0,T;(W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))'). \end{aligned}$$

(3.26)

Since

$$\begin{aligned} \begin{aligned} \mathbb {P}_{{{\,\textrm{div}\,}}}\left( \widetilde{\rho v}^N(t)\right) -\mathbb {P}_{{{\,\textrm{div}\,}}}\left( \rho ^N(t)v^N(t)\right) =(t-(k+1)h)\partial _{t}\mathbb {P}_{{{\,\textrm{div}\,}}}\left( \widetilde{\rho v}^{N}\right) (t) \end{aligned} \end{aligned}$$

for all $t\in [kh,(k+1)h)$, $k\in \mathbb {N}_0$ and $|t-(k+1)h|\leqslant h$ we infer using (3.24) that

$$\begin{aligned} \begin{array}{l} \mathbb {P}_{{{\,\textrm{div}\,}}}\left( \widetilde{\rho v}^N\right) -\mathbb {P}_{{{\,\textrm{div}\,}}}\left( \rho ^Nv^N\right) \rightarrow 0\text { in } L^\frac{8}{7}(0,T;(V(\Omega ))'). \end{array} \end{aligned}$$

(3.27)

Next we consider the following combination of interpolation and duality

$$\begin{aligned} (W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))'=\left( (L^2_{{{\,\textrm{div}\,}}}(\Omega ), V(\Omega ))_{\frac{1}{2},2}\right) ' =\left( L^2_{{{\,\textrm{div}\,}}}(\Omega ), (V(\Omega ))'\right) _{\frac{1}{2},2}, \end{aligned}$$

where the first equality is a special case of [1, (5.2.17)] and the second one follows by [14, Theorem 3.7.1]. Employing the inequality that corresponds to the latter interpolation we obtain

$$\begin{aligned} \begin{aligned}&\left\| \mathbb {P}_{{{\,\textrm{div}\,}}}(\widetilde{\rho v}^{N})-\mathbb {P}_{{{\,\textrm{div}\,}}}(\rho ^{N}v^{N})\right\| _{L^{2}(0,T;(W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))')}\\&\quad \leqslant C\left\| \mathbb {P}_{{{\,\textrm{div}\,}}}(\widetilde{\rho v}^{N})-\mathbb {P}_{{{\,\textrm{div}\,}}}(\rho ^{N}v^{N})\right\| _{L^{\infty }(0,T;L^{2}_{{{\,\textrm{div}\,}}}(\Omega ))}^{\frac{1}{2}}\\&\qquad \times \left\| \mathbb {P}_{{{\,\textrm{div}\,}}}(\widetilde{\rho v}^{N})-\mathbb {P}_{{{\,\textrm{div}\,}}}(\rho ^{N}v^{N})\right\| _{L^{1}(0,T;(V(\Omega ))')}^{\frac{1}{2}}. \end{aligned} \end{aligned}$$

Combining the latter inequality with the bound on $\left( \mathbb {P}_{{{\,\textrm{div}\,}}}(\widetilde{\rho v}^{N})-\mathbb {P}_{{{\,\textrm{div}\,}}}(\rho ^{N}v^{N})\right)$ in $L^{\infty }(0,T;L^{2}(\Omega ))$ and the convergence (3.27) we furnish that

$$\begin{aligned} \begin{array}{l} \mathbb {P}_{{{\,\textrm{div}\,}}}(\widetilde{\rho v}^N)-\mathbb {P}_{{{\,\textrm{div}\,}}}(\rho ^Nv^N)\rightarrow 0\,\,\text{ in }\,\, L^{2}(0,T;((W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))'). \end{array} \end{aligned}$$

(3.28)

The convergences (3.26) and (3.28) together furnish that

$$\begin{aligned} \mathbb {P}_{{{\,\textrm{div}\,}}}(\rho ^Nv^N)\rightarrow \mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v)\text { in }L^2(0,T; (W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))'). \end{aligned}$$

(3.29)

Next, by (3.22) and (3.17)$_3$ we conclude $(\rho ^N)^\frac{1}{2}v^N\rightharpoonup \rho ^\frac{1}{2}v$ in $L^2(Q_T)$. Furthermore, combining (3.29) and (3.17)$_1$ it follows that

$$\begin{aligned} \begin{array}{ll} \displaystyle \int _0^T\int _\Omega \rho ^N|v^N|^2&{}\displaystyle =\int _0^T\langle \mathbb {P}_{{{\,\textrm{div}\,}}}(\rho ^Nv^N),v^N\rangle _{(W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))',W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )}\\ &{} \displaystyle \rightarrow \int _0^T\langle \mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v),v\rangle _{(W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega ))',W^{1,2}_{0,{{\,\textrm{div}\,}}}(\Omega )} \displaystyle =\int _0^T\int _\Omega \rho |v|^2 \end{array} \end{aligned}$$

(3.30)

implying $\Vert (\rho ^N)^\frac{1}{2}v^N\Vert _{L^2(Q_T)}\rightarrow \Vert \rho ^\frac{1}{2}v\Vert _{L^2(Q_T)}$. Hence passing to a nonrelabeled subsequence one has

$$\begin{aligned} (\rho ^N)^\frac{1}{2}v^N\rightarrow \rho ^\frac{1}{2}v\text { in }L^2(Q_T) \text { and a.e. in }Q_T. \end{aligned}$$

(3.31)

Moreover, (3.22) and the existence of a positive lower bound on $\{\rho ^N\}$, obtained by a similar argument as in Remark 1.2, imply $(\rho ^N)^{-\frac{1}{2}}\rightarrow \rho ^{-\frac{1}{2}}$ a.e. in $Q_T$. As (3.12)$_2$ ensures the equiintegrability of the sequence $\{|v^N|^q\}$, $q\in [1,\frac{10}{3})$ we conclude by the Vitali convergence theorem

$$\begin{aligned} v^N=(\rho ^N)^{-\frac{1}{2}}(\rho ^N)^\frac{1}{2}v^N\rightarrow v\text { in }L^q(Q_T),\ q\in [1,\tfrac{10}{3}). \end{aligned}$$

(3.32)

In particular the strong convergence $v^{N}\rightarrow v$ in $L^{2}(Q_{T})$ (as a consequence of (3.32)), the boundedness of $(v^{N}-v)$ in $L^{2}(0,T;L^{6}(\Omega ))$ and the following interpolation inequality

$$\begin{aligned} \displaystyle \Vert v^{N}-v\Vert _{L^{2}(0,T;L^{4}(\Omega ))}\leqslant C\Vert v^{N}-v\Vert _{L^{2}(0,T;L^{6}(\Omega ))}^{\frac{3}{4}}\Vert v^{N}-v\Vert _{L^{2}(Q_{T})}^{\frac{1}{4}} \end{aligned}$$

provides

$$\begin{aligned} \begin{array}{l} v^N\rightarrow v\,\,\text { in }\,\,L^{2}(0,T;L^{4}(\Omega )). \end{array} \end{aligned}$$

(3.33)

The last important convergence is

$$\begin{aligned} M^N\rightarrow M\text { in }L^2(0,T;W^{1,2}(\Omega )). \end{aligned}$$

(3.34)

The convergence (3.34) is crucial in order to pass to the limit in the term containing $|\nabla M^{N}|^{2}$ in (3.10) and the term comprising of $\nabla M^{N}$ in the momentum equation (3.7). The proof of (3.34) relies on the monotone structure of ${{\,\textrm{div}\,}}(\xi (\phi ^{N}_{h})\nabla M^{N})$ and can be done by following the arguments used to show [36, (5.41),  Section 5.1.2]. In order to do so, the strong convergence of $M^{N}$ to M in $L^{8}(0,T;L^{4}(\Omega ))$ and $v^{N}$ to v in $L^{2}(0,T;L^{4}(\Omega ))$ are used. Since these convergences are available in the present scenario (we refer to (3.19) and (3.33)), we face no particular difficulty to follow line by line the proof of [36, (5.41),  Section 5.1.2].

3.1.2 Some uniform estimates on $M^{N}$ and $\phi ^{N}$

In this section we will obtain further uniform estimates that involve the integrability of $\nabla M^N$ w.r.t. spatial variables for an exponent greater than 2, the integrability of the second gradient of $\phi ^{N}$ and of ${\widetilde{\Psi }}'_{0}(\phi ^{N})$ w.r.t. spatial variables for an exponent greater than 1 depending only on the energy estimate (3.11) for the interpolants. These improved estimates will aid in recovering the weak formulation of Cahn–Hilliard equations. In that direction we will make use of an abstract elliptic regularity result from [32]. The central result of this section is Lemma 3.2 which will be proved by using the following result.

Lemma 3.1

Let $\Omega$ be a bounded domain of class $C^{1}$ in $\mathbb {R}^{d},$ $d\geqslant 2$. Let ${\widetilde{\xi }}:\Omega \rightarrow \mathbb {R}^{+}$ be a bounded, measurable function satisfying

$$\begin{aligned} 0<c_1\leqslant {\widetilde{\xi }}(\cdot )\leqslant c_2\text { on }\Omega , \text { for some}\,\, c_{1},c_{2}>0. \end{aligned}$$

(3.35)

Then there is $2<p<3$ such that any solution $M=(M_{1},M_{2},M_{3})\in W^{1,2}(\Omega )$ of the following elliptic problem with homogeneous Neumann boundary condition

$$\begin{aligned} \begin{aligned} {{\,\textrm{div}\,}}({{\widetilde{\xi }}\nabla M})&=g{} & {} \text { in }\Omega ,\\ \partial _{n}M&=0{} & {} \text { on }\partial \Omega , \end{aligned} \end{aligned}$$

(3.36)

where $g\in L^s(\Omega )$ with $s\geqslant \frac{p'd}{d(p'-1)+p'}$ ($p'$ being the conjugate exponent of p) satisfies

$$\begin{aligned} \begin{array}{l} \Vert M\Vert _{W^{1,p}(\Omega )}\leqslant C\left( \Vert g\Vert _{L^s(\Omega )}+\Vert M\Vert _{W^{1,2}(\Omega )}\right) . \end{array} \end{aligned}$$

(3.37)

The constant C depends on d,  $c_{1}$ and $c_{2}$ and the domain $\Omega .$

Proof

The result stated in Lemma 3.1 is a special case of the more general result in [32, Remark 13]. For the convenience of the readers we present the proof. First, we rewrite (3.36) component wise in the form

$$\begin{aligned} \begin{aligned} -\sum \limits _{i=1}^{d}\partial _{i}({{\widetilde{\xi }}}\partial _{i}M_{k})+M_{k}=&-g_{k}+M_{k}{} & {} \text { in }\Omega ,\\ \partial _{n}\,M_k=&0{} & {} \text { on }\partial \Omega , \end{aligned} \end{aligned}$$

(3.38)

for $k\in \{1,2,3\}$. The operator $A:W^{1,2}(\Omega )\rightarrow (W^{1,2}(\Omega ))^{'}$ appearing in the weak form of the latter problem is defined as

$$\begin{aligned} \langle Au,w\rangle =\int _\Omega a_{ij}\partial _i u\partial _j w+uw\text { for }u,w\in W^{1,2}(\Omega ) \end{aligned}$$

(3.39)

with the matrix $(a_{ij})_{i,j\in \{1,..,d\}}$ given by

$$\begin{aligned} a_{ij}(x)=\left\{ \begin{array}{ll} {\widetilde{\xi }}(x)&{}\quad \text{ when }\quad i=j,\\ 0&{}\quad \text{ when }\quad i\ne j. \end{array}\right. \end{aligned}$$

Taking into account (3.35) we infer $a_{ij}(\cdot )\in L^{\infty }(\Omega )$ and the ellipticity condition

$$\begin{aligned} \sum _{i,j=1}^{d}a_{ij}(x)\theta _{i}\theta _{j}\geqslant c_{1}|\theta |^{2}\text { for any }\theta =(\theta _{1},\ldots ,\theta _{d})\in \mathbb {R}^d. \end{aligned}$$

In order to apply [32, Theorem 1], we note that as $\Omega$ is of class $C^1$, assumptions of [32, Theorem 1] are fulfilled, cf. [32, Remark 1 and 7]. Hence we conclude the existence of some $2<p<3$ such that the operator A maps $W^{1,p}(\Omega )$ onto $(W^{1,p'}(\Omega ))^{'}$, indeed $p'$ is the conjugate exponent of p. The upper bound on p can be found in [32, Theorem 1], from which one also infers the inequality

$$\begin{aligned} \Vert A^{-1}f\Vert _{W^{1,p}(\Omega )}\leqslant c\Vert f\Vert _{{(W^{1,p'}(\Omega ))^{'}}} \end{aligned}$$

(3.40)

due to the linearity of $A^{-1}$. We notice that the constant c depends on $\Omega$ and constants from (3.35). In order to conclude (3.37), we employ (3.40) with f being the right hand side of (3.38)$_1$ and use the embedding $L^s(\Omega )\hookrightarrow (W^{1,p'}(\Omega ))'$ where $s\geqslant \frac{p'd}{d(p'-1)+p'}$ for $2<p<3.$ $\square$

Lemma 3.2

Let $(v^{N},M^{N},\phi ^{N},\mu ^{N})$ be the interpolants defined by (3.2)–(3.3), satisfying (3.7)–(3.10) and the energy estimate (3.11). Then $M^{N}$ satisfies

$$\begin{aligned} \begin{array}{l} \Vert M^{N}\Vert _{L^{2}(0,T;W^{1,p}(\Omega ))}\leqslant CE^{+}_{tot}(v_0,M_0,\phi _0)^{\frac{3}{2}}, \end{array} \end{aligned}$$

(3.41)

where $E^{+}_{tot}(v_0,M_0,\phi _0)=\left( E_{tot}(v_0,M_0,\phi _0)+1\right) ,$ for some $2<p<3$ and the positive constant C might depend on $\alpha$, $c_{1},$ $c_{2},$ $c_{3}$, cf. (1.6), Sobolev embedding constants and the domain $\Omega .$

Further $\phi ^N$ and ${\widetilde{\Psi }}'_{0}(\phi ^N)$ satisfy

$$\begin{aligned} \Vert \phi ^{N}\Vert _{L^{2}(0,T;W^{2,q}(\Omega ))}+\Vert {\widetilde{\Psi }}'_{0}(\phi ^{N})\Vert _{L^2(0,T;L^q(\Omega ))}\leqslant CE^{+}_{tot}(v_0,M_0,\phi _0)^3. \end{aligned}$$

(3.42)

where $1<q=\frac{2p}{p+2}<2$ and the positive constant C in (3.42) depends on $\alpha$, $c_{1},$ $c_{2},$ $c_{3},$ Sobolev embedding constants and the domain $\Omega .$

Proof

In order to prove (3.41) we consider (3.8) for a.e. $t\in (0,T)$ as the elliptic problem

$$\begin{aligned} \begin{aligned} {{\,\textrm{div}\,}}(\xi (\phi ^N_h)\nabla M^N){=}&\partial ^{-}_{t,h}M^N\!+\!(v^N\!\cdot \!\nabla )M^N\!{+}\!\frac{\xi (\phi ^N_h)}{\alpha ^2}(|M^N|^2M^N\!{-}\!M^N_h){} & {} \text { in }\Omega ,\\ \partial _{n}M^{N}{=}&0{} & {} \text { on }\partial \Omega , \end{aligned}\nonumber \\ \end{aligned}$$

(3.43)

which is possible since, in view of (3.12), all the terms involved in (3.43) are defined a.e. We focus on the estimate of the right hand side in (3.43)$_1$. Using the Hölder inequality, the Sobolev embedding and (3.12) we get

$$\begin{aligned} \begin{aligned}&\Vert (v^N\cdot \nabla ) M^N\Vert _{L^2(0,T;L^\frac{3}{2}(\Omega ))}\\&\quad \leqslant \Vert v^N\Vert _{L^2(0,T;L^6(\Omega ))}\Vert \nabla M^N\Vert _{L^\infty (0,T;L^2(\Omega ))}\\&\quad \leqslant C\Vert v^N\Vert _{L^2(0,T;W^{1,2}(\Omega ))}\Vert \nabla M^N\Vert _{L^\infty (0,T;L^2(\Omega ))} \\&\quad \leqslant C E^{+}_{tot}(v_0,M_0,\phi _0). \end{aligned} \end{aligned}$$

(3.44)

Since $L^{2}(0,T;W^{1,2}(\Omega ))$ is dense in $L^{2}(0,T;L^{3}(\Omega )),$ one can choose test functions $\psi _{2}\in L^{2}(0,T;L^{3}(\Omega ))$ in (3.8) and use (3.12)$_{1,6}$ to compute the following

$$\begin{aligned} \begin{aligned}&\Vert \partial ^{-}_{t,h}M^N\Vert _{L^2(0,T;L^\frac{3}{2}(\Omega ))}\\&\quad \displaystyle =\sup _{\{\psi _{2}\in L^{2}(0,T;L^{3}(\Omega ))\;|\;\Vert \psi _{2}\Vert _{L^{2}(0,T;L^{3}(\Omega ))}\leqslant 1\}}\int _{0}^{T}\int _{\Omega }\partial ^{-}_{t,h}M^N\cdot \psi _{2}\\&\quad \displaystyle \leqslant C\bigg (\Vert v^N\Vert _{L^2(0,T;W^{1,2}(\Omega ))}\Vert \nabla M^N\Vert _{L^\infty (0,T;L^2(\Omega ))}\\&\qquad +\left\| {{\,\textrm{div}\,}}(\xi (\phi ^N_h)\nabla M^N)-\frac{\xi (\phi ^N_h)}{\alpha ^2}(|M^N|^2M^N-M^N_h)\right\| _{L^{2}(Q_{T})}\bigg )\\&\quad \displaystyle \leqslant C\left( E^{+}_{tot}(v_0,M_0,\phi _0)+E^{+}_{tot}(v_0,M_0,\phi _0)^\frac{1}{2}\right) . \end{aligned} \end{aligned}$$

(3.45)

Moreover, using bounds on $\xi$ in (1.6), the Hölder inequality, the Sobolev embedding and the definition of $M^N_h$ we arrive at

$$\begin{aligned} \begin{aligned}&\left\| \frac{\xi (\phi ^N_h)}{\alpha ^2}(|M^N|^2M^N-M^N_h)\right\| _{L^2(0,T;L^\frac{3}{2}(\Omega ))}\\&\quad \leqslant C\left( \Vert M^N\Vert ^3_{L^6(0,T;L^\frac{9}{2}(\Omega ))}+\Vert M^N_h\Vert _{L^2(0,T;L^\frac{3}{2}(\Omega ))}\right) \\&\quad \leqslant C\left( \Vert M^N\Vert ^3_{L^\infty (0,T;W^{1,2}(\Omega ))}+\Vert M^N_h\Vert _{L^\infty (0,T;W^{1,2}(\Omega ))}\right) \\&\quad \leqslant C\left( E^{+}_{tot}(v_0,M_0,\phi _0)^\frac{3}{2}+E^{+}_{tot}(v_0,M_0,\phi _0)^\frac{1}{2}\right) . \end{aligned} \end{aligned}$$

(3.46)

Applying Lemma 3.1 with $s=\frac{3}{2}$ ( the value $s=\frac{3}{2}$ is admissible since $p<3$ implying $p'>\frac{3}{2}$ $\Rightarrow$ $\frac{3p'}{3(p'-1)+p'}<\frac{3}{2},$ where the dimension $d=3.$ Similar justification holds when $d=2.$ ) to (3.43) we obtain

$$\begin{aligned} \begin{aligned} \Vert M^N\Vert _{W^{1,p}(\Omega )}\leqslant&C\Biggl (\Vert \partial ^{-}_{t,h}M^N\Vert _{L^\frac{3}{2}(\Omega )}+\Vert (v^N\cdot \nabla ) M^N\Vert _{L^\frac{3}{2}(\Omega )}\\&+\left\| \frac{\xi (\phi ^N_h)}{\alpha ^2}(|M^N|^2M^N-M^N_h)\right\| _{L^\frac{3}{2}(\Omega )}+\Vert M^N\Vert _{W^{1,2}(\Omega )}\Biggr ) \end{aligned} \end{aligned}$$

a.e. in (0, T) with some $2<p<3$ and the constant C independent of the time variable. We combine the latter inequality with (3.44), (3.45), (3.46) and the Young inequality to conclude (3.41).

To show (3.42), we will use (3.10). Using (3.12)$_{4},$ (3.13) and (3.14)$_{3}$ one can bound the first two terms appearing in the left hand side of (3.10) in $L^{2}(0,T;L^{6}(\Omega ))$ by a constant multiple of $E_{tot}(v_0,M^N_0,\phi ^N_0)^{\frac{1}{2}}.$ Further since $\Vert H_{0}(\phi ^N,\phi ^N_h)\Vert _{L^{\infty }(Q_{T})}\leqslant c_{3}$ and by (3.12)$_3$ along with (3.41), one has

$$\begin{aligned} \begin{aligned}&\tfrac{1}{2}\left\| H_{0}(\phi ^N,\phi ^N_h)|\nabla M^N|^{2}\right\| _{L^{2}(0,T;L^{\frac{2p}{p+2}}(\Omega ))}\\&\quad \leqslant C\Vert \nabla M^{N}\Vert _{L^{2}(0,T;L^{p}(\Omega ))}\Vert \nabla M^{N}\Vert _{L^{\infty }(0,T;L^{2}(\Omega ))}\\&\quad \leqslant CE^{+}_{tot}(v_0,M_0,\phi _0), \end{aligned} \end{aligned}$$

(3.47)

and

$$\begin{aligned} \begin{aligned}&\tfrac{1}{4\alpha ^2}\left\| H_{0}(\phi ^N,\phi ^N_h)(|M^N|^{2}-1)^{2}\right\| _{L^{2}(0,T;L^{\frac{2p}{p+2}}(\Omega ))} \\&\quad \leqslant C\left( \Vert M^{N}\Vert ^{3}_{L^{\infty }(0,T;L^{6}(\Omega ))}\Vert M^{N}\Vert _{L^{2}(0,T;L^{p}(\Omega ))}+1\right) \\&\quad \leqslant C\left( E^{+}_{tot}(v_0,M_0,\phi _0)^{\frac{3}{2}}E^{+}_{tot}(v_0,M_0,\phi _0)^{\frac{3}{2}}+1\right) \\&\quad \leqslant CE^{+}_{tot}(v_0,M_0,\phi _0)^{3}. \end{aligned} \end{aligned}$$

(3.48)

In both of (3.47) and (3.48) the positive constant C might depend on $\alpha$, $c_{1},$ $c_{2},$ $c_{3},$ Sobolev embedding constants and $|\Omega |.$ From the discussion above (in particular the inequalities (3.47) and (3.48)) one infers from (3.10)

$$\begin{aligned} \begin{aligned} -\eta \Delta \phi ^N+{\widetilde{\Psi }}'_{0}(\phi ^N)=&f(M^{N},\phi ^{N},\phi ^{N}_{h},\mu ^{N}){} & {} \text{ in } \Omega ,\\ \partial _{n}\phi ^{N}=&0{} & {} \text{ on } \partial \Omega \end{aligned} \end{aligned}$$

(3.49)

a.e. in (0, T) with

$$\begin{aligned} \left\| f(M^{N},\phi ^{N},\phi ^{N}_{h},\mu ^{N})\right\| _{L^{2}(0,T;L^{\frac{2p}{p+2}}(\Omega ))}\leqslant CE^{+}_{tot}(v_0,M_0,\phi _0)^3, \end{aligned}$$

where C might depend on $\alpha$, $c_{1},$ $c_{2},$ $c_{3},$ Sobolev embedding constants and the domain $\Omega .$

Finally, applying the inequality (2.8) from Proposition 2.1 to (3.49) and the Young inequality again we obtain (3.42). $\square$

3.1.3 Additional convergences of $\{M^N\}$, $\{\phi ^N\}$, $\{\Psi '(\phi ^N)\}$

In view of estimate (3.41), we immediately obtain that for some $p>2$ we have up to a nonrelabeled subsequence

$$\begin{aligned} M^N\rightharpoonup M\text { in }L^2(0,T;W^{1,p}(\Omega )), \end{aligned}$$

where M comes from (3.17). This concludes (1.15)$_1$. Similarly, by (3.42) we have up to a nonrelabeled subsequence

$$\begin{aligned} \phi ^N\rightharpoonup \phi \text { in }L^2(0,T;W^{2,\frac{2p}{p+2}}(\Omega )), \end{aligned}$$

(3.50)

proving (1.15)$_2$. The next task is to show that up to a nonrelabeled subsequence

$$\begin{aligned} {\widetilde{\Psi }}'_0(\phi ^N)\rightharpoonup {\widetilde{\Psi }}'_0(\phi )\text { in }L^2(0,T;L^{\frac{2p}{p+2}}(\Omega )), \end{aligned}$$

(3.51)

from which (1.15)$_3$ follows. We observe that the estimate of ${\widetilde{\Psi }}'_0(\phi ^N)$ in (3.42) implies that

$$\begin{aligned} {\widetilde{\Psi }}'_0(\phi ^N)\rightharpoonup \zeta \text { in }L^2(0,T;L^{\frac{2p}{p+2}}(\Omega )). \end{aligned}$$

Hence we have to identify $\zeta$. Let us begin with showing that

$$\begin{aligned} {{\widetilde{\Psi }}}'_0(\phi ^N)\rightarrow {{\widetilde{\Psi }}}'_0(\phi )\text { a.e. in }Q_T. \end{aligned}$$

(3.52)

To this end we adopt arguments devised in the context of the Cahn–Hillard equations with a logarithmic free energy, see [20, p. 1510], and developed for the case of the Navier–Stokes–Cahn–Hilliard system with a singular potential, see [25, p. 285]. We define for arbitrary but fixed $\delta \in (0,1)$ the quantity $a_\delta =\min \left\{ {{\widetilde{\Psi }}}'_0(1-\delta ),-{{\widetilde{\Psi }}}'_0(-1+\delta )\right\}$. Then we have $a_\delta \leqslant |{{\widetilde{\Psi }}}'_0(s)|$ for $1>|s|>1-\delta$ as ${{\widetilde{\Psi }}}'_0$ is nondecreasing. Hence we obtain

$$\begin{aligned} a_\delta \left| \left\{ (t,x)\in Q_T: 1>|\phi ^N(t,x)|>1-\delta \right\} \right| \leqslant \int _{Q_T}|{{\widetilde{\Psi }}}'_0(\phi ^N)|\leqslant c \end{aligned}$$

by (3.42) and Hölder’s inequality. Combining the pointwise convergence $\phi ^N\rightarrow \phi$ from (3.21) and the Fatou Lemma with the latter inequality we conclude

$$\begin{aligned} \begin{aligned}&|\{(t,x)\in Q_T:1\geqslant |\phi (t,x)|\geqslant 1-\delta \}|\\&\quad \leqslant \liminf _{N\rightarrow \infty }\left| \left\{ (t,x)\in Q_T:1>|\phi ^N(t,x)|>1-\delta \right\} \right| \leqslant ca_\delta ^{-1} \end{aligned} \end{aligned}$$

(3.53)

for any $\delta \in (0,1)$. Taking into account assumption (1.8)$_1$ it follows that $a_\delta \rightarrow \infty$ as $\delta \rightarrow 0_+$. Hence the limit passage $\delta \rightarrow 0_+$ in (3.53) yields

$$\begin{aligned} |\{(t,x)\in Q_T:|\phi (t,x)|=1\}|=0, \end{aligned}$$

in other words $|\phi |<1$ a.e. in $Q_T$. This bound, the pointwise convergence $\phi ^N\rightarrow \phi$ from (3.21) and the assumed regularity ${{\widetilde{\Psi }}}'_0\in C^1((-1,1))$, cf. Assumption 1.1, imply (3.52). Having (3.52) and the bound on $\{{{\widetilde{\Psi }}}'_0(\phi ^N)\}$ from (3.42) at hand we apply the Vitali convergence theorem to conclude that ${{\widetilde{\Psi }}}'_0(\phi ^N)\rightarrow {{\widetilde{\Psi }}}'_0(\phi )$ in $L^1(Q_T)$. Hence we have $\zeta ={{\widetilde{\Psi }}}'_0(\phi )$ and (3.51) is proved.

The convergence (3.51) along with the fact that $\phi \in [-1,1]$ and (2.1)-(2.2) in particular imply that $\Psi '(\phi )\in L^{2}(0,T;L^{\frac{2p}{p+2}}(\Omega )).$

3.1.4 The energy inequality for the weak solution

This section is devoted to the proof of the fact that the quadruple $(v,M,\phi ,\mu )$ obtained as limits of interpolants (we refer to (3.17)) satisfies (1.16) for all $t\in (0,T),$ where $E_{tot}$ is as defined in (1.17). To this end we take into account (3.17), (3.18), (3.19), (3.22), (3.33) and (3.34) and select subsequences that will not be relabeled such that for a.e. $t\in (0,T)$

$$\begin{aligned} \begin{aligned} \rho ^N(t)&\rightarrow \rho (t){} & {} \text { in }L^2(\Omega ),\\ v^N(t)&\rightarrow v(t){} & {} \text { in }L^4(\Omega ),\\ \nabla M^N(t)&\rightarrow \nabla M(t){} & {} \text { in }L^2(\Omega ),\\ M^N(t)&\rightarrow M(t){} & {} \text { in }L^4(\Omega ),\\ \nabla \phi ^N(t)&\rightharpoonup \nabla \phi (t){} & {} \text { in }L^2(\Omega ),\\ \phi ^N(t)&\rightarrow \phi (t){} & {} \text { in }L^2(\Omega )\text { and a.e. in }\Omega . \end{aligned} \end{aligned}$$

(3.54)

We want to show that for a.e. $t\in (0,T)$

$$\begin{aligned} E_{tot}(v(t),M(t),\phi (t))\leqslant \liminf _{N\rightarrow \infty }E_{tot}(v^N(t),M^N(t),\phi ^N(t)). \end{aligned}$$

(3.55)

We argue as in [36, Section 5.2, pp. 28–29] and focus only on the terms from $E_{tot}(v^N(t),M^N(t),\phi ^N(t))$ that are not treated in [36]. We fix $t\in (0,T)$, in which convergences from (3.54) are available. By (3.54)$_{1,2}$ we get

$$\begin{aligned} \lim _{N\rightarrow \infty }\frac{1}{2}\int _\Omega \rho ^N(t)|v^N(t)|^2=\frac{1}{2}\int _\Omega \rho (t)|v(t)|^2. \end{aligned}$$

Since ${\widetilde{\Psi }}_{0}\in C([-1,1]),$ we obtain by using (3.54)$_{6}$ and the Lebesgue dominated convergence theorem

$$\begin{aligned} \begin{aligned} \lim _{N\rightarrow \infty }\int _\Omega {{\widetilde{\Psi }}}(\phi ^N(t))=&\lim _{N\rightarrow \infty }\int _\Omega \left( {{\widetilde{\Psi }}}_0(\phi ^N(t))-\frac{\kappa }{2}(\phi ^N(t))^2\right) \\ =&\int _\Omega \left( \widetilde{\Psi }_0(\phi (t))-\frac{\kappa }{2}(\phi (t))^2\right) =\int _\Omega \widetilde{\Psi }(\phi (t)). \end{aligned} \end{aligned}$$

The remaining details for the proof of (3.55) can be found in [36, Section 5.2, pp. 28–29]. Applying the convergences from (3.1) we conclude $E_{tot}(v_0,M^N_0,\phi ^N_0)$ $\rightarrow E_{tot}(v_0,M_0,\phi _0)$ in a straightforward way. Hence to conclude (1.16) it suffices to combine (3.55), the fact that

$$\begin{aligned} \displaystyle \sqrt{2\nu (\phi ^{N}_{h})}\mathbb {D}(v^{N})\rightharpoonup \sqrt{2\nu (\phi )}\mathbb {D}(v)\,\,\text{ in }\quad L^{2}(Q_{T}), \end{aligned}$$

which can be proved as in [36, eq. (5.57)], the weak lower semicontinuity of norms with (3.17)$_{6}$ and (3.20).

3.1.5 Continuity with respect to time of $v,M,\phi$

This section aims to show that some of the limit functions obtained in previous sections are continuous w.r.t. time variable in a certain sense. Namely, we show

$$\begin{aligned} \begin{aligned} \rho v\in&C_w([0,T];L^2(\Omega )),\\ v\in&C_{w}([0,T];L^2(\Omega )),\\ M\in&C_w([0,T];W^{1,2}(\Omega )),\\ M\in&C([0,T];L^2(\Omega )),\\ \phi \in&C_w([0,T];W^{1,2}(\Omega )),\\ \phi \in&C([0,T];L^2(\Omega )). \end{aligned} \end{aligned}$$

(3.56)

First, for the proof of (3.56)$_{3,4,5,6}$ we refer to [36, Section 5.3, p. 29]. Let us next prove (3.56)$_{1}.$ As $\rho \in L^\infty (Q_T)$ and $v\in L^\infty (0,T;L^2(\Omega ))$, we have

$$\begin{aligned} \begin{array}{l} \displaystyle \rho v\in L^\infty (0,T;L^2(\Omega )). \end{array} \end{aligned}$$

(3.57)

Next in view of (3.24) one has up to a nonrelabeled subsequence

$$\begin{aligned} \begin{array}{l} \displaystyle \partial _t\mathbb {P}_{{{\,\textrm{div}\,}}}(\widetilde{\rho v}^N)\rightharpoonup \partial _t\mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v)\,\,\text{ in }\,\, L^\frac{8}{7}(0,T;(V(\Omega ))') \end{array} \end{aligned}$$

(where the identification of the limit follows from (3.25)). Then the fact that $\partial _t\mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v)\in L^\frac{8}{7}(0,T;(V(\Omega ))')$ implies $\mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v)\in C([0,T];(V(\Omega ))')$. This along with (3.57) renders

$$\begin{aligned} \begin{array}{l} \mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v)\in C_{w}([0,T];L^{2}_{{{\,\textrm{div}\,}}}(\Omega )) \end{array} \end{aligned}$$

(3.58)

by using [48, Ch. III, Lemma 1.4].

Next using the definition (1.4)–(1.5) of the Leray projector $\mathbb {P}_{{{\,\textrm{div}\,}}}$ we write

$$\begin{aligned} \begin{array}{l} \displaystyle \rho v=\mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v)+\nabla p, \end{array} \end{aligned}$$

(3.59)

where $p(t)\in W^{1,2}(\Omega ),$ $\displaystyle \int _{\Omega } p(t)=0$ and p(t) solves the weak Neumann problem (1.5). Now one can follow the arguments used in [6, Section 5.2, pp. 475–476] to show that $\nabla p\in C_{w}([0,T];L^{2}(\Omega )).$ This along with (3.58) furnishes the proof of (3.56)$_{1}.$

Finally, we wish to show (3.56)$_{2}.$ By definition one needs to prove $v(\cdot ,t_{n})\rightharpoonup v(\cdot ,t)$ in $L^{2}(\Omega )$ for any sequence $\{t_n\}\subset [0,T]$ such that $t_{n}\rightarrow t.$ In view of the non-degeneracy of $\rho ,$ one first infers from (3.59)

$$\begin{aligned} v(\cdot ,t)=\frac{1}{\rho (\cdot ,t)}\mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v)(\cdot ,t)+\frac{1}{\rho (\cdot ,t)}\nabla p(\cdot ,t), \end{aligned}$$

(with this definition one also defines v in a set of measure zero, so that v is defined everywhere in [0, T]) uses (3.56)$_{1},$ $\nabla p\in C_{w}([0,T];L^{2}(\Omega ))$ and the fact that $\rho \in C([0,T];L^{2}(\Omega ))$ (which follows from (3.56)$_{6}$) to show that $v(\cdot ,t_{n})\rightharpoonup v(\cdot ,t)$ in $L^{1}(\Omega ).$ Finally, since $v(\cdot ,t_{n})$ is uniformly bounded in $L^{2}(\Omega ),$ one concludes that $v(\cdot ,t_{n})\rightharpoonup v(\cdot ,t)$ in $L^{2}(\Omega )$ and thereby finishing the proof of (3.56)$_{2}.$

3.2 Recovering the weak formulations

In this section we verify that the quadruple $(v,M,\phi ,\mu )$ satisfies the formulation of the problem in the sense of Definition 1.4 by performing the limit passage $N\rightarrow \infty$ in (3.7)–(3.10). We start with the momentum equation. We consider a fixed $\psi _1\in C^1_c([0,T);V(\Omega ))$ in (3.7). Since $\widetilde{\rho v}^{N}$ is bounded in $L^{\infty }(0,T;L^{2}(\Omega ))$, which follows from (3.12)$_{1}$ and (3.16), we have

$$\begin{aligned} \mathbb {P}_{{{\,\textrm{div}\,}}}(\widetilde{\rho v}^N(t))\rightharpoonup \mathbb {P}_{{{\,\textrm{div}\,}}}(\rho v(t))\text { in }L^2(\Omega )\text { for a.e. }t\in (0,T), \end{aligned}$$

(3.60)

where the weak limit in (3.60) is identified by using (3.25). Fixing $\tau \in (0,T)$ such that (3.60) holds we take into consideration that $\partial ^-_{t,h}\rho ^N v^N=\partial _t\widetilde{\rho v}^N$ by (3.6)$_1$ and integrate by parts with respect to time in (3.7) to obtain

$$\begin{aligned} \begin{aligned}&\int _\Omega \widetilde{\rho v}^N(\tau )\psi _1(\tau )-\int _\Omega \widetilde{\rho v}^N(0)\psi _1(0)+\int _0^\tau \left( \int _{\Omega }- \widetilde{\rho v}^N\cdot \partial _t\psi _{1}\right. \\&\qquad \left. -\int _{\Omega }(\rho ^N_hv^N\otimes v^N)\cdot \nabla {\psi }_{1}-\int _{\Omega }v^N\otimes J^N\cdot \nabla \psi _{1}\right. \\&\qquad \left. +\int _{\Omega }\left( {{\,\textrm{div}\,}}(\xi (\phi ^N_h)\nabla M^N)-\frac{{\xi (\phi ^N_h)}}{\alpha ^{2}}(|M^N|^{2}M^N-M^N_h)\right) \nabla M^N\cdot \psi _{1}\right) \\&\quad =\int _0^\tau \left( -2\int _{\Omega }\nu (\phi ^{N}_{h})\mathbb {D} v^N\cdot \mathbb {D}\psi _{1}-\int _{\Omega }\nabla \mu ^N\phi ^N_h\cdot \psi _{1}\right) . \end{aligned} \end{aligned}$$

Thanks to (3.60) and the definition (1.4)–(1.5) of the Leray projector, we pass to the limit in the first term on the left hand side of the latter identity. By the definition of $\widetilde{\rho v}^N(0)$ we have, employing also (3.1)$_2$,

$$\begin{aligned} \begin{aligned} \widetilde{\rho v}^N(0)=\rho ^N(-h)v^N(-h)&\displaystyle =\tfrac{1}{2}\left( {{\widetilde{\rho }}}_1+{{\widetilde{\rho }}}_2+({{\widetilde{\rho }}}_2-{{\widetilde{\rho }}}_1)\phi ^N_0\right) v_0\\&\rightarrow \tfrac{1}{2}\left( \widetilde{\rho }_1+{{\widetilde{\rho }}}_2+({{\widetilde{\rho }}}_2-{{\widetilde{\rho }}}_1)\phi _0\right) v_0=\rho _0v_0\text { in }L^1(\Omega ), \end{aligned} \end{aligned}$$

which allows us to perform the passage in the second term. To pass to the limit in the third term we use (3.25) and the definition (1.4)–(1.5) of the Leray projector. We perform the limit passage in the fourth term with the help of (3.22) and (3.32). We recall that $\displaystyle J^N=-\frac{{{\widetilde{\rho }}}_2-{{\widetilde{\rho }}}_1}{2}\nabla \mu ^N$. Hence combining the convergences (3.17)$_6$ and (3.32) ensures the limit passage in the fifth term. For the limit passage in the last term on the left hand side we use (3.20) and (3.34). The limit passage on the right hand side is ensured by

$$\begin{aligned} \displaystyle {\nu (\phi ^{N}_{h})}\mathbb {D}(v^{N})\rightharpoonup {\nu (\phi )}\mathbb {D}(v)\,\,\text{ in }\quad L^{2}(Q_{T}), \end{aligned}$$

whose proof can be found in [36, eq. (5.56)] and (3.17)$_6$ combined with (3.18). We arrive at

$$\begin{aligned} \begin{aligned}&\int _\Omega \rho v(\tau )\psi _1(\tau )-\int _\Omega \rho v(0)\psi _1(0)+\int _0^\tau \left( \int _{\Omega }-\rho v\cdot \partial _t\psi _{1}\right. \\&\qquad \left. -\int _{\Omega }(\rho v\otimes v)\cdot \nabla {\psi }_{1} -\int _{\Omega }v\otimes J\cdot \nabla \psi _{1}\right. \\&\qquad \left. +\int _{\Omega }\left( {{\,\textrm{div}\,}}(\xi (\phi )\nabla M)-\frac{{\xi (\phi )}}{\alpha ^{2}}(|M|^{2}M-M)\right) \nabla M\cdot \psi _{1}\right) \\&\quad =\int _0^\tau \left( -2\int _{\Omega }\nu (\phi )\mathbb {D} v\cdot \mathbb {D}\psi _{1}-\int _{\Omega }\nabla \mu \phi \cdot \psi _{1}\right) . \end{aligned} \end{aligned}$$

(3.61)

Next we consider $t\in (0,T)$ and a sequence $\{\tau ^k\},$ s.t. $\tau ^k\rightarrow t$ and the latter identity holds for $\tau =\tau ^k$. Employing (3.56)$_1$ and the fact that all terms under the integration sign over the time interval are integrable with respect to time we conclude (1.13)$_1$ by the limit passage $k\rightarrow \infty$.

We note that the validity of identities (1.13)$_{2,3}$ (by the limit passage in (3.8) and (3.9)) can be proved by following line by line the arguments used to show [36, $(2.4)_{2}$ and $(2.4)_{3}$] in [36, Section 5.4, p. 31]. In order to verify that (1.13)$_4$ is fulfilled, we pass to the limit $N\rightarrow \infty$ in (3.10). In view of the convergences (3.17)$_{6}$, (3.18), (3.34) and (3.19) we conclude

$$\begin{aligned}&\mu ^N+\frac{\kappa }{2}(\phi ^N+\phi ^N_h)-H_0(\phi ^N,\phi ^N_h)\frac{|\nabla M^N|^{2}}{2}-\frac{H_0(\phi ^N,\phi ^N_h)}{4\alpha ^{2}}(|M^N|^{2}-1)^{2}\\&\quad \rightharpoonup \mu +\kappa \phi -\xi '(\phi )\frac{|\nabla M|^{2}}{2}-\frac{\xi '(\phi )}{4\alpha ^{2}}(|M|^{2}-1)^{2} \text { in }L^1(Q_T). \end{aligned}$$

Indeed, the passage to the limit in the first two terms is straightforward and the $L^1$ weak convergence of the remaining two terms is explained in detail in [36, Section 5.4, p. 32]. For the limit passage in the terms on the right hand side of (3.10) we use the convergence

$$\begin{aligned} -\eta \Delta \phi ^N+{{\widetilde{\Psi }}}'_0(\phi ^N)\rightharpoonup -\eta \Delta \phi +{{\widetilde{\Psi }}}'_0(\phi )\text { in }L^1(Q_T), \end{aligned}$$

which follows by (3.50) and (3.51). Thus we arrive at

$$\begin{aligned} \begin{aligned}&\int _0^T\int _\Omega \left( \mu +\kappa \phi -\xi '(\phi )\frac{|\nabla M|^{2}}{2}-\frac{\xi '(\phi )}{4\alpha ^{2}}(|M|^{2}-1)^{2}\right) \psi _4\\&\quad =\int _0^T\int _\Omega \left( -\eta \Delta \phi +{\widetilde{\Psi }}'_{0}(\phi )\right) \psi _4 \end{aligned} \end{aligned}$$

for all $\psi _4\in L^\infty (0,T;L^\infty (\Omega ))$. Hence it follows that identity (1.13)$_4$ is fulfilled.

3.3 The attainment of initial data $v_0, M_0, \phi _0$

In this section, we prove (1.14) with the help of (1.13), which we proved in the previous section. First we show the following identities

$$\begin{aligned} \begin{aligned} v(0)=&v_0{} & {} \text { a.e. in }\Omega ,\\ M(0)=&M_0{} & {} \text { a.e. in }\Omega ,\\ \phi (0)=&\phi _0{} & {} \text { a.e. in }\Omega . \end{aligned} \end{aligned}$$

(3.62)

Setting $\psi _3(t,x)=\theta (t)\vartheta (x)$ in (1.13)$_3$, where $\theta \in C^1_c([0,T))$ with $\theta (0)>0$ and $\vartheta \in C^{\infty }_{c}(\Omega )$ are arbitrary but fixed, we obtain using (3.56)$_{5}$

$$\begin{aligned} \int _\Omega \phi _{0}\theta (0)\vartheta =\lim _{t\rightarrow 0_+}\int _\Omega \phi (t)\theta (t)\vartheta =\int _\Omega \phi (0)\theta (0)\vartheta , \end{aligned}$$

which implies (3.62)$_3.$ Setting $\psi _1(t,x)=\theta (t)\omega (x)$ in (1.13)$_1$, where $\theta \in C^1_c([0,T))$ with $\theta (0)>0$ and $\omega \in V(\Omega )$ are arbitrary but fixed, yields

$$\begin{aligned} \int _\Omega \rho _0v_0\cdot \theta (0)\omega =\lim _{t\rightarrow 0_+}\int _\Omega \rho (t)v(t)\cdot \theta (t)\omega =\int _\Omega \rho (0)v(0)\cdot \theta (0)\omega , \end{aligned}$$

where the second equality follows by (3.56)$_{1}$ and (3.62)$_3$ implies $\rho (0)=\rho _0$. Setting in the latter identity $\omega =v_0-v(0)$, which is allowed due to the density of $V(\Omega )$ in $L^2_{{{\,\textrm{div}\,}}}(\Omega )$, implies (3.62)$_1$. We note that the fact that $\rho _0$ has a positive lower bound was also used. Finally, we repeat the above arguments to justify (3.62)$_2$.

With the help of (3.56) we will show that the energy inequality (1.16) holds for all $t\in [0,T]$. We start by considering an arbitrary $t\in [0,T]$ and a sequence $\{t^k\}$ such that $t^k\geqslant t$, $t^k\rightarrow t$ as $k\rightarrow \infty$ and

$$\begin{aligned} \begin{aligned} \rho (t^k)v(t^k)&\rightharpoonup \rho (t)v(t){} & {} \text{ in } L^2(\Omega ),\\ v(t^k)&\rightharpoonup v(t){} & {} \text{ in } L^2(\Omega ),\\ M(t^{k})&\rightharpoonup M(t){} & {} \text{ in } W^{1,2}(\Omega ),\\ M(t^{k})&\rightarrow M(t){} & {} \text{ in } L^{2}(\Omega ),\\ \phi (t^{k})&\rightharpoonup \phi (t){} & {} \text{ in } W^{1,2}(\Omega ),\\ \phi (t^{k})&\rightarrow \phi (t){} & {} \text{ in } L^{2}(\Omega )\text { and a.e. in }\Omega ,\\ \rho (t^k)&\rightarrow \rho (t){} & {} \text{ in } L^{2}(\Omega )\text { and a.e. in }\Omega \end{aligned} \end{aligned}$$

(3.63)

and (1.16) holds for each $t^k$. The convergence (3.63)$_{7}$ follows from (3.63)$_{6}$ by using the definition (1.9) of $\rho .$ The existence of such a sequence $\{t^k\}$ is ensured by (3.56). Because of the convexity of $|\cdot |^2$ (i.e. the inequality $|A|^{2}-|B|^{2}\geqslant 2B\cdot (A-B),$ for all $A,B\in \mathbb {R}^{m},$ $m\geqslant 1$) and convergences (3.63)$_{1,2,7}$ it follows that

$$\begin{aligned} \begin{aligned}&\liminf _{k\rightarrow \infty }\frac{1}{2}\int _\Omega \rho (t^k)|v(t^k)|^2 \\&\quad \geqslant \liminf _{k\rightarrow \infty }\int _\Omega \left( \frac{1}{2}\rho (t^k)|v(t)|^2+\rho (t^k)v(t)\cdot (v(t^k)-v(t))\right) \\&\quad =\frac{1}{2}\int _\Omega \rho (t)|v(t)|^2. \end{aligned} \end{aligned}$$

(3.64)

For the passage to the limit $k\rightarrow \infty$ in both terms we have used that $\rho (t^k)v(t)\rightarrow \rho (t)v(t)$ in $L^2(\Omega )$ (which follows by using Lebesgue’s dominated convergence theorem and the fact that $\rho$ is bounded) and also (3.63)$_{2}$ in the second term. Due to the weak lower semicontinutiy of convex functionals, the fact that ${\widetilde{\Psi }}\in C([-1,1])$ and (3.63)$_{5,6}$, we obtain

$$\begin{aligned} \begin{aligned}&\liminf _{k\rightarrow \infty }\int _\Omega \left( \frac{\eta }{2}|\nabla \phi (t^k)|^2+{\widetilde{\Psi }}(\phi (t^k)) \right) \geqslant \int _\Omega \left( \frac{\eta }{2}|\nabla \phi (t)|^2+{\widetilde{\Psi }}(\phi (t)\right) . \end{aligned} \end{aligned}$$

(3.65)

Moreover, we obtain

$$\begin{aligned} \begin{aligned}&\liminf _{k\rightarrow \infty } \int _{\Omega }\left( \xi (\phi (t^k))|\nabla M(t^k)|^2+\frac{\xi (\phi (t^k))}{\alpha ^2}(|M(t^k)|^2-1)^2\right) \\&\quad \geqslant \int _{\Omega }\left( \xi (\phi (t))|\nabla M(t)|^2+\frac{\xi (\phi (t))}{\alpha ^2}(|M(t)|^2-1)^2\right) , \end{aligned} \end{aligned}$$

(3.66)

by arguing as in [36, (5.77)]. Altogether, (3.64), (3.65), (3.66) and the absolute continuity of the map

$$\begin{aligned} \begin{aligned}&t\mapsto \int _0^t\biggl (\Vert \sqrt{2\nu }\mathbb {D} v\Vert ^2_{L^2(\Omega )}+\Vert \nabla \mu \Vert ^2_{L^2(\Omega )}\\&\quad +\left\| {{\,\textrm{div}\,}}(\xi (\phi )\nabla M)-\frac{\xi (\phi )}{\alpha ^2}M(|M|^2-1)\right\| ^2_{L^2(\Omega )}\biggr ) \end{aligned} \end{aligned}$$

(3.67)

imply that (1.16) holds for all $t\in [0,T]$. Hence in particular it follows that

$$\begin{aligned} \limsup _{t\rightarrow 0_+}E_{tot}(v(t),M(t),\phi (t))\leqslant E_{tot}(v_0,M_0,\phi _0). \end{aligned}$$

(3.68)

Employing again (3.56) along with (3.62) (similarly as we have obtained (3.64)–(3.66)) we deduce

$$\begin{aligned} \liminf _{t\rightarrow 0_+}E_{tot}(v(t),M(t),\phi (t))\geqslant E_{tot}(v(0),M(0),\phi (0))=E_{tot}(v_0,M_0,\phi _0), \end{aligned}$$

which along with (3.68) infers

$$\begin{aligned} \lim _{t\rightarrow 0_+}E_{tot}(v(t),M(t),\phi (t))= E_{tot}(v_0,M_0,\phi _0). \end{aligned}$$

(3.69)

Taking into account the definition of $E_{tot},$ employing the inequalities $|A|^2-|B|^2\geqslant 2B\cdot (A-B)+2|A-B|^2$ (which follows from the strong convexity of $|\cdot |^2$) and $|A|^4-|B|^4\geqslant 4|B|^{2}B\cdot (A-B)$ (which follows from the convexity of $|\cdot |^4$) for all $A,B\in \mathbb {R}^m,$ one obtains the following for each $t\in (0,T)$

$$\begin{aligned} \begin{aligned}&E_{tot}(v(t),M(t),\phi (t))-E_{tot}(v_0,M_0,\phi _0)\\&\quad \geqslant \frac{1}{2}\int _\Omega (\rho (t)-\rho _0)|v_0|^2+ \int _\Omega \rho (t)v_0\cdot (v(t)-v_0)+\int _\Omega \rho (t)|v(t)-v_0|^2\\&\qquad +\int _\Omega \left( \xi (\phi (t))-\xi (\phi _0)\right) |\nabla M_0|^2+\int _\Omega 2\xi (\phi (t))\nabla M_0\cdot (\nabla M(t)-\nabla M_0)\\&\qquad +2\int _\Omega \xi (\phi (t))|\nabla M(t)-\nabla M_0|^2+\frac{1}{4\alpha ^2}\int _\Omega \left( \xi (\phi (t))-\xi (\phi _0)\right) |M_0|^4 \\&\qquad +\frac{1}{\alpha ^2}\int _\Omega \xi (\phi (t))|M_0|^2M_0\cdot \left( M(t)-M_0\right) \\&\qquad -\frac{1}{2\alpha ^2}\int _\Omega \left( \xi (\phi (t))|M(t)|^2-\xi (\phi _0)|M_0|^2\right) +\frac{1}{4\alpha ^2}\int _\Omega \left( \xi (\phi (t))-\xi (\phi _0)\right) \\&\qquad {+}\!\eta \int _\Omega \nabla \phi _0\cdot \left( \nabla \phi (t){-}\nabla \phi _0\right) \!{+}\!\eta \int _\Omega |\nabla \phi (t){-}\nabla \phi _0|^2\!{+}\!\int _{\Omega }\left( {\widetilde{\Psi }}_{0}(\phi )-{\widetilde{\Psi }}_{0}(\phi _{0})\right) \\&\qquad -\frac{\kappa }{2}\int _\Omega \left( \phi ^2(t)-\phi _0^2\right) =\sum _{m=1}^{14}I_m(t). \end{aligned}\nonumber \\ \end{aligned}$$

(3.70)

The task now is to prove (1.14) by taking the limsup $t\rightarrow 0_+$ on both sides of the inequality (3.70). To this end we consider an arbitrary sequence $\{t^k\}$ such that $t^k\rightarrow 0_+$ as $k\rightarrow \infty$. The sequence $\{t^k\}$ has a subsequence $\{t^{k'}\}$ such that the following hold

$$\begin{aligned} \phi (t^{k'})\rightarrow \phi _0,\ \text { a.e. in }\Omega \text { as }k'\rightarrow \infty \end{aligned}$$

(3.71)

by (3.56)$_6$ and (3.62)$_3$. Accordingly, we have

$$\begin{aligned} \rho (t^{k'})\rightarrow \rho _0,\ \text { a.e. in }\Omega \text { as }k'\rightarrow \infty \end{aligned}$$

(3.72)

by (1.9). For the proof of

$$\begin{aligned} \lim _{k^{'}\rightarrow \infty }I_m(t^{k'})=0\text { for }m=4,5,7,8,9,10,11 \end{aligned}$$

(3.73)

we refer to [36, (5.82)-(5.85)]. Next we deal with $I_1$, $I_2$, $I_{13}$ and $I_{14}$. Convergence (3.72) and the fact that $\rho$ is a bounded function imply

$$\begin{aligned} \lim _{k'\rightarrow \infty }I_1(t^{k'})=0 \end{aligned}$$

(3.74)

by the Lebesgue dominated convergence theorem. Moreover, we have that $\rho (t^{k'})v_0\rightarrow \rho _0v_0$ in $L^2(\Omega )$, which along with (3.56)$_{2}$ and (3.62)$_1$ yields

$$\begin{aligned} \lim _{k'\rightarrow \infty }I_2(t^{k'})=0. \end{aligned}$$

(3.75)

Since ${\widetilde{\Psi }}_{0}\in C([-1,1]),$ the following

$$\begin{aligned} \lim _{k'\rightarrow \infty }I_{13}(t^{k'})=0 \end{aligned}$$

(3.76)

is obtained as an immediate consequence of (3.71) and the Lebesgue dominated convergence theorem.

Finally, by (3.56)$_{6}$ and (3.62)$_3$ we obtain

$$\begin{aligned} \lim _{k'\rightarrow \infty }I_{14}(t^{k'})=0. \end{aligned}$$

(3.77)

Hence

$$\begin{aligned} \begin{aligned}&\limsup _{k '\rightarrow \infty } \left( \underline{\rho }\Vert v(t^{k'})-v_0\Vert ^2_{L^2(\Omega )}+c_1\Vert \nabla M(t^{k'})-\nabla M_0\Vert ^2_{L^2(\Omega )}\right. \\&\quad \left. +\eta \Vert \nabla \phi (t^{k'})-\nabla \phi _0\Vert ^2_{L^2(\Omega )}\right) \leqslant 0 \end{aligned} \end{aligned}$$

(3.78)

follows from (3.70) by (3.73)–(3.77) provided that we apply (1.6)$_2$ and take into consideration that there is a positive lower bound on $\rho$, which we denote by $\underline{\rho }.$

The inequality (3.78) along with (3.56)$_{4,6}$ infer

$$\begin{aligned} \begin{aligned}&\lim _{k '\rightarrow \infty }\left( \Vert v(t^{k'})-v_0\Vert _{L^2(\Omega )}+\Vert M(t^{k'})- M_0\Vert _{W^{1,2}(\Omega )}\right. \\&\quad \left. +\Vert \phi (t^{k'})-\phi _0\Vert _{W^{1,2}(\Omega )}\right) =0. \end{aligned} \end{aligned}$$

(3.79)

Since $\{t^{k}\}$ is an arbitrarily sequence possessing a subsequence satisfying (3.79), one concludes the proof of (1.14).

3.4 Attainment of the boundary condition and some regularity results for M in Lebesgue spaces

In this section we discuss the proofs of the items (ii) and (iii) of Theorem 1.5. For the proof of the item (ii) we refer the readers to [36, Section 6.1]. The item (iii) was formally commented in [36, Section 6.2] but one needs to suitably regularize the magnetization equation to make the arguments concrete. Here we provide the details for the proof of item (iii).

In the direction of proving item (iii) of Theorem 1.5, we first show that for given v and $\phi$ in the functional settings (1.12)$_{1,3}-$(1.15)$_{2}$ there is a unique M satisfying (1.12)$_{2}-$(1.15)$_{1}$ and solving the weak formulation (1.13)$_{2}$ of the magnetization equation. Since $\partial _{t}M\in L^{2}(0,T;L^{\frac{3}{2}}(\Omega )),$ equation (1.13)$_{2}$ can be rewritten as:

$$\begin{aligned} \begin{aligned} \int _0^t\int _\Omega \biggl (\partial _{t}M+(v\cdot \nabla ) M\biggl )\cdot \psi _{2}=&-\int _0^t\int _\Omega \xi (\phi )\nabla M\cdot \nabla \psi _{2}\\&-\int _0^t\int _\Omega \frac{1}{\alpha ^2}\bigl (\xi (\phi )(|M|^2-1)M\bigr )\cdot \psi _{2} \end{aligned} \end{aligned}$$

(3.80)

for $t\in (0,T)$ and $\psi _{2}\in C^{1}_{c}(0,T;W^{1,2}(\Omega ))$ or equivalently

$$\begin{aligned} \begin{aligned} \displaystyle \int _\Omega \partial _{t}M\cdot \psi _{2}+\int _\Omega (v\cdot \nabla ) M\cdot \psi _{2}=&-\int _\Omega \xi (\phi )\nabla M\cdot \nabla \psi _{2}\\&-\int _\Omega \frac{1}{\alpha ^2}\bigl (\xi (\phi )(|M|^2-1)M\bigr )\cdot \psi _{2} \end{aligned} \end{aligned}$$

(3.81)

for a.e. $t\in (0,T)$ and $\psi _{2}\in W^{1,2}(\Omega ).$

Let $M_{1}$ and $M_{2}$ belong to (1.12)$_{2}-$(1.15)$_{1}$ and solve (3.80) with v and $\phi$ in the framework (1.12)$_{1,3}-$ (1.15)$_{2}.$ One can now take the difference of the equations solved by $M_{1}$ and $M_{2}$ and consider $(M_{1}-M_{2})$ as a test function, which is possible since $C^{1}_{c}(0,T;W^{1,2}(\Omega ))$ is dense in $L^{2}(0,T;W^{1,2}(\Omega )).$ Consequently using the incompressibility of v and the inequality $\bigl (|M_1|^2M_1-|M_2|^2M_2\bigr )\cdot (M_1-M_2)\geqslant 0$ (since the map $\alpha \mapsto |\alpha |^{2}\alpha$ is monotone) one furnishes

$$\begin{aligned} \frac{1}{2}\Vert (M_{1}-M_{2})(t)\Vert ^{2}_{L^{2}(\Omega )}\leqslant C\int _{0}^{t}\Vert M_{1}-M_{2}\Vert ^{2}_{L^{2}(\Omega )}, \end{aligned}$$

for a.e. $t\in (0,T).$ Hence by the Grönwall inequality one at once renders that $M_{1}=M_{2}$ a.e. in $Q_{T}.$

Now we plan to use test functions of the form $|M|^{r-2}M$ with $r>2$ in (3.81). But due to the lack of regularity (particularly one needs for a.e. $t\in (0,T),$ $M\in L^{r-1}(\Omega )$ for arbitrary $r>2$) this does not qualify as a test function. Instead we consider a regularized magnetization equation, i.e. we first take a sequence $\{\phi ^{m}\}_{m}$ in $L^{2}(0,T;C^{\infty }(\overline{\Omega }))$ such that

$$\begin{aligned} \phi ^{m}\rightarrow \phi \,\,\text{ in }\,\, L^{2}(Q_{T}) \end{aligned}$$

(such a sequence can easily be constructed by a suitable argument involving cut-off and convolution by mollifiers). Now let $M^{m}$ be the weak solution to (3.80) or (3.81) corresponding to $\phi ^{m}$ with boundary condition $\partial _{n}M^{m}\mid _{\Sigma _{T}}=0$ and initial condition $M^{m}(\cdot ,0)=M_{0}\in W^{1,2}(\Omega ).$ Our idea is to consider $|M^{m}|^{r-2}M^{m}$ as a test function in the equation solved by $(\phi ^{m},M^{m})$ thereby proving an uniform estimate of $M^{m}$ in $L^{r}(\Omega )$ and next pass $m\rightarrow \infty$ to construct a weak solution M corresponding to $\phi$ for (3.80) or equivalently (3.81) which also solves the desired $L^{r}(\Omega )$ estimate. Of course, because of the uniqueness of the solution of the magnetization equation corresponding to the fixed pair $(v,\phi )$ and the initial data $M_{0},$ which we have already proved, this process will give the same M solving (1.13)$_{2}.$

With the help of a time discretization scheme one can prove the existence of a weak solution $M^{m}\in L^{\infty }(0,T;W^{1,2}(\Omega ))\cap W^{1,2}(0,T;L^{\frac{3}{2}}(\Omega ))$ of (3.80) or equivalently (3.81) corresponding to a vector field v (satisfying (1.12)$_{1}$) and $\phi ^{m}.$ Moreover we notice that, in a strong form, this $M^{m}$ solves

$$\begin{aligned} \begin{aligned} \Delta M^{m}=&\frac{1}{\xi (\phi ^{m})}\Bigl (\partial _{t}M^{m}+(v\cdot \nabla )M^{m}-\xi '(\phi ^{m})\nabla M^{m}\cdot \nabla \phi ^{m}\\&+\frac{\xi (\phi ^{m})}{\alpha ^{2}}\bigl (|M^{m}|^{2}-1\bigr )M^{m}\Bigr ){} & {} \text { in }\Omega ,\\ \partial _{n}M^{m}=&0{} & {} \text { on }\partial \Omega .\\ \end{aligned} \end{aligned}$$

(3.82)

In view of the fact that $M^{m}\in L^{\infty }(0,T;W^{1,2}(\Omega ))\cap W^{1,2}(0,T;L^{\frac{3}{2}}(\Omega ))$ the right hand of (3.82)$_{1}$ can be estimated in $L^{2}(0,T;L^{\frac{3}{2}}(\Omega ))$ and hence by standard elliptic regularity $M^{m}\in L^{2}(0,T;W^{2,\frac{3}{2}}(\Omega ))\hookrightarrow L^{2}(0,T;L^{r}(\Omega ))$ for any $0<r<\infty .$ Hence for a.e. $t\in (0,T),$ $|M^{m}|^{r-1}M^{m}(t),$ $r>2$ can be used as a test function in (3.81). Consequently

$$\begin{aligned} \begin{aligned}&\frac{1}{r}\partial _{t}\Vert M^{m}\Vert ^{r}_{L^{r}(\Omega )}+\int _{\Omega }(v\cdot \nabla )M^{m}|M^{m}|^{r-2}M^{m}\\&\quad +\int _{\Omega } \xi (\phi ^{m})(r-1)|M^{m}|^{r-2}|\nabla M^{m}|^{2}+\frac{1}{\alpha ^{2}}\int _{\Omega }\xi (\phi ^{m})|M^{m}|^{r+2}\\&\quad -\frac{1}{\alpha ^{2}}\int _{\Omega }\xi (\phi ^{m})|M^{m}|^{r}=0. \end{aligned} \end{aligned}$$

(3.83)

Once again integrating by parts the second term and using that $\text{ div }\,v=0$ on $\Omega$ one concludes from (3.83) that:

$$\begin{aligned} \begin{array}{l} \displaystyle \partial _{t}\Vert M^{m}\Vert ^{r}_{L^{r}(\Omega )}\leqslant \frac{c_{2}r}{\alpha ^{2}}\Vert M^{m}\Vert ^{r}_{L^{r}(\Omega )}, \end{array} \end{aligned}$$

(3.84)

where $c_{2}>0$ is the constant appearing in the assumption (1.6). Now if one assumes $M_{0}\in W^{1,2}(\Omega )\cap L^{r}(\Omega ),$ $r>6,$ using Gronwall’s inequality one has the following from (3.84):

$$\begin{aligned} \begin{array}{l} \displaystyle \Vert M^{m}(t)\Vert _{L^{r}(\Omega )}\leqslant \Vert M_{0}\Vert _{L^{r}(\Omega )}e^{\frac{c_{2}}{\alpha ^{2}}t}\quad \text{ for } \text{ all }\quad t\in [0,T]. \end{array} \end{aligned}$$

(3.85)

Additionally if $M_{0}\in L^{\infty }(\Omega ),$ one can take the limit $r\rightarrow \infty$ in (3.85) to conclude that:

$$\begin{aligned} \begin{array}{l} \displaystyle \Vert M^{m}(t)\Vert _{L^{\infty }(\Omega )}\leqslant \Vert M_{0}\Vert _{L^{\infty }(\Omega )}e^{\frac{c_{2}}{\alpha ^{2}}t}\quad \text{ for } \text{ all }\quad t\in [0,T]. \end{array} \end{aligned}$$

(3.86)

Now we let $m\rightarrow \infty$ in the equation solved by $(\phi ^{m},M^{m}),$ i.e. (3.80) with $(\phi ,M)$ replaced by $(\phi ^{m},M^{m}).$ The limit passage in the equation is obtained in a standard way (roughly it consists in showing the weak compactness of $M^{m}$ in $L^{2}(0,T;W^{1,2}(\Omega ))$ and next using Aubin-Lions to achieve the strong compactness in $L^{2}(0,T;L^{4}(\Omega ))$). Finally in view of the estimates (3.85) and (3.86), which are independent of m,  one concludes (1.20) and (1.21).

3.5 Summary of the proof of Theorem 1.5

For the sake of the readers we summarize the proof of Theorem 1.5 with exact references to the sections.

• For the obtainment of the regularities (1.12) with the exception of $M\in W^{1,2}(0,T;L^{\frac{3}{2}}(\Omega )),$ $\phi \in L^{2}(0,T;W^{2,1}(\Omega ))$ and $\Psi '(\phi )\in L^{1}(Q_{T}),$ we refer the readers to (3.17) and (3.56). One can obtain the $W^{1,2}(0,T;L^{\frac{3}{2}}(\Omega ))$ regularity of M simply by estimating $\partial _{t}M\in L^{2}(0,T;L^{\frac{3}{2}}(\Omega ))$ by using (1.1)$_{3}$ and the available regularities for v and M. More precisely $(v\cdot \nabla )M$ can be estimated in $L^{2}(0,T;L^{\frac{3}{2}}(\Omega ))$ as in (3.44) and the boundedness of ${{\,\textrm{div}\,}}(\xi (\phi )\nabla M)-\frac{\xi (\phi )}{\alpha ^{2}}(|M|^{2}-1)M$ in $L^{2}(Q_{T})$ follows from (1.16). The additional $p-$regularities (1.15) of M,  $\phi$ and $\Psi '(\phi )$ can be found in Sect. 3.1.3 and they are of course stronger than $\phi \in L^{2}(0,T;W^{2,1}(\Omega ))$ and $\Psi '(\phi )\in L^{1}(Q_{T})$ (stated as a part of (1.12)).

• The weak formulation (1.13) solved by $(v,M,\phi ,\mu )$ is proved in Sect. 3.2.

• The energy estimate (1.16) is obtained in Sect. 3.1.4.

• The attainment of the initial data in the sense of (1.14) is obtained in Sect. 3.3.

• The items (ii) and (iii) of Theorem 1.5 corresponding to the attainment of boundary condition for M and some regularity in Lebesgue spaces are proved in Sect. 3.4.

In view of the above items we finally conclude the proof of Theorem 1.5. $\square$