1 sec:Introduction↩︎

The Marangoni effect, named after Carlo Marangoni, describes the fluid motion generated along an interface due to gradients in surface tension [1]. This inhomogeneity of surface tension typically arises from spatial variations in temperature (thermocapillarity) or solute concentration (solutocapillarity) [2], [3]. When thermal effects are considered, for most fluids, the surface tension decreases with increasing temperature [4], and consequently, the resulting tangential stress induces a convective flow from the warmer regions (lower surface tension) to the cooler regions (higher surface tension). This ubiquitous phenomenon was recognized as a primary driving mechanism that contributes to fluid motion alongside buoyancy effects in the classic Bénard experiments [5]. Marangoni convection can influence the spreading behavior near interfaces and lead to the appearance of dynamic organized patterns. It has become a fundamental subject of research, with diverse applications in crystal growth/welding, electron beam melting, nanotechnology, and biology (see, e.g., [6]–[9]).

In this study, we consider the evolution of a non-isothermal fluid mixture with two incompressible immiscible viscous Newtonian fluids under the influence of thermally induced surface tension gradients. In the sharp interface formulation (see, e.g., [10], [11]), the transition region between immiscible constituents is idealized as a zero-thickness hypersurface separating two distinct bulk domains. Although this treatment provides a clear mathematical representation of phase boundaries, it becomes challenging during topological changes such as droplet pinch-off and coalescence, where vanishing length scales can generate singular interfacial dynamics, see [12], [13]. This difficulty has motivated alternative interfacial descriptions that can regularize the sharp-boundary idealization and provide a more flexible framework for describing interfacial deformation and topological transitions. Phase-field models have emerged as an efficient method for studying interfacial dynamics and phase transitions in multi-component fluids, see [14], [15] and the references therein. In this approach, the sharp interface is regularized as a thin transition layer where physical variables vary continuously but steeply, i.e., the so-called diffuse interface. The diffuse interface regularization avoids explicit interface tracking while naturally incorporates thermodynamic driving forces, facilitating both the mathematical analysis and numerical implementation for multi-phase flows [14], [16], [17].

Let \(\Omega \subset \mathbb{R}^d\) (\(d = 2\) or \(3\)) be a bounded domain with a \(C^3\)-boundary \(\partial \Omega\). We consider the following hydrodynamic phase-field system: \[\begin{align} &\partial_t(\rho(\phi)\boldsymbol{u}) + \mathrm{div}(\boldsymbol{u}\otimes(\rho(\phi)\boldsymbol{u}+\mathbf{J})) - \mathrm{div}\,(2 \nu(\phi,\theta) D \boldsymbol{u}) + \nabla p \notag \\ &\qquad = - \mathrm{div}\, \boldsymbol{\sigma} + \boldsymbol{f}_{\mathrm{b}}(\phi, \theta), \tag{1}\\ &\mathrm{div}\, \boldsymbol{u}= 0, \tag{2}\\ &\partial_t \phi +\boldsymbol{u}\cdot \nabla \phi = \mathrm{div} (m(\phi,\theta) \nabla \mu), \tag{3}\\ &\mu = -\Delta \phi + W^{\prime}(\phi), \tag{4}\\ &\partial_t \theta +\boldsymbol{u}\cdot \nabla \theta= \mathrm{div}\,(\kappa(\phi,\theta) \nabla \theta), \tag{5} \end{align}\] in \(\Omega \times (0,\infty)\), subject to the boundary and initial conditions: \[\begin{align} {4} &\boldsymbol{u}=\mathbf{0}, \quad &&\partial_{\mathbf{n}} \phi = \partial_{\mathbf{n}} \mu =0, \quad &&\theta=\theta_{\mathrm{b}}, \quad &&\text{on} ~ \partial \Omega \times (0,\infty), \tag{6}\\ &\boldsymbol{u}|_{t=0}=\boldsymbol{u}_0, \quad &&\phi|_{t=0}=\phi_0, \quad &&\theta|_{t=0} = \theta_0, \quad &&\text{in} ~ \Omega. \tag{7} \end{align}\] In 6 , \(\mathbf{n}=\mathbf{n}(x)\) denotes the outward unit normal vector on \(\partial\Omega\) and \(\partial_{\mathbf{n}}\) denotes the outward normal derivative on the boundary, \(\theta_{\mathrm{b}}=\theta_{\mathrm{b}}(x)\) is a given function on \(\partial\Omega\). The state variables of the system 1 –5 are denoted by \((\boldsymbol{u}, p, \phi, \mu, \theta)\). More precisely, \(\boldsymbol{u}:\Omega\times(0,\infty)\to\mathbb{R}^d\) denotes the volume-averaged velocity of the binary fluid mixture and \(p:\Omega\times(0,\infty)\to\mathbb{R}\) the pressure, which are governed by a modified Navier–Stokes system 1 –2 ; \(\phi:\Omega\times(0,\infty)\to [-1,1]\) is the phase-field variable and \(\mu:\Omega\times(0,\infty)\to\mathbb{R}\) the chemical potential, together forming the convective Cahn–Hilliard equation 3 –4 ; and \(\theta:\Omega\times(0,\infty)\to\mathbb{R}\) denotes the relative temperature (with respect to a given constant ambient temperature), which satisfies the convective heat transport equation 5 .

The coupled system 1 –7 presents a nontrivial coupling between fluid dynamics, phase separation process, and thermal convection, leading to complex interfacial phenomena [3], [18]. It not only generalizes the phase-field model for thermo-induced Marangoni effects derived in [19], [20] via the energetic variational approach, but also incorporates the thermodynamically consistent framework proposed in [21] for the general scenario with unmatched densities, based on the utilization of a volume-averaged velocity \(\boldsymbol{u}\) that keeps the binary fluid mixture incompressible (cf. 2 ). The phase-field variable \(\phi\) serves as a conserved order parameter that represents the difference in volume fractions of the fluid mixture, such that \(\{\phi=-1\}\) represents fluid \(1\) and \(\{\phi=1\}\) represents fluid \(2\). Then the average density \(\rho\) is assumed to be the typical linear form \[\begin{align} \rho(\phi)=\frac{\rho_2-\rho_1}{2}\phi+\frac{\rho_1+\rho_2}{2}, \label{density} \end{align}\tag{8}\] where the positive constants \(\rho_1\), \(\rho_2\) denote the homogeneous positive density of the unmixed components of the fluid. For simplicity, the density variation with respect to the temperature is explicitly considered in the buoyancy force \(\boldsymbol{f}_{\mathrm{b}}\) and assumed to satisfy the linearized thermal expansion approximation (see e.g., [20]), that is, \[\begin{align} \boldsymbol{f}_{\mathrm{b}}(\phi, \theta) = -\rho(\phi)(1-\alpha\theta)g \mathbf{e}_d. \label{bouy} \end{align}\tag{9}\] Here, \(\alpha\) is the coefficient of thermal expansion, \(g\) is the gravitational acceleration, and \(\mathbf{e}_d\) denotes the unit upward vector, i.e., \(\mathbf{e}_2=(0,1)^T\), \(\mathbf{e}_3=(0,0,1)^T\). In the equation of momentum balance 1 , the matrix-valued function \(D \boldsymbol{u}:= \frac{1}{2}(\nabla \boldsymbol{u}+ \nabla \boldsymbol{u}^T)\) denotes the symmetric gradient of \(\boldsymbol{u}\), while the relative flux related to the diffusion of the fluid components is given by \[\begin{align} \mathbf{J}=-\frac{\rho_2-\rho_1}{2}m(\phi,\theta)\nabla \mu. \label{Jflux} \end{align}\tag{10}\] The Cauchy stress tensor \(\boldsymbol{\sigma}\) is defined as \[\begin{align} \label{Cauchy95stress95tensor} \boldsymbol{\sigma} = \lambda(\theta) ( \nabla \phi \otimes \nabla \phi ) + \lambda(\theta) \left( \frac{1}{2}|\nabla \phi|^2 + W(\phi)\right) \mathbb{I}_d, \end{align}\tag{11}\] where \(\mathbb{I}_d\) denotes the \(d\)-dimensional unit matrix, and the temperature-dependent surface tension coefficient \(\lambda\) follows the empirical Eötvös law: \[\begin{align} \lambda(\theta) = \lambda_0(a - b \theta), \label{Eotvos} \end{align}\tag{12}\] with constant coefficients \(\lambda_0,a,b>0\). Since the surface tension depends on temperature, its spatial variation along the interface induces tangential stresses that drive interfacial fluid motion, i.e., the Marangoni convection. The mixing energy of binary fluids is given by the classical Ginzburg–Landau type free energy (cf. [14], [22]) \[\begin{align} \label{mix} E_{\mathrm{mix}}(\phi)= \int_\Omega \frac{1}{2}|\nabla \phi|^2+ W(\phi)\,\mathrm{d}x, \end{align}\tag{13}\] in which the gradient term contributes to the free-energy excess of the interfacial region and the nonlinear function \(W=W(\phi)\) denotes the bulk energy density with a double well structure. Here, we set the width of the diffuse interface to \(1\) for simplicity, since we do not consider the sharp interface limit in this study. A physically relevant example of \(W\) is the logarithmic (Flory–Huggins) potential [22]: \[\begin{align} \label{Wphi} W(\phi) = \frac{A}{2} \left[ (1+\phi)\ln(1+\phi) + (1-\phi)\ln(1-\phi) \right] - \frac{A_c}{2} \phi^2, \quad \phi \in (-1,1), \end{align}\tag{14}\] with constant coefficients \(A, A_c\) satisfying \(0 < A < A_c\). The chemical potential \(\mu\) in the Cahn–Hilliard equation 3 –4 is defined as the variational derivative of the mixing energy \(E_{\mathrm{mix}}\) with respect to \(\phi\) (subject to the homogeneous Neumann boundary condition \(\partial_\mathbf{n}\phi=0\)). Its gradient provides a driving force for the diffusion mechanism [16], [23]. The heat equation 5 gives a simplification of the thermal energy equation. In general, temperature dependent extensions of the Cahn–Hilliard equation for phase separation processes are rather involved, since different transport properties of the temperature and relations of the free energy can lead to different dynamical models, see [24]–[27] and the references therein.

In this study, we allow physical coefficients such as fluid viscosity \(\nu\), diffusion mobility \(m\) and thermal diffusivity \(\kappa\) to depend both on the order parameter \(\phi\) and temperature \(\theta\). The temperature dependence of the fluid viscosity and thermal diffusivity is physically important in the investigation of the detailed motion in certain non-isothermal flows, see, for instance, [28] and the references therein. On the other hand, for multi-phase fluids, the dependence of structural coefficients on the mixture composition is physically meaningful and should be taken into account in order to accurately model the complex fluid dynamics, since different components may have distinct physical properties, see, e.g., [19], [29].

Several simplified versions of the coupled system 1 –7 have been studied in the literature. For instance, assuming matched densities, constant mobility, and only temperature dependence of the viscosity and thermal diffusivity, in [30] the author proved the existence and uniqueness of global weak/strong solutions to the initial boundary value problem in two dimensions. When the Marangoni effect is neglected, still assuming matched densities and constant coefficients, problem 1 –7 reduces to the so-called Cahn–Hilliard–Boussinesq system, which has been extensively analyzed. We refer to [31]–[34] for results on well-posedness, regularity, and long-time behavior of solutions. Research on related models such as the Cahn–Hilliard–Oberbeck–Boussinesq system can be found in [35], where weak and very weak solutions in two dimensions were investigated. When the Cahn–Hilliard equation is replaced by a second order Allen–Cahn equation (with a regular potential like \(W(\phi)=\frac{1}{4}(\phi^{2}-1)^2\)), existence of global weak solutions, existence and uniqueness of strong solutions and long-time behavior of globally bounded solutions have been analyzed under suitable assumptions on the structural data, see [36]–[39] for detailed discussions. In addition, the local-in-time existence of weak solutions was established in [40] for the Allen–Cahn variant with a temperature-dependent interfacial thickness. Recently, the authors in [41] investigated a non-isothermal Navier–Stokes–Allen–Cahn system for incompressible two-phase flows of equal densities. Despite the high nonlinearity in the thermodynamically consistent heat equation, they established local well-posedness of the initial-boundary value problem and the existence of global “entropic” weak solutions. In addition, they analyzed the sharp interface limit, showing that an entropic weak solution to the phase-field model converges to a distributional (or BV) solution to a non-isothermal Navier–Stokes/mean curvature flow.

All of the aforementioned works assumed matched densities for the fluid components and sometimes adopted the Boussinesq approximation for the thermally induced variation of the density. On the other hand, when the thermal effect is neglected, problem 1 –7 reduces to the Abels–Garcke–Grün system (AGG for short, see [21]), which is a thermodynamically consistent and frame-indifferent diffuse interface model for isothermal viscous incompressible two-phase flows with unmatched densities. We refer to [42] for the existence of global weak solutions in two and three dimensions for a non-degenerate mobility \(m=m(\phi)\). The corresponding result in the case of degenerate mobility was established in [43]. Concerning local strong well-posedness in three dimensions, we refer to [44], see also [45] for the case with a regular potential. When the spatial dimension is two, the existence of strong solutions locally in time for bounded domains and globally in time for periodic boundary conditions was obtained in [46]. Finally, we refer to [47] for global regularity and asymptotic stabilization of global weak solutions in three dimensions as well as the existence of global strong solutions in two-dimensional bounded domains. Extensions to multi-component systems can be found in the recent work [48].

To the best of our knowledge, there is no analytical result on the hydrodynamic system 1 –7 , which extends the isothermal AGG model by coupling with the energy transport equation 5 and incorporating thermal effects on both density (buoyancy forces) and surface tension (Marangoni effects). In this study, we establish the following results for problem 1 –7 in the general scenario with unmatched densities and variable viscosity, mobility, as well as thermal diffusivity:

• the existence of global weak solutions in both two and three dimensions, which are uniformly bounded in time (see Theorem 1);

• the uniqueness of global weak solutions in two dimensions, under suitable assumptions on the structural data (see Theorem 2).

The problem 1 –7 satisfies two fundamental properties, i.e., mass conservation and energy balance, which serve as the basis for the subsequent analysis. Integrating 3 over \(\Omega\), using the boundary conditions 6 and the incompressibility condition 2 , after integration by parts, we obtain \[\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} \int_\Omega\phi \,\mathrm{d}x =0,\quad \forall\, t>0,\label{mass} \end{align}\tag{15}\] which implies that the mass of the binary fluids is conserved for all time. Next, for sufficiently smooth solutions, multiplying 1 by \(\boldsymbol{u}\), 3 by \(\lambda_0 a \mu\) and 5 by \(\theta\), integrating over \(\Omega\) (assuming for simplicity here \(\theta_\mathrm{b}=0\)), and adding the resultants together, we formally arrive at the following energy identity \[\begin{align} &\frac{\mathrm{d}}{\mathrm{d} t} \mathcal{E}_{\mathrm{tot}} (\boldsymbol{u}, \phi, \theta) + \int_\Omega \big( 2\nu(\phi, \theta)|D\boldsymbol{u}|^2 + \lambda_0 am(\phi,\theta)|\nabla \mu|^2 + \kappa(\phi, \theta)|\nabla \theta|^2\big)\,\mathrm{d}x \notag \\ &\quad = \int_\Omega \mathrm{div}\big[\lambda_0 b\theta(\nabla \phi\otimes \nabla \phi)\big]\cdot \boldsymbol{u}\,\mathrm{d}x + \int_\Omega \boldsymbol{f}_{\mathrm{b}}(\phi, \theta)\cdot \boldsymbol{u}\,\mathrm{d}x,\quad \forall\, t>0, \label{energy-a} \end{align}\tag{16}\] where the total energy \(\mathcal{E}_{\mathrm{tot}}\) is given by \[\begin{align} \label{energy-b} \mathcal{E}_{\mathrm{tot}} (\boldsymbol{u}, \phi, \theta) =\int_\Omega \left[\frac{1}{2}\rho(\phi)|\boldsymbol{u}|^2 + \lambda_0 a\left(\frac{1}{2}|\nabla \phi|^2+ W(\phi)\right) + \frac{1}{2} |\theta|^2\right]\,\mathrm{d}x. \end{align}\tag{17}\] We observe that the Marangoni-driven capillary term \(\mathrm{div} (\lambda(\theta) \nabla \phi \otimes \nabla \phi)\) introduces a highly nonlinear coupling between the temperature and the phase-field variable, which significantly complicates the analysis. In view of 16 , the Marangoni effect and the buoyancy force cause the loss of a standard energy dissipation structure for the system, compared to the corresponding isothermal case [21]. This poses significant difficulties in obtaining global weak solutions that are expected to be uniformly bounded in time. In addition, though physically meaningful, the nonhomogeneous Dirichlet boundary condition for the temperature introduces additional obstacles to the derivation of the necessary estimates. To achieve our goal, we adopt an implicit-explicit time-discretization scheme that maintains the uniform boundedness of the temperature along time steps (see 33 –37 ). This discretization scheme is inspired by [42] for the isothermal AGG model. By working directly with a singular potential like 14 , it also keeps the physical range of the phase-field variable, i.e., \(\phi\in [-1,1]\). The uniform \(L^{\infty}\)-bounds of the approximate solutions \((\phi,\theta)\) help us to handle some troublesome terms that are highly nonlinear. Furthermore, by exploring the coupling structure of the system, we derive a suitable discrete energy inequality that enables us to control the energy growth and obtain uniform-in-time estimates for the approximate solutions (see ?? ). Our analysis does not require any assumption of smallness on the initial temperature, thereby improving the conditions in [36], [37] for the related Navier–Stokes–Allen–Cahn system (with matched densities). In the two-dimensional setting, we obtain further regularity property of \(\theta\) using the Hölder estimate 88 (cf. [30]). Based on this improved regularity property, we can establish the uniqueness of global weak solutions by assuming that the fluid mixture has matched densities, the mobility depends only on the concentration, and the thermal diffusivity depends solely on the temperature. Whether uniqueness holds under more general assumptions remains an open question; we refer the readers to Remarks 6 and 7 for detailed discussions.

Building on the existence of global weak solutions, it is natural to study the long-time behavior. We expect that in three dimensions, results similar to those in [47] can be achieved (i.e., the longtime separation property, eventual regularity, weak-strong uniqueness, and convergence to a single equilibrium), provided that the mobility is assumed to be a positive constant. Moreover, in light of the recent contribution [49], we can apply the method therein to prove that every global weak solution converges to a single equilibrium as \(t\to \infty\) in the general case of a non-degenerate and non-constant mobility. These issues will be addressed in a forthcoming study.

Plan of the paper. In Section 2, we present the notation, mathematical tools, and main results. Section 3 is devoted to the proof of the existence of a global weak solution in both two and three dimensions. In Section 4, we prove the existence and uniqueness of global weak solutions with improved regularity for the temperature in two dimensions.

2 Main Results↩︎

2.1 Notation↩︎

We first introduce the notation and conventions used throughout this work.

Let \(B\) be a Banach space whose norm is denoted by \(\|\cdot\|_{B}\). We denote its dual by \(B^{\prime}\) and the dual product by \(\langle \cdot, \cdot \rangle_{B}\). If \(B\) is a Hilbert space, we write \((\cdot,\cdot)_B\) for the associated inner product. The bold letter \(\boldsymbol{B}\) denotes the generic space of vectors or matrices, with each component belonging to \(B\). Let \(I \subset [0, \infty]\) be an open set. We denote by \(L^p(I; B)\) the space that consists of Bochner measurable \(p\)-integrable functions (if \(p \in [1,\infty)\)) or essentially bounded functions (if \(p =\infty\)). In the case \(I=(t_1, t_2)\), we simply write \(L^p(t_1,t_2;B)\). Moreover, \(L^p_{\mathrm{uloc}}([t_0,\infty);B)\) denotes the uniformly local variant of \(L^p(t_0,\infty;B)\) consisting of all strongly measurable functions \(f:[t_0,\infty)\to B\) such that \[\|f\|_{L^p_{\mathrm{uloc}}([t_0,\infty);B)}=\sup_{t\geqslant t_0}\|f\|_{L^p(t,t+1;B)}<\infty.\] For any \(T\in (t_0,\infty)\), we simply have \(L^p_{\mathrm{uloc}}([t_0,T);B)= L^p(t_0,T;B)\). For \(p\in [1,\infty]\), \(T\in (0,\infty)\), \(f\in W^{1,p}(0,T;B)\) if and only if \(f\), \(\frac{\mathrm{d}}{\mathrm{d}t} f\in L^{p}(0,T;B)\), where \(\frac{\mathrm{d}}{\mathrm{d}t} f\) denotes the vector-valued distributional derivative of \(f\). The uniformly local space \(W^{1,p}_{\mathrm{uloc}}([t_0,\infty);B)\) is defined by replacing \(L^p(0,T; B)\) with \(L^p_{\mathrm{uloc}}([t_0,\infty);B)\). If \(p=2\), we simply set \(H^1(0,T;B) = W^{1,2}(0,T;B)\) and \(H^1_{\mathrm{uloc}}([t_0,\infty);B) = W^{1,2}_{\mathrm{uloc}}([t_0,\infty);B)\). Let \(I=[0,T]\) if \(T \in (0,\infty)\) or \(I=[0,\infty)\) if \(T=\infty\). Then \(C(I;B)\) denotes the Banach space of all continuous functions \(f:I \to B\) equipped with the supremum norm. We also denote by \(C_{w}(I;B)\) the topological vector space of all weakly continuous functions \(f:I \to B\). In addition, we denote by \(C^\infty_0(0, T; B)\) the abstract vector space of all smooth functions \(f:(0, T)\to B\) with \(\mathrm{supp} f\subset\subset (0,T)\). Given any two matrices \(X=(X_{ij})_{i,j=1}^{d}\), \(Y=(Y_{ij})_{i,j=1}^{d}\in \mathbb{R}^{d\times d}\), we have \((XY)_{ik}=\sum_{j=1}^dX_{ij}Y_{jk}\). Then we denote the Frobenius inner product by \(X:Y=\mathrm{trace}(Y^TX)=\sum_{i,k=1}^dX_{ik}Y_{ik}\) and its associated norm \(|X|=\sqrt{X:X}\).

Let \(\Omega \subset \mathbb{R}^d\), \(d\in \{2,3\}\), be a bounded domain with sufficiently smooth boundary \(\partial \Omega\). We write \(C(\overline{\Omega})\) (or \(C^{\alpha}(\overline{\Omega})\) with \(\alpha\in (0,1)\)) for the set of continuous (or \(\alpha\)-Hölder continuous) functions defined in the closure of \(\Omega\). The symbol \(C_0^{\infty}(\Omega)\) refers to the space of functions that are differentiable infinitely many times and compactly supported in \(\Omega\). For the standard Lebesgue and Sobolev spaces in \(\Omega\), we use the notation \(L^{p}(\Omega)\), \(W^{k,p}(\Omega)\) for any \(p \in [1,\infty]\) and \(k\in \mathbb{N}\), equipped with the corresponding norms \(\|\cdot\|_{L^{p}(\Omega)}\), \(\|\cdot\|_{W^{k,p}(\Omega)}\), respectively. Besides, \(W^{k,p}_0(\Omega)\) represents the closure of \(C^\infty_0(\Omega)\) in \(W^{k,p}(\Omega)\), and \(W^{-k,p^\prime}(\Omega) = (W^{k,p}_0(\Omega))^\prime\) refers to the corresponding dual space, where \(p^\prime\) is the dual exponent of \(p\). When \(p = 2\), these spaces are Hilbert spaces and we use the standard conventions \(H^{k}(\Omega) := W^{k,2}(\Omega)\), \(H^{k}_0(\Omega) := W^{k,2}_0(\Omega)\) and \(H^{-1}(\Omega)=(H^1_0(\Omega))'\). For \(s\geqslant 0\) and \(p \in [1,\infty)\), we denote by \(H^{s,p}(\Omega)\) the Bessel-potential spaces and by \(W^{s,p}(\Omega)\) the Slobodeckij spaces. We have \(H^{s,2}(\Omega) = W^{s,2}(\Omega)\) for all \(s\), but for \(p\neq 2\) the identity \(H^{s,p}(\Omega) = W^{s,p}(\Omega)\) is only true if \(s\in \mathbb{N}\). In addition, for \(s\in \mathbb{N}\), \(H^{s,p}(\Omega)\) and \(W^{s,p}(\Omega)\) coincide with the usual Sobolev spaces. The corresponding function spaces on the boundary \(\partial\Omega\) are defined using local charts. For convenience, we define \[H_{\mathbf{n}}^2(\Omega) = \{ f \in H^2(\Omega) : \partial_{ \mathbf{n}} f = 0 \text{ on } \partial \Omega \},\] and denote the inner product and norm of \(L^2(\Omega)\) by \((\cdot, \cdot)\) and \(\|\cdot\|\) without ambiguity.

For every \(f\in (H^1(\Omega))^{\prime}\), its generalized mean value over \(\Omega\) is defined as \(\overline{f}=|\Omega|^{-1}\langle f,1\rangle_{H^1(\Omega)}\). If \(f\in L^1(\Omega)\), then the mean value is given by \(\overline{f}=|\Omega|^{-1}\int_\Omega f \,\mathrm{d}x\). For any \(M\in \mathbb{R}\), we set \[L^2_{(M)}(\Omega):= \left\{f\in L^2(\Omega):\overline{f} =M\right\}.\] The orthogonal projection \(P_0: L^2(\Omega)\to L^2_{(0)}(\Omega)\) is defined as \(P_0 f= f- \overline{f}\). Besides, we introduce the linear spaces \[\begin{align} & V_{(0)} := \left\{f \in H^1(\Omega) : \overline{f} = 0 \right\}, \quad V_{(0)}^{\prime} := \{f \in (H^1(\Omega))^{\prime} : \overline{f} =0\}, \end{align}\] and recall the well-known Poincaré–Wirtinger inequality: \[\left\|f-\overline{f}\right\|\leqslant C_{PW} \|\nabla f\|,\quad \forall\, f\in H^1(\Omega),\] where the positive constant \(C_{PW}\) depends only on \(\Omega\). Consider the Neumann problem \[\begin{align} \begin{cases} -\Delta u = f, \quad &\text{in}~\Omega, \\ \partial_\mathbf{n} u = 0, \quad &\text{on}~\partial \Omega. \end{cases} \end{align}\] We define the solution operator \(\mathcal{G}: V_{(0)}^{\prime} \to V_{(0)}\) in the following sense: for every \(f \in V_{(0)}^{\prime}\), \(\mathcal{G}f \in V_{(0)}\) is the unique function satisfying \[\begin{align} \label{Pre-G} \langle f , v \rangle_{H^1(\Omega)} = ( \nabla \mathcal{G}f, \nabla v), \quad \forall\, v \in V_{(0)}. \end{align}\tag{18}\] For \(f \in V_{(0)}^{\prime}\), it is easy to check that \(\| \nabla \mathcal{G}f \|\) is a norm on \(V_{(0)}^{\prime}\) equivalent to the standard one. Thus, we also write it as \(\|\cdot\|_{V_{(0)}^{\prime}}\). In order to handle the nonconstant mobility, we also consider the following Neumann problem \[\begin{align} \label{Gq} \begin{cases} -\mathrm{div} ( m(q) \nabla u) = f, \quad &\text{in}~\Omega, \\ m(q) \partial_\mathbf{n} u = 0, \quad &\text{on}~\partial \Omega. \end{cases} \end{align}\tag{19}\] Here, \(q: \Omega \to [-1,1]\) is a measurable function. Define the solution operator \(\mathcal{G}_q\) as \[\begin{align} \label{Pre-Gq} \langle f , v \rangle_{H^1(\Omega)} = ( m(q) \nabla \mathcal{G}_q f, \nabla v), \quad \forall\, v \in V_{(0)}. \end{align}\tag{20}\] As it has been shown in [50], for \(f \in V_{(0)}^{\prime}\), \(\| \nabla \mathcal{G}_q f \|\) and \(\| \nabla \mathcal{G}f \|\) are equivalent norms to the \((H^1(\Omega))^{\prime}\)-norm. Further properties of the operator \(\mathcal{G}_q\) can be found in, e.g., [50]. Suppose that \(m \in C^1([-1,1])\) and \(q \in H^2(\Omega)\). Then for any \(f \in L^2_{(0)}(\Omega)\), we find \(\mathcal{G}_q f\in H_{\mathbf{n}}^2(\Omega)\) and the following estimate when the spatial dimension is two (see [50]): \[\begin{align} \label{GqH2} \| \mathcal{G}_q f \|_{H^2(\Omega)} \leqslant C \big( \| \nabla q \| \| q \|_{H^2(\Omega)} \| \nabla \mathcal{G}_q f \| + \| f \| \big), \quad d=2. \end{align}\tag{21}\] Here, the positive constant \(C\) depends only on \(\Omega\) and the function \(m\).

Next, we recall some function spaces for the Navier–Stokes equations. Let \(\boldsymbol{C}^\infty_{0,\sigma}(\Omega)\) be the space of divergence-free vector fields in \((C^\infty_0(\Omega))^d\). We define \(\boldsymbol{H}_{\sigma}\) and \(\boldsymbol{V}_{\sigma}\) as the closure of \(\boldsymbol{C}^\infty_{0,\sigma}(\Omega)\) with respect to the \(\boldsymbol{L}^2\) and \(\boldsymbol{H}^1\) norms, respectively. The space \(\boldsymbol{V}_{\sigma}\) is equipped with the scalar product \((\boldsymbol{u},\boldsymbol{v})_{\boldsymbol{V}_{\sigma}}:=(\nabla \boldsymbol{u},\nabla \boldsymbol{v})\) for all \(\boldsymbol{u},\, \boldsymbol{v} \in \boldsymbol{V}_{\sigma}\) and the norm \(\|\boldsymbol{u}\|_{\boldsymbol{V}_{\sigma}}=\|\nabla \boldsymbol{u}\|\). We recall the classical Korn’s inequality \[\begin{align} \|\nabla \boldsymbol{u}\| \leqslant \sqrt{2} \|D\boldsymbol{u}\|\leqslant \sqrt{2}\|\nabla \boldsymbol{u}\|,\quad \forall\, \boldsymbol{u}\in \boldsymbol{V}_{\sigma}. \label{Korn} \end{align}\tag{22}\] Besides, Poincaré’s inequality entails that \[\begin{align} \|\boldsymbol{u}\| \leqslant C_P \|\nabla \boldsymbol{u}\|,\quad \forall\, \boldsymbol{u}\in \boldsymbol{V}_{\sigma}, \label{Poin} \end{align}\tag{23}\] where the positive constant \(C_P\) depends only on \(\Omega\). It is well known that \(\boldsymbol{L}^2(\Omega)\) can be decomposed into \(\boldsymbol{H}_{\sigma}\oplus\boldsymbol{G}(\Omega)\), where \(\boldsymbol{G}(\Omega):=\{\boldsymbol{f}\in\boldsymbol{L}^2(\Omega): \exists\, g\in H^1(\Omega),\;\boldsymbol{f}=\nabla g\}\). Then we introduce the Helmholtz–Leray projection in the space of divergence-free functions \(\boldsymbol{P}:\boldsymbol{L}^2(\Omega)\to \boldsymbol{H}_{\sigma}\). Recall the Stokes operator \(\boldsymbol{S}: \boldsymbol{V}_{\sigma}\cap\boldsymbol{H}^2(\Omega)\to\boldsymbol{H}_{\sigma}\) such that \[(\boldsymbol{S}\boldsymbol{u},\boldsymbol{\zeta})=(\nabla \boldsymbol{u},\nabla\boldsymbol{\zeta}),\quad \forall\, \boldsymbol{\zeta} \in \boldsymbol{V}_{\sigma}, \nonumber\] with the domain \(\mathcal{D}(\boldsymbol{S})= \boldsymbol{V}_{\sigma}\cap\boldsymbol{H}^2(\Omega)\) (see, e.g., [51]). The operator \(\boldsymbol{S}\) is a canonical isomorphism from \(\boldsymbol{V}_{\sigma}\) to \(\boldsymbol{V}_{\sigma}^{\prime}\) with the inverse \(\boldsymbol{S}^{-1}:\boldsymbol{V}_{\sigma}^{\prime}\to\boldsymbol{V}_{\sigma}\). For any \(\boldsymbol{f}\in \boldsymbol{V}_{\sigma}^{\prime}\), there is a unique \(\boldsymbol{u}=\boldsymbol{S}^{-1}\boldsymbol{f}\in\boldsymbol{V}_{\sigma}\) such that \[(\nabla(\boldsymbol{S}^{-1}\boldsymbol{f}),\nabla \boldsymbol{\zeta})=\langle\boldsymbol{f},\boldsymbol{\zeta}\rangle_{\boldsymbol{V}_{\sigma}}, \quad \forall\, \boldsymbol{\zeta} \in \boldsymbol{V}_{\sigma}, \nonumber\] and \(\|\nabla(\boldsymbol{S}^{-1}\boldsymbol{f})\|=\langle\boldsymbol{f},\boldsymbol{S}^{-1}\boldsymbol{f} \rangle_{\boldsymbol{V}_{\sigma}}^{\frac{1}{2}}\) is an equivalent norm on \(\boldsymbol{V}_{\sigma}^{\prime}\). In addition, there exists a positive constant \(C\) such that \(\|\boldsymbol{u}\|_{\boldsymbol{H}^2(\Omega)}\leqslant C\|\boldsymbol{S} \boldsymbol{u}\|\) for any \(\boldsymbol{u}\in \mathcal{D}(\boldsymbol{S})\).

Throughout this paper, the symbols \(C\), \(C_i\), \(i\in \mathbb{N}\), denote generic positive constants that may depend on coefficients of the system, norms of the initial data, \(\Omega\) and time. Their values may change from line to line, and specific dependencies will be pointed out when necessary.

2.2 Statement of results↩︎

First, we introduce the assumptions used in the subsequent analysis.

• \(\Omega \subset \mathbb{R}^d\), \(d\in\{2,3\}\), is a bounded domain with \(C^3\)-boundary \(\partial \Omega\).

• The viscosity \(\nu: \mathbb{R}^2 \to \mathbb{R}^+\) belongs to \(C^1(\mathbb{R}^2)\) and \[\begin{align} \text{there exists}\;\underline{\nu} >0\;\text{such that}\;\; \nu(s_1,s_2)\geqslant \underline{\nu}, \quad \forall\, (s_1,s_2) \in \mathbb{R}^2. \end{align}\]

• The mobility \(m:\mathbb{R}^2\to\mathbb{R}^+\) belongs to \(C^1(\mathbb{R}^2)\) and \[\begin{align} \text{there exists}\;\underline{m}>0\;\text{such that}\;\;m(s_1,s_2)\geqslant \underline{m}, \quad \forall\, (s_1,s_2) \in \mathbb{R}^2. \end{align}\]

• The thermal diffusivity \(\kappa:\mathbb{R}^2\to\mathbb{R}^+\) belongs to \(C^2(\mathbb{R}^2)\) and \[\begin{align} \text{there exists}\;\underline{\kappa}>0\;\text{such that}\;\;\kappa(s_1,s_2)\geqslant \underline{\kappa}, \quad \forall\, (s_1,s_2) \in \mathbb{R}^2. \end{align}\]

• The singular potential \(W\) belongs to the class of functions \(C([-1,1])\cap C^{2}(-1,1)\) and it can be decomposed into the following form \[W(s)=W_1(s)+W_2(s).\nonumber\] The singular convex part \(W_1\) fulfills \[\lim_{s\to \pm 1} W_1^{\prime}(s)=\pm \infty ,\quad \text{and}\;\; W_1^{\prime \prime}(s)\geqslant c_0,\quad \forall\, s\in (-1,1),\nonumber\] for some strictly positive constant \(c_0\). Besides, we make the extension \(W_1(s)=+\infty\) for any \(s\notin[-1,1]\). Concerning the regular (possibly concave) part \(W_2\), we assume that \[\text{there exists}\;c_W >0 \;\;\text{such that}\;\;|W^{\prime \prime}_2(s)|\leqslant c_W,\quad \forall\, s\in\mathbb{R}.\]

• The averaged density \(\rho\) is given by 8 , the relative flux \(\mathbf{J}\) satisfies 10 and the buoyancy force \(\boldsymbol{f}_{\mathrm{b}}\) takes the form of 9 . The Cauchy stress tensor \(\boldsymbol{\sigma}\) is defined as in 11 with the surface tension coefficient \(\lambda\) determined by 12 . \(\rho_1\), \(\rho_2\), \(\lambda_0\), \(a\), \(b\), \(\alpha\), \(g\) are given positive constants.

We are now in a position to state the main results.

Next, in the two-dimensional case, we establish the uniqueness of global weak solutions under the assumptions that the fluid mixture has matched densities, the thermal diffusivity depends only on the temperature, the mobility depends only on the concentration, and the initial temperature is more regular.

3 Existence of Global Weak Solutions↩︎

In this section, we prove Theorem 1 on the existence of global weak solutions to problem 1 –7 . The proof is based on an implicit-explicit time-discretization scheme in the spirit of [42], [54], with suitable modifications according to new structures of the coupled system. Without loss of generality, we present the argument for the three-dimensional case while pointing out necessary changes in two dimensions.

3.1 sec:Preliminaries↩︎

To begin with, let us recall some properties of subgradients related to the Ginzburg–Landau free energy (see [55], [56]). Concerning the convex part of the free energy (cf. Remark 1) \[\begin{align} E(\phi) = \int_{\Omega} \left( \frac{1}{2} |\nabla \phi|^{2} + F(\phi) \right) \mathrm{d}x, \label{G-L} \end{align}\tag{24}\] we define \[\begin{align} \label{convex-part-operator} \left.\widetilde{E}(\phi)= \left\{\begin{array}{ccc} E(\phi), &\text{for } \phi \in \mathcal{D}(\widetilde{E}),\\ +\infty, &\text{else}, \end{array}\right. \right. \end{align}\tag{25}\] with the effective domain \[\mathcal{D}(\widetilde{E}) := \left\{ \phi \in H^1(\Omega)\;:\;-1 \leqslant \phi \leqslant 1 \text{ a.e. in } \Omega \right\}.\] Then we define the subgradient \[\partial \widetilde{E}(\phi) = - \Delta \phi + F^{\prime}(\phi),\] with the domain \[\mathcal{D}(\partial\widetilde{E}) =\left\{\phi \in H^{2}_\mathbf{n}(\Omega)\;: \; F^{\prime}(\phi)\in L^{2}(\Omega),\quad F^{\prime\prime}(\phi)|\nabla \phi|^{2}\in L^{1}(\Omega)\right\}.\] Similarly, for any given \(M\in (-1,1)\), we define \[\begin{align} \label{convex-part-operator-M} \left.\widetilde{E}_M(\phi)= \left\{\begin{array}{ccc} E(\phi), &\text{for } \phi \in \mathcal{D}(\widetilde{E}_M),\\ +\infty, &\text{else}, \end{array}\right. \right. \end{align}\tag{26}\] with the effective domain \[\mathcal{D}(\widetilde{E}_M) := \left\{ \phi \in H^1(\Omega)\cap L^2_{(M)}(\Omega)\;:\;-1 \leqslant \phi \leqslant 1 \text{ a.e. in } \Omega \right\},\] as well as the subgradient \[\partial \widetilde{E}_M(\phi) = - \Delta \phi + P_0(F^{\prime}(\phi)),\] with the domain \[\mathcal{D}(\partial\widetilde{E}_M) =\left\{\phi \in H^{2}_\mathbf{n}(\Omega)\cap L^2_{(M)}(\Omega)\;: \; F^{\prime}(\phi)\in L^{2}(\Omega),\quad F^{\prime\prime}(\phi)|\nabla \phi|^{2}\in L^{1}(\Omega)\right\}.\] According to [55], the following estimates hold \[\begin{align} & \| \phi \|_{H^2(\Omega)}^2 + \| F^{\prime}(\phi) \|^2 + \int_{\Omega} F^{\prime \prime}(\phi) | \nabla \phi |^2 \: \mathrm{d}x\leqslant C \left( \| \partial \widetilde{E}(\phi) \|^2 + \| \phi \|^2 + 1 \right), \label{es-pE-1} \end{align}\tag{27}\] for all \(\phi\in \mathcal{D}(\partial \widetilde{E})\), and \[\begin{align} & \| \phi \|_{H^2(\Omega)}^2 + \| F^{\prime}(\phi) \|^2 + \int_{\Omega} F^{\prime \prime}(\phi) | \nabla \phi |^2 \: \mathrm{d}x\leqslant C \left( \| \partial \widetilde{E}_M(\phi) \|^2 + \| \phi \|^2 + 1 \right), \label{es-pE-2} \end{align}\tag{28}\] for all \(\phi\in \mathcal{D}(\partial \widetilde{E}_M)\). The positive constant \(C\) in 27 (resp. 28 ) is independent of \(\phi\in \mathcal{D}(\partial\widetilde{E})\) (resp. \(\phi\in \mathcal{D}(\partial\widetilde{E}_M)\)).

Next, we recall some results for the following Neumann problem with singular nonlinearity, which will be useful in the subsequent analysis (see [55], [57]–[59]): \[\begin{align} \begin{cases} - \Delta \phi + F^{\prime}(\phi) = f,\quad \text{in}\;\Omega,\\ \partial_{\mathbf{n}} \phi=0,\qquad \qquad\quad \;\text{on}\;\partial\Omega. \end{cases} \label{sing-phi} \end{align}\tag{29}\]

3.2 The implicit-explicit time-discretization scheme↩︎

Let us now present the implicit-explicit time-discretization scheme for problem 1 –7 . In order to handle the nonhomogeneous Dirichlet boundary condition for \(\theta\), we choose \(\Theta_{\mathrm{b}}\) as the harmonic extension of the boundary datum \(\theta_{\mathrm{b}}\), which is inspired by e.g., [32] and [52]: \[\begin{cases} -\Delta \Theta_{\mathrm{b}}=0,\quad \text{in}\;\Omega,\\ \Theta_{\mathrm{b}}=\theta_{\mathrm{b}},\quad\quad \, \text{on}\;\partial\Omega. \end{cases} \label{Theta95b}\tag{30}\] Since \(\theta_{\mathrm{b}}\in H^\frac{5}{2}(\partial\Omega)\), we find \[\Theta_{\mathrm{b}}\in H^3(\Omega)\quad \text{and}\quad \|\Theta_{\mathrm{b}}\|_{H^3(\Omega)}\leqslant C\|\theta_{\mathrm{b}}\|_{H^\frac{5}{2}(\partial\Omega)},\] for some constant \(C>0\) depending only on \(\Omega\). Introducing the new variable \[\vartheta=\theta - \Theta_{\mathrm{b}},\] we can reformulate the heat equation 5 as follows (cf. [52]): \[\begin{align} & \partial_t \vartheta + \boldsymbol{u}\cdot \nabla \vartheta - \mathrm{div}\,(\kappa(\phi,\vartheta+ \Theta_{\mathrm{b}}) \nabla \vartheta) = -\boldsymbol{u}\cdot \nabla \Theta_{\mathrm{b}} + \mathrm{div}\,(\kappa(\phi,\vartheta+ \Theta_{\mathrm{b}}) \nabla \Theta_{\mathrm{b}}). \label{Eq:vartheta} \end{align}\tag{31}\] Moreover, \(\vartheta\) satisfies the following homogeneous Dirichlet boundary condition and the initial condition \[\begin{align} \vartheta=0\quad \text{on}\;\partial\Omega\times (0,T),\qquad \vartheta|_{t=0}=\vartheta_0:=\theta_0- \Theta_{\mathrm{b}}\quad \text{in}\;\Omega. \label{vartheta:bdini} \end{align}\tag{32}\] Due to the Sobolev embedding theorem \(H^2(\Omega)\hookrightarrow L^\infty(\Omega)\) for \(d=2,3\), we have \(\vartheta_0\in L^\infty(\Omega)\) since \(\theta_0\in L^\infty(\Omega)\) and \(\Theta_{\mathrm{b}}\in H^2(\Omega)\).

For any given positive integer \(N\), we take the time step as \(h=\frac{1}{N}\). Next, for every \(k \in \mathbb{N}\), let \(\boldsymbol{u}^k\in \boldsymbol{H}_\sigma\), \(\phi^k\in H^2_\mathbf{n}(\Omega)\) with \(\phi^k\in [-1,1]\) almost everywhere in \(\Omega\), \(\overline{\phi^k}\in (-1,1)\), and \(\vartheta^k\in H^2(\Omega) \cap H^1_0(\Omega)\). Define \(\rho^k=\rho(\phi^k)\) by 8 . Then, we look for a quadruple \[(\boldsymbol{u}, \phi, \mu, \vartheta) =\big(\boldsymbol{u}^{k+1}, \phi^{k+1}, \mu^{k+1}, \vartheta^{k+1}\big) \in \boldsymbol{V}_\sigma \times \mathcal{D}(\partial\widetilde{E}) \times H^2_{\mathbf{n}}(\Omega) \times ( H^2(\Omega)\cap H^1_0(\Omega))\] as a solution of the following time-discrete system at the time step \(k+1\): \[\begin{align} &\left( \frac{\rho\boldsymbol{u}- \rho^k\boldsymbol{u}^k}{h} , \boldsymbol{v}\right) + \big(\mathrm{div}(\rho^k \boldsymbol{u}\otimes \boldsymbol{u}), \boldsymbol{v}\big) + (\mathrm{div}(\boldsymbol{u}\otimes \mathbf{J}), \boldsymbol{v}) \notag \\ &\qquad + \big(2 \nu(\phi^k, \vartheta^k+ \Theta_{\mathrm{b}}) D \boldsymbol{u}, D \boldsymbol{v}\big) - \lambda_0 a (\mu \nabla \phi^k, \boldsymbol{v}) \notag \\ &\quad = - \lambda_0 b \big((\vartheta^k+ \Theta_{\mathrm{b}}) (\nabla \phi \otimes \nabla \phi), \nabla \boldsymbol{v}\big) + \big(\boldsymbol{f}_{\mathrm{b}}(\phi^k, \vartheta^k+ \Theta_{\mathrm{b}}), \boldsymbol{v}\big), \quad &&\forall\, \boldsymbol{v}\in \boldsymbol{V}_{\sigma}, \tag{33}\\ & \mathbf{J}= -\frac{\rho_2-\rho_1}{2}m(\phi^k,\vartheta^k+ \Theta_{\mathrm{b}})\nabla \mu, \quad &&\text{a.e. in } \Omega, \tag{34}\\ & \frac{\phi - \phi^k}{h} + \boldsymbol{u}\cdot \nabla \phi^k = \mathrm{div}(m(\phi^k,\vartheta^k+ \Theta_{\mathrm{b}}) \nabla \mu), \quad &&\text{a.e. in } \Omega, \tag{35} \\ & \mu + c_W (\phi+\phi^k) = - \Delta \phi + F^{\prime}(\phi), \quad &&\text{a.e. in } \Omega, \tag{36} \\ & \frac{\vartheta - \vartheta^k}{h} + \boldsymbol{u}\cdot \nabla \vartheta - \mathrm{div} ( \kappa(\phi^k,\vartheta^k+ \Theta_{\mathrm{b}}) \nabla \vartheta) \notag \\ &\quad = - \boldsymbol{u}\cdot \nabla \Theta_{\mathrm{b}} + \mathrm{div} ( \kappa(\phi^k,\vartheta^k+ \Theta_{\mathrm{b}}) \nabla \Theta_{\mathrm{b}}), \quad &&\text{a.e. in } \Omega. \tag{37} \end{align}\]

Applying the Stampacchia truncation method for equation 37 , we obtain the following result:

Next, we show that the discrete solution \(\phi\) satisfies the property of mass conservation.

The following lemma gives a refinement of Lemma 1 under the specific choice \(f=\mu+c_W(\phi+\phi^k)\) as in the discrete approximation 36 .

Let \(N\) be an arbitrarily given positive integer. We establish the existence of a solution to the discrete problem 33 –37 .

3.3 Existence of weak solutions to the continuous problem↩︎

The proof of Theorem 1 consists of several steps.

Step 1. Construction of the initial data for induction. Let \(N\) be a given positive integer and \(h=\frac{1}{N}\). We first construct the initial data \((\boldsymbol{u}^0,\phi^0,\vartheta^0)\) for the induction associated with the discrete problem 33 –37 .

Given \(\phi_0 \in H^{1}(\Omega)\) satisfying \(\|\phi_0\|_{L^{\infty}(\Omega)} \leqslant 1\), \(|\overline{\phi_0}|<1\), as in [42], we take \[\phi_0^N=\varphi\left(\frac{1}{N}\right),\] where \(\varphi\) is the (unique) solution of the following linear parabolic problem \[\begin{cases} \partial_t \varphi = \Delta \varphi, \qquad \;\,\text{in}\;\Omega\times (0,\infty),\\ \partial_\mathbf{n}\varphi=0,\qquad\;\;\;\; \text{on}\;\partial\Omega\times (0,\infty),\\ \varphi|_{t=0}=\phi_0,\qquad \;\text{in}\;\Omega. \end{cases} \notag\] By the classical theory of linear parabolic equations, we easily find \(\varphi\in C([0,\infty);H^1(\Omega))\cap C((0,\infty);H^2_\mathbf{n}(\Omega))\) such that \(\|\varphi(t)\|_{L^\infty(\Omega)}\leqslant 1\), \(\overline{\varphi(t)}=\overline{\phi_0}\in (-1,1)\) for all \(t\geqslant 0\). As a consequence, we have \[\phi_0^N\to \phi_0\quad \text{in}\;H^1(\Omega)\;\;\text{as}\;N\to \infty,\] and due to the convexity of \(F\), \[F(\phi_0^N)\to F(\phi_0) \quad \text{in}\;\;L^1(\Omega).\]

Given \(\theta_0\in L^\infty(\Omega)\), we choose \[\theta_0^N=\varpi^N, \qquad \vartheta_0^N=\theta_0^N- \Theta_{\mathrm{b}},\] where \(\varpi^N\) is the (unique) solution to the following linear elliptic problem \[\begin{cases} -\dfrac{1}{N} \Delta \varpi^N +\varpi^N = \theta_0,\quad \text{in}\;\Omega,\\ \varpi^N=\theta_\mathrm{b},\qquad \qquad \qquad \;\;\text{on}\;\partial \Omega, \end{cases} \notag\] and the function \(\Theta_\mathrm{b}\) is given by 30 . Applying the classical theory for linear elliptic equations, we have \(\varpi^N-\Theta_\mathrm{b}\in H^2(\Omega)\cap H^1_0(\Omega)\). In addition, an argument similar to that for Lemma 2 yields \[\min\Big\{\operatorname*{ess\,inf}_{\Omega}\theta_0,\,\min_{\partial\Omega}\theta_{\mathrm{b}}\Big\} \leqslant \varpi^N(x) \leqslant \max\Big\{\operatorname*{ess\,sup}_{\Omega}\theta_0,\,\max_{\partial\Omega}\theta_{\mathrm{b}}\Big\},\quad \text{a.e. in}\;\overline{\Omega},\] which is consistent with ?? . Testing the equation for \(\varpi^N\) by \(\varpi^N-\Theta_\mathrm{b}\), using the Cauchy–Schwarz inequality and Poincaré’s inequality, we obtain \[\begin{align} \frac{1}{N}\|\nabla (\varpi^N-\Theta_\mathrm{b})\|^2+ \| \varpi^N-\Theta_\mathrm{b}\|^2 \leqslant \|\theta_0-\Theta_\mathrm{b}\|\|\varpi^N-\Theta_\mathrm{b}\| \leqslant C\|\theta_0-\Theta_\mathrm{b}\|\|\nabla (\varpi^N-\Theta_\mathrm{b})\|, \end{align}\] which implies \[\| \varpi^N-\Theta_\mathrm{b}\|\leqslant \|\theta_0-\Theta_\mathrm{b}\|,\qquad \frac{1}{N}\|\nabla (\varpi^N-\Theta_\mathrm{b})\| \leqslant C\|\theta_0-\Theta_\mathrm{b}\|.\] Next, testing the equation for \(\varpi^N\) by \((-\Delta_D)^{-1}(\varpi^N-\theta_0)\in H^2(\Omega)\cap H^1_0(\Omega)\) (here \((-\Delta_D)^{-1}\) denotes the inverse Laplacian subject to the homogeneous Dirichlet boundary condition), we infer from the fact \(\Delta \Theta_{\mathrm{b}} = 0\) in \(\Omega\), Cauchy–Schwarz inequality and the above estimates that \[\begin{align} \|\varpi^N-\theta_0\|_{H^{-1}(\Omega)}^2 & = - \frac{1}{N} (\varpi^N-\Theta_\mathrm{b},\varpi^N-\theta_0) \\ &\leqslant \frac{1}{N} \|\varpi^N-\Theta_\mathrm{b}\| (\|\varpi^N-\Theta_\mathrm{b}\|+ \|\Theta_\mathrm{b}-\theta_0\|)\\ &\leqslant \frac{2}{N}\|\theta_0-\Theta_\mathrm{b}\|^2. \end{align}\] As a consequence, it holds \[\varpi^N\to \theta_0\quad \text{in}\;H^{-1}(\Omega)\;\;\text{as}\;\;N\to \infty.\] By the uniform \(L^2\)-bound for \(\varpi^N\) and the uniqueness of the limit, we also have \[\varpi^N\rightharpoonup \theta_0\quad \text{in}\;L^2(\Omega)\;\;\text{as}\;\;N\to \infty.\] Hence, we observe that \[\begin{align} \|\varpi^N-\theta_0\|^2 &= (\varpi^N-\Theta_\mathrm{b}, \varpi^N-\theta_0)- (\theta_0-\Theta_\mathrm{b}, \varpi^N-\theta_0) \\ & = \frac{1}{N}(\varpi^N-\Theta_\mathrm{b}, \Delta (\varpi^N-\Theta_\mathrm{b}))- (\theta_0-\Theta_\mathrm{b}, \varpi^N-\theta_0) \\ & = \underbrace{-\frac{1}{N}\|\nabla (\varpi^N-\Theta_\mathrm{b})\|^2}_{\leqslant 0} - (\theta_0-\Theta_\mathrm{b}, \varpi^N-\theta_0). \end{align}\] This together with the weak convergence of \(\varpi^N\) in \(L^2(\Omega)\) yields \[\begin{align} 0\leqslant \limsup_{N\to \infty} \|\varpi^N-\theta_0\|^2 \leqslant \limsup_{N\to \infty} \big[- (\theta_0-\Theta_\mathrm{b}, \varpi^N-\theta_0)\big] =0, \notag \end{align}\] that is, the strong convergence of \(\varpi^N\) in \(L^2(\Omega)\). Hence, we have \[\theta^N_0\to \theta_0 \quad \text{and}\quad \vartheta^N_0\to \theta_0-\Theta_\mathrm{b}\quad \text{in}\;L^2(\Omega)\;\;\text{as}\;\;N\to \infty.\] Therefore, for the discrete problem 33 –37 , we set the initial data \[(\boldsymbol{u}^0,\phi^0,\vartheta^0)=(\boldsymbol{u}_0, \phi^N_0,\vartheta^N_0).\]

Step 2. Construction of approximate solutions and uniform estimates. Like in [42], we define the piecewise constant interpolant \(f^{N}(t)\) on \([-h,\infty )\) through \(f^{N}(t)=f^{k}\) for \(t\in [(k-1)h,kh)\), where \(k\in \mathbb{N}\) and \(f\in \{\boldsymbol{u},\phi,\mu,\vartheta\}\). We also define \(\rho^{N}:=\rho( \phi^{N})\). Due to the implicit formulation for \(\mu\), we simply set \(\mu^N(t)=0\) for \(t\in [-h,0)\). Moreover, we introduce the time shifts \(f_{h}:=f(t-h)\) and the notation for time differences as well as difference quotients: \[\begin{align} \left( \Delta _{h}^{+}f\right) (t):=f(t+h)-f(t),\qquad & \partial_{t,h}^{+}f(t):=\frac{1}{h}\left( \Delta _{h}^{+}f\right)(t), \\ \left( \Delta _{h}^{-}f\right) (t):=f(t)-f(t-h),\qquad & \partial_{t,h}^{-}f(t):=\frac{1}{h}\left( \Delta _{h}^{-}f\right)(t). \end{align}\]

From 33 –37 , we can derive the corresponding time-continuous equations. For an arbitrary vector \(\boldsymbol{v}\in C_{0}^{\infty }\big( 0,\infty;\boldsymbol{C}_{0,\sigma}^{\infty}(\Omega)\big)\), we choose \(\widetilde{\boldsymbol{v}}:=\int_{kh}^{(k+1)h}\boldsymbol{v}\,\mathrm{d}t\) as a test function in the weak formulation 33 and sum over \(k\in \mathbb{Z}^+\) to get \[\begin{align} &\int_0^\infty\!\!\int_\Omega \left[- (\rho^N\boldsymbol{u}^N) \cdot \partial_{t,h}^{+} \boldsymbol{v} - (\rho_h^N \boldsymbol{u}^N \otimes \boldsymbol{u}^N):\nabla \boldsymbol{v} - (\boldsymbol{u}^N\otimes \mathbf{J}^N): \nabla \boldsymbol{v}\right]\,\mathrm{d}x\mathrm{d}t \notag \\ &\qquad + \int_0^\infty\!\!\int_\Omega 2\nu(\phi_h^N, \vartheta_h^N+ \Theta_{\mathrm{b}}) D \boldsymbol{u}^N : D \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t - \int_0^\infty\!\!\int_\Omega \lambda_0 a \mu^N \nabla \phi_h^N\cdot \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t \notag \\ &\quad = - \int_0^\infty\!\!\int_\Omega \lambda_0 b (\vartheta_h^N+ \Theta_{\mathrm{b}}) (\nabla \phi^N \otimes \nabla \phi^N): \nabla \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t \notag\\ &\qquad + \int_0^\infty\!\!\int_\Omega \boldsymbol{f}_{\mathrm{b}}(\phi_h^N, \vartheta_h^N+ \Theta_{\mathrm{b}})\cdot \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t, \label{time-conti-u} \end{align}\tag{66}\] where \[\begin{align} & \mathbf{J}^N= -\frac{\rho_2-\rho_1}{2}m(\phi_h^N, \vartheta_h^N+ \Theta_{\mathrm{b}})\nabla \mu^N,\quad &&\text{a.e. in } \Omega\times (0,\infty), \tag{67}\\ & \mu^N + c_W (\phi^N+\phi_h^N) = - \Delta \phi^N + F^{\prime}(\phi^N), \quad &&\text{a.e. in } \Omega\times (0,\infty). \tag{68} \end{align}\] Concerning the first term of 66 , we have also used the following integration by parts with respect to time \[\begin{align} \int_0^\infty\!\!\int_\Omega \partial_{t,h}^{-} (\rho^N\boldsymbol{u}^N )\cdot \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t= -\int_0^\infty\!\!\int_\Omega (\rho^N\boldsymbol{u}^N) \cdot \partial_{t,h}^{+} \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t. \label{rhou-intgra} \end{align}\tag{69}\] Analogously, we have \[\begin{align} & \int_0^\infty\!\!\int_\Omega \left(\partial_{t,h}^{-} \phi^N \zeta - (\phi_h^N \boldsymbol{u}^N) \cdot \nabla \zeta \right)\,\mathrm{d}x\mathrm{d}t + \int_0^\infty\!\!\int_\Omega m(\phi_h^N,\vartheta_h^N + \Theta_{\mathrm{b}}) \nabla \mu^N \cdot \nabla \zeta\,\mathrm{d}x\mathrm{d}t = 0, \tag{70} \\ & \int_0^\infty\!\!\int_\Omega \left(\partial_{t,h}^{-} \vartheta^N \xi - (\vartheta^N \boldsymbol{u}^N) \cdot \nabla \xi \right)\,\mathrm{d}x\mathrm{d}t + \int_0^\infty\!\!\int_\Omega \kappa(\phi_h^N,\vartheta_h^N + \Theta_{\mathrm{b}}) \nabla \vartheta^N \cdot \nabla \xi\,\mathrm{d}x\mathrm{d}t \notag \\ &\quad = \int_0^\infty\!\!\int_\Omega \left[ (\Theta_{\mathrm{b}} \boldsymbol{u}^N) \cdot \nabla \xi - \kappa(\phi_h^N,\vartheta_h^N+ \Theta_{\mathrm{b}}) \nabla \Theta_{\mathrm{b}} \cdot \nabla \xi \right]\,\mathrm{d}x\mathrm{d}t, \tag{71} \end{align}\] for all \(\zeta\in C_0((0,\infty); C^\infty(\overline{\Omega}))\), \(\xi \in C_0((0,\infty);C_0^\infty(\overline{\Omega}))\).

Define \[\widehat{E}_{\mathrm{tot}}(\boldsymbol{u}^{k},\phi^{k},\vartheta^{k})= (1+C_1 h)^k E_{\mathrm{tot}}(\boldsymbol{u}^{k},\phi^{k},\vartheta^{k}),\quad k\in \mathbb{N}.\] From the discrete energy inequality ?? , we obtain \[\begin{align} &\frac{\widehat{E}_{\mathrm{tot}}(\boldsymbol{u}^{k+1},\phi^{k+1},\vartheta^{k+1})- \widehat{E}_{\mathrm{tot}}(\boldsymbol{u}^{k},\phi^{k},\vartheta^{k})}{h} \notag \\ &\quad = (1+C_1 h)^k\left( \frac{E_{\mathrm{tot}}(\boldsymbol{u}^{k+1},\phi^{k+1},\vartheta^{k+1})- E_{\mathrm{tot}}(\boldsymbol{u}^{k},\phi^{k},\vartheta^{k})}{h} + C_1E_{\mathrm{tot}}(\boldsymbol{u}^{k+1},\phi^{k+1},\vartheta^{k+1})\right) \notag\\ &\quad \leqslant C_2(1+C_1 h)^k. \label{discre-ineq-1} \end{align}\tag{72}\] Summation over \(k\) yields the following result \[\begin{align} \widehat{E}_{\mathrm{tot}}(\boldsymbol{u}^{k},\phi^{k},\vartheta^{k}) &\leqslant \widehat{E}_{\mathrm{tot}}(\boldsymbol{u}^{0},\phi^{0},\vartheta^{0}) + C_2h \sum_{j=0}^{k-1} (1+C_1 h)^j,\quad \forall\, k\in \mathbb{Z}^+, \notag \end{align}\] which implies \[\begin{align} E_{\mathrm{tot}}(\boldsymbol{u}^{k},\phi^{k},\vartheta^{k}) &\leqslant (1+C_1 h)^{-k} E_{\mathrm{tot}}(\boldsymbol{u}^{0},\phi^{0},\vartheta^{0}) + \frac{C_2}{C_1}\left[1-(1+C_1 h)^{-k}\right]\notag \\ &\leqslant (1+C_1 h)^{-k} E_{\mathrm{tot}}(\boldsymbol{u}^{0},\phi^{0},\vartheta^{0}) + \frac{C_2}{C_1},\quad \forall\, k\in \mathbb{Z}^+. \label{discre-ineq-2} \end{align}\tag{73}\] Let \(E^N(t)\) be the piecewise linear interpolant of \(E_{\mathrm{tot}}(\boldsymbol{u}^k, \phi^k, \vartheta^k)\) at \(t_k= kh\) given by \[E^N(t)= \frac{(k+1)h-t}{h} E_{\mathrm{tot}}(\boldsymbol{u}^k, \phi^k, \vartheta^k) + \frac{t-kh}{h} E_{\mathrm{tot}}(\boldsymbol{u}^{k+1}, \phi^{k+1}, \vartheta^{k+1})\] for \(t\in [kh,(k+1)h)\), \(k\in \mathbb{N}\). For \(t\in (kh, (k+1)h)\), we also define the dissipation function \[D^N(t)= \frac{1}{4} \int_\Omega \left( \underline{\nu} | \nabla \boldsymbol{u}^{k+1} |^2 + \lambda_0 a \underline{m} | \nabla \mu^{k+1} |^2 + \underline{\kappa} | \nabla \vartheta^{k+1} |^2 \right)\,\mathrm{d}x.\] Then it follows from ?? that \[\begin{align} \frac{\mathrm{d}}{\mathrm{d} t} E^N(t) + D^N(t) & \leqslant - C_1 E_{\mathrm{tot}}(\boldsymbol{u}^{k+1},\phi^{k+1}, \vartheta^{k+1}) + C_2, \label{conti-energy-inequality} \end{align}\tag{74}\] for all \(t\in (kh,(k+1)h)\). Recalling that for all \(k\in \mathbb{N}\), \(E_{\mathrm{tot}}(\boldsymbol{u}^{k+1},\phi^{k+1}, \vartheta^{k+1})\geqslant -C_4\) for some positive constant \(C_4\) that depends only on \(\Omega\), we infer from 74 that \[\begin{align} E_{\mathrm{tot}}(\boldsymbol{u}^{j},\phi^{j}, \vartheta^{j}) + \int_{t_i}^{t_j} D^N(t)\,\mathrm{d}\tau \leqslant E_{\mathrm{tot}}(\boldsymbol{u}^{i},\phi^{i}, \vartheta^{i}) + (C_1C_4+C_2)(t_j-t_i), \label{discre-ineq-3} \end{align}\tag{75}\] for all \(t_i=ih\), \(t_j=jh\) with \(i, j\in \mathbb{N}\) and \(j\geqslant i\).

From 73 , 75 , the construction of the initial data \((\boldsymbol{u}^{0},\phi^{0}, \vartheta^{0})\) and the \(L^\infty\)-estimates for \(\phi^k\), \(\theta^k\), we can deduce that \[\begin{align} &\boldsymbol{u}^N\;\;\text{is bounded in}\;\;L^\infty(0,\infty; \boldsymbol{H}_\sigma)\cap L^2_{\mathrm{uloc}}([0,\infty);\boldsymbol{V}_{\sigma}), \\ &\phi^N \;\;\text{is bounded in}\;\;L^\infty(0,\infty;H^1(\Omega)),\quad \|\phi^N\|_{L^\infty(\Omega)}\leqslant 1\;\;\text{for all}\;t\geqslant 0, \\ &F(\phi^N) \;\;\text{is bounded in}\;L^\infty(0,\infty; L^1(\Omega)), \\ &\nabla \mu^N \;\;\text{is bounded in}\;\;L^2_{\mathrm{uloc}}([0,\infty);\boldsymbol{L}^2(\Omega)), \\ &\mathbf{J}^N \;\;\text{is bounded in}\;\;L^2_{\mathrm{uloc}}([0,\infty);\boldsymbol{L}^2(\Omega)), \\ & \vartheta^N\;\;\text{is bounded in}\;\;L^\infty(\Omega\times(0,\infty))\cap L^2_{\mathrm{uloc}}([0,\infty);H^1(\Omega)). \end{align}\] Recalling ?? , for \(\theta^N:=\vartheta^N+\Theta_\mathrm{b}\), we have \[\begin{align} & \theta^N\;\;\text{is bounded in}\;\;L^\infty(\Omega\times (0,\infty))\cap L^2_{\mathrm{uloc}}([0,\infty);H^1(\Omega)), \end{align}\] and moreover, \[\begin{align} & \min\Big\{\operatorname*{ess\,inf}_{\Omega}\theta^0,\,\min_{\partial\Omega}\theta_{\mathrm{b}}\Big\} \leqslant \theta^N(x) \leqslant \max\Big\{\operatorname*{ess\,sup}_{\Omega}\theta^0,\,\max_{\partial\Omega}\theta_{\mathrm{b}}\Big\},\quad \text{a.e. in}\;\overline{\Omega}. \label{maximumprinciple-theta-N} \end{align}\tag{76}\] Applying Lemma 4, we can further deduce that \[\begin{align} &\overline{\mu^N}\;\;\text{is bounded in}\;\;L^2_{\mathrm{uloc}}([0,\infty)), \\ & \mu^N \;\;\text{is bounded in}\;\;L^2_{\mathrm{uloc}}([0,\infty);H^1(\Omega)), \\ &\phi^N \;\;\text{is bounded in}\;\;L^4_{\mathrm{uloc}}([0,\infty);H^2(\Omega))\cap L^2_{\mathrm{uloc}}([0,\infty);W^{2,p}(\Omega)), \\ &F'(\phi^N) \;\;\text{is bounded in}\;\;L^2_{\mathrm{uloc}}([0,\infty);L^p(\Omega)), \end{align}\] for \(p\in [2,6]\) if \(d=3\), \(p\in[2,\infty)\) if \(d=2\).

In order to pass to the limit for all nonlinearities, we need to show strong convergence of \(\boldsymbol{u}^N\), \(\phi^N\), \(\vartheta^N\) (and also \(\theta^N\)), by investigating time derivatives of suitable approximate solutions.

Let \(\widetilde{\phi}^{N}\) be the piecewise linear interpolant of \(\phi^{N}(t^{k})\), where \(t^{k}=kh\), \(k\in \mathbb{N}\), that is, \(\widetilde{\phi}^{N}=\frac{1}{h}\chi _{\lbrack 0,h]}\ast _{t}\phi^{N}\), where the convolution is taken only with respect to the time variable \(t\). Similarly, we define \(\widetilde{\vartheta}^{N}=\frac{1}{h}\chi _{\lbrack 0,h]}\ast _{t}\vartheta^{N}\). Then it follows that \[\partial _{t}\widetilde{\phi }^{N}=\partial _{t,h}^{-}\phi ^{N},\quad \partial _{t}\widetilde{\vartheta }^{N}=\partial _{t,h}^{-}\vartheta^{N}, \quad \text{for almost all}\;t\in (0,\infty),\]and \[\| \widetilde{\phi }^{N}-\phi ^{N}\|_{(H^{1}(\Omega ))'} \leqslant h\| \partial _{t}\widetilde{\phi }^{N}\|_{(H^{1}(\Omega))'},\quad \| \widetilde{\vartheta }^{N}-\vartheta ^{N}\|_{H^{-1}(\Omega)} \leqslant h\| \partial _{t}\widetilde{\vartheta}^{N}\|_{H^{-1}(\Omega)}. \label{tilde-error}\tag{77}\] From 70 , 71 , we can obtain \[\begin{align} &\partial _{t,h}^{-}\phi ^{N} \;\;\text{is bounded in}\;\;L^2_{\mathrm{uloc}}([0,\infty);(H^1(\Omega))'), \\ &\partial _{t,h}^{-}\vartheta^{N} \;\;\text{is bounded in}\;\; L^2_{\mathrm{uloc}}([0,\infty); H^{-1}(\Omega)), \end{align}\] the first follows from the boundedness of \(\phi^N\boldsymbol{u}^N\), \(\nabla \mu^N\) in \(L^2_{\mathrm{uloc}}([0,\infty);\boldsymbol{L}^2(\Omega))\), while the second follows from the boundedness of \(\vartheta^N\boldsymbol{u}^N\) and \(\nabla \vartheta^N\) in \(L^2_{\mathrm{uloc}}([0,\infty);\boldsymbol{L}^2(\Omega))\). By definition, we easily find \[\begin{align} &\partial _{t}\widetilde{\phi }^{N} \;\;\text{is bounded in}\;\;L^2_{\mathrm{uloc}}([0,\infty);(H^1(\Omega))'), \\ &\partial _{t}\widetilde{\vartheta }^{N} \;\;\text{is bounded in}\;\; L^2_{\mathrm{uloc}}([0,\infty); H^{-1}(\Omega)). \end{align}\]

Next, let \(\widetilde{\rho \boldsymbol{u}}^N\) be the piecewise linear interpolant of \(\rho^N\boldsymbol{u}^N(t^k)\), where \(t^{k}=kh\), \(k\in \mathbb{N}\). We can check that \(\boldsymbol{P}(\widetilde{\rho \boldsymbol{u}}^N)\) is bounded in \(L^\frac{4}{3}_{\mathrm{uloc}}([0,\infty);\boldsymbol{V}_{\sigma})\), which is a consequence of the boundedness of \(\boldsymbol{u}^N\) in \(L^\infty(0,\infty;\boldsymbol{H}_\sigma)\cap L^2_{\mathrm{uloc}}([0,\infty); \boldsymbol{V}_{\sigma})\) and the boundedness of \(\phi^N\) in \(L^2_{\mathrm{uloc}}([0,\infty); W^{2,6}(\Omega))\cap L^\infty(\Omega\times (0,\infty))\). Recalling the estimates (cf. [42]) \[\begin{align} &\rho_h^N \boldsymbol{u}^N\otimes\boldsymbol{u}^N \;\;\text{is bounded in}\;\;L^2_{\mathrm{uloc}}([0,\infty);\boldsymbol{L}^\frac{3}{2}(\Omega)), \\ & \boldsymbol{u}^N\otimes \mathbf{J}^N \;\;\text{is bounded in}\;\;L^\frac{8}{7}_{\mathrm{uloc}}([0,\infty);\boldsymbol{L}^\frac{4}{3}(\Omega)), \\ & \mu^N\nabla \phi_h^N \;\;\text{is bounded in}\;\;L^2_{\mathrm{uloc}}([0,\infty);\boldsymbol{L}^\frac{3}{2}(\Omega)), \end{align}\] and the following fact \[\begin{align} \|(\vartheta_h^N+ \Theta_{\mathrm{b}})(\nabla \phi^N\otimes \nabla \phi^N)\|_{L^2(\Omega)} & \leqslant \|\theta_h^N\|_{L^\infty(\Omega)}\|\nabla \phi^N\|_{L^\infty(\Omega)}\| \nabla \phi^N\|\notag \\ &\leqslant C\|\theta_h^N\|_{L^\infty(\Omega)}\|\phi^N\|_{W^{2,4}(\Omega)}\|\phi^N\|_{H^1(\Omega)},\notag \end{align}\] which implies \[\begin{align} (\vartheta_h^N+ \Theta_{\mathrm{b}})(\nabla \phi^N\otimes \nabla \phi^N)\;\;\text{is bounded in}\;\;L^2_{\mathrm{uloc}}([0,\infty);\boldsymbol{L}^2(\Omega)), \end{align}\] then from 66 and 69 , we can deduce that \[\begin{align} \partial_t(\widetilde{\rho \boldsymbol{u}}^N)=\partial^{-}_{t,h}(\rho^N\boldsymbol{u}^N)\;\;\text{is bounded in}\;\;L^\frac{8}{7}_{\mathrm{uloc}}([0,\infty);\boldsymbol{W}^{-1,4}(\Omega)). \end{align}\]

Step 3. Passage to the limit as \(N\to \infty\). In Step 2, we have collected all the ingredients necessary to draw the conclusion of Theorem 1. In particular, all the estimates obtained above are independent of the parameters \(N\), \(h\) and \(k\). Hence, applying the same compactness argument as in [42], [54] with minor modifications, we can pass to the limit as \(N \to \infty\) (always understood in the sense of a convergent subsequence) to obtain a global weak solution \((\boldsymbol{u},\phi,\mu,\vartheta)\) (and thus \(\theta=\vartheta+\Theta_\mathrm{b}\)) of problem 1 –7 on \([0,\infty)\) with required regularity properties.

In the following, we will not repeat the limiting process as \(N\to \infty\) mentioned above. To finish the proof, we just sketch the convergence for the nonlinear term \(\theta_h^N(\nabla \phi^N\otimes \nabla \phi^N)\) due to the Marangoni effect. Below the convergence as \(N\to\infty\) will be understood in the sense of a subsequence, and the convergent subsequence will not be relabeled for simplicity.

With the uniform estimates obtained in Step 2, we recall that the following weak and strong convergence results can be obtained for \(\phi^N\) as \(N\to \infty\) (see [42]) \[\begin{align} & \phi^N \to \phi, \qquad \text{weakly star in}\;\;L^\infty(0,\infty;H^1(\Omega)),\\ & \phi^N \to \phi, \qquad \text{strongly in}\;\;L^2(0,T;H^1(\Omega)), \;\;\text{for any} \;T\in(0,\infty), \end{align}\] for some function \(\phi\in L^\infty(0,\infty;H^1(\Omega))\cap C_w([0,\infty);H^1(\Omega))\). Moreover, the boundedness of \(\phi^N\) in \(L^4_{\mathrm{uloc}}([0,\infty);H^2(\Omega))\cap L^2_{\mathrm{uloc}}([0,\infty);W^{2,6}(\Omega))\) also implies \[\begin{align} & \phi^N \to \phi, \qquad \text{weakly in}\;\;L^4(0,T;H^2(\Omega))\cap L^2(0,T;W^{2,6}(\Omega)), \;\;\text{for any} \;T\in(0,\infty), \end{align}\] with \(\phi\in L^4_{\mathrm{uloc}}([0,\infty);H^2(\Omega))\cap L^2_{\mathrm{uloc}}([0,\infty);W^{2,6}(\Omega))\). The linear dependence of \(\rho\) on \(\phi\) enables us to derive analogous results for \(\rho^N\). Similar arguments also apply to \(\vartheta^N\). First, we easily find that as \(N\to \infty\), it holds \[\begin{align} & \vartheta^N \to \vartheta, \qquad \text{weakly star in}\;\;L^\infty(\Omega \times (0,\infty) ),\\ & \vartheta^N \to \vartheta, \qquad \text{weakly in}\;\;L^2(0,T;H^1(\Omega)),\;\;\text{for any} \;T\in(0,\infty), \end{align}\] for some function \(\vartheta\in L^\infty(\Omega \times (0,\infty)) \cap L^2_{\mathrm{uloc}}([0,\infty);H^1(\Omega))\). Next, by the uniform estimates for \(\widetilde{\vartheta}^N\) and the Aubin–Lions compactness lemma, we get \[\begin{align} & \widetilde{\vartheta}^N \to \widetilde{\vartheta}, \qquad \text{strongly in}\;\;L^2(0,T;L^2(\Omega)), \;\;\text{for any} \;T\in(0,\infty), \end{align}\] for some function \(\widetilde{\vartheta}\in L^2_{\mathrm{uloc}}([0,\infty);H^1(\Omega)) \cap H^1_{\mathrm{uloc}}([0,\infty); H^{-1}(\Omega))\). Then it follows from the Lions–Magenes lemma (see e.g., [62], or [63]) that the limit \(\widetilde{\vartheta}\) satisfies \(\widetilde{\vartheta} \in C([0,\infty);L^2(\Omega))\cap L^\infty(0,\infty;L^2(\Omega))\). On the other hand, 77 implies that as \(N\to \infty\), \[\vartheta^N-\widetilde{\vartheta}^N \to 0 \qquad \text{strongly in}\;\;L^2(0,T;H^{-1}(\Omega)), \;\;\text{for any} \;T\in(0,\infty),\] since \(\partial _{t}\widetilde{\vartheta }^{N}\) is uniformly bounded in \(L^2_{\mathrm{uloc}}([0,\infty); H^{-1}(\Omega))\), and \(h \to 0\). As a result, we have \(\vartheta=\widetilde{\vartheta}\). From the fact \[\begin{align} \|\vartheta^N -\vartheta\| &\leqslant \| \vartheta^N-\widetilde{\vartheta}^N\|+\|\widetilde{\vartheta}^N - \widetilde{\vartheta}\|\\ &\leqslant \| \vartheta^N-\widetilde{\vartheta}^N\|_{H^1(\Omega)}^\frac{1}{2} \| \vartheta^N-\widetilde{\vartheta}^N\|_{H^{-1}(\Omega)}^\frac{1}{2}+\|\widetilde{\vartheta}^N - \widetilde{\vartheta}\|, \end{align}\] we arrive at \[\vartheta^N \to \vartheta \qquad \text{strongly in}\;\;L^2(0,T;L^2(\Omega)), \;\;\text{for any} \;T\in(0,\infty).\] By the \(L^\infty\)-boundedness of \(\vartheta^N\), \(\vartheta\) and the interpolation, for any \(q\in (2,\infty)\), we get \[\vartheta^N \to \vartheta \qquad \text{strongly in}\;\;L^q(0,T;L^q(\Omega)), \;\;\text{for any} \;T\in(0,\infty).\] Recalling the definition of \(\vartheta_h^N\), we also find that \[\vartheta_h^N \to \vartheta \qquad \text{strongly in}\;\;L^q(0,T;L^q(\Omega)), \;\;\text{for any} \;T\in(0,\infty).\]

Let \(T\in (0,\infty)\) be given. With the above convergence results, for any \(\boldsymbol{v}\in L^4(0,T;\boldsymbol{W}^{1,3}(\Omega))\), we consider the difference \[\begin{align} & \int_0^T \!\! \int_\Omega \big[\theta_h^N(\nabla \phi^N\otimes \nabla \phi^N)- \theta(\nabla \phi\otimes \nabla \phi)\big] : \nabla \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t \\ &\quad = \int_0^T\!\!\int_\Omega (\vartheta_h^N-\vartheta) (\nabla \phi^N\otimes \nabla \phi^N): \nabla \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t + \int_0^T\!\!\int_\Omega \big[\theta \nabla \phi^N\otimes (\nabla \phi^N-\nabla \phi)\big]: \nabla \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t\\ &\qquad + \int_0^T\!\!\int_\Omega \big[ \theta( \nabla \phi^N -\nabla \phi)\otimes \nabla \phi\big]: \nabla \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t. \end{align}\] The first term on the right-hand side can be estimated as follows \[\begin{align} &\left|\int_0^T\!\!\int_\Omega (\vartheta_h^N-\vartheta) (\nabla \phi^N\otimes \nabla \phi^N): \nabla \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t\right| \\ &\quad \leqslant \int_0^T \|\vartheta_h^N-\vartheta\|_{L^6(\Omega)} \|\nabla \phi^N\|_{L^4(\Omega)}^2\|\nabla \boldsymbol{v}\|_{L^3(\Omega)}\,\mathrm{d}t\\ &\quad \leqslant C \int_0^T \|\vartheta_h^N-\vartheta\|_{L^6(\Omega)} \|\phi^N\|_{H^2(\Omega)}^\frac{3}{2}\|\phi^N\|_{H^1(\Omega)}^\frac{1}{2}\|\nabla \boldsymbol{v}\|_{L^3(\Omega)}\,\mathrm{d}t \\ &\quad \leqslant \|\vartheta_h^N-\vartheta\|_{L^6(0,T;L^6(\Omega))} \|\phi^N\|_{L^4(0,T;H^2(\Omega))}^\frac{3}{2} \|\phi^N\|_{L^\infty(0,T;H^1(\Omega))}^\frac{1}{2} \|\nabla \boldsymbol{v}\|_{L^\frac{24}{11}(0,T;L^3(\Omega))} \\ &\quad \to 0\quad \text{as} \;N\to \infty. \end{align}\] Concerning the second term, we see that \[\begin{align} & \left| \int_0^T\!\!\int_\Omega \big[\theta \nabla \phi^N\otimes (\nabla \phi^N-\nabla \phi)\big]: \nabla \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t\right| \\ &\quad \leqslant \|\theta\|_{L^\infty(\Omega\times(0,T))} \int_0^T \|\nabla \phi^N\|_{L^6(\Omega)}\|\nabla \phi^N-\nabla \phi\|\|\nabla \boldsymbol{v}\|_{L^3(\Omega)}\,\mathrm{d}t\\ &\quad \leqslant C \|\theta\|_{L^\infty(\Omega\times(0,T))} \|\phi^N\|_{L^4(0,T;H^2(\Omega))}\|\phi^N-\phi\|_{L^2(0,T;H^1(\Omega))}\|\nabla \boldsymbol{v}\|_{L^4(0,T;L^3(\Omega))}\\ &\quad \to 0\quad \text{as} \;N\to \infty. \end{align}\] Finally, noticing that \(\nabla \phi^N\to \nabla \phi\) weakly in \(L^4(0,T; L^6(\Omega))\) and \[\begin{align} \|\theta\partial_i\phi\partial_j\boldsymbol{v}\|_{L^\frac{4}{3}(0,T;L^\frac{6}{5}(\Omega))} & \leqslant \|\theta\|_{L^\infty(\Omega\times(0,T))} \|\partial_i\phi\|_{L^4(0,T;L^2(\Omega))} \|\partial_j\boldsymbol{v}\|_{L^2(0,T;L^3(\Omega))}, \end{align}\] we obtain the convergence for the third term \[\begin{align} \int_0^T\!\!\int_\Omega \big[ \theta( \nabla \phi^N -\nabla \phi)\otimes \nabla \phi\big]: \nabla \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t \to 0\quad \text{as} \;N\to \infty. \end{align}\] Hence, for any test function \(\boldsymbol{v}\in L^4(0,T;\boldsymbol{W}^{1,3}(\Omega))\), it holds \[\begin{align} & \int_0^T \!\! \int_\Omega \big[\theta_h^N(\nabla \phi^N\otimes \nabla \phi^N)- \theta(\nabla \phi\otimes \nabla \phi)\big] : \nabla \boldsymbol{v}\,\mathrm{d}x\mathrm{d}t \to 0\quad \text{as} \;N\to \infty. \end{align}\]

Using the convergence results above, we can pass to the limit as \(N\to \infty\) in 66 –71 to obtain ?? –?? . In addition, following the argument in [42], we see that the initial data \((\boldsymbol{u}_0,\phi_0,\theta_0)\) can be attained.

The proof of Theorem 1 is complete. 0◻

4 Uniqueness in Two Dimensions↩︎

This section is devoted to the proof of Theorem 2. We first investigate an auxiliary problem for the convective heat equation with a pair of given data \((\boldsymbol{u}, \phi)\). By analyzing the regularity of solutions to the auxiliary problem, we establish the existence of global weak solutions to the two-dimensional problem 1 –7 with improved regularity properties for \(\theta\), provided that the initial temperature satisfies \(\theta_0 \in C^{\gamma}(\overline{\Omega}) \cap H_0^1(\Omega)\) for some \(\gamma \in (0,1)\). After that, we prove the uniqueness result in two dimensions under additional structural assumptions.

4.1 An auxiliary problem↩︎

Let \((\boldsymbol{u},\phi)\) be given functions. We consider the following initial boundary value problem of a convective heat equation: \[\begin{align} &\partial_t \theta + \boldsymbol{u}\cdot \nabla \theta = \mathrm{div} ( \kappa(\phi, \theta) \nabla \theta ), \quad \,\text{in}~ \Omega \times (0,\infty), \tag{78} \\ &\theta |_{\partial \Omega} = \theta_{\mathrm{b}},\qquad\qquad \qquad \qquad \qquad\, \text{on}\;\partial\Omega\times (0,\infty), \tag{79}\\ &\theta |_{t = 0} = \theta_0, \qquad \qquad \qquad \qquad \quad \;\;\; \text{in}\;\Omega. \tag{80} \end{align}\]

4.2 Proof of Theorem 2↩︎

We are ready to prove Theorem 2.

Step 1. Existence. Let \((\boldsymbol{u}^*, \phi^*, \mu^*, \theta^*)\) be a global weak solution to problem 1 –7 constructed in Theorem 1. Under the additional assumption \(\theta_0 \in C^{\gamma}(\overline{\Omega}) \cap H^1(\Omega)\) for some \(\gamma \in (0,1)\), \(\theta_0|_{\partial\Omega}=\theta_{\mathrm{b}} \in H^{\frac{3}{2}}(\partial \Omega)\), we consider the auxiliary problem 78 –80 with the given data \(\boldsymbol{u}=\boldsymbol{u}^*\), \(\phi=\phi^*\). According to Proposition 2-(ii), (iii), problem 78 –80 admits a unique global strong solution \(\theta\). Noticing that \(\theta^*\) can be regarded as a global weak solution of problem 78 –80 with the same structural data, then by the weak-strong uniqueness result Proposition 2-(iii), we have \(\theta^*(t)=\theta(t)\) for \(t\geqslant 0\).

Therefore, under the additional assumption on \(\theta_0\) mentioned above, \((\boldsymbol{u}^*, \phi^*, \mu^*, \theta^*)\) indeed gives a global weak solution to problem 1 –7 with improved regularity properties for the temperature as stated in Theorem 2.

Step 2. Uniqueness. In the case of unmatched densities \(\rho_1\neq \rho_2\), the uniqueness of weak solutions is beyond reach even in the absence of temperature coupling (recall Remark 6). On the other hand, if the mobility and thermal diffusivity depend both on the order parameter and on the temperature, additional regularity assumptions on \(\theta_1\) seem necessary to guarantee the uniqueness (see Remark 7), as will become apparent in the subsequent proof. Therefore, in the following, we focus on the special case such that \[\rho_1=\rho_2, \quad m(\phi, \theta) \equiv m(\phi), \quad \kappa(\phi, \theta) \equiv \kappa(\theta).\]

Let \((\boldsymbol{u}_i,\phi_i,\mu_i,\theta_i)\), \(i=1,2\), be two global weak solutions to problem 1 –7 constructed in Step 1 satisfying the same boundary and initial conditions. For simplicity, we denote their difference by \[\boldsymbol{u}= \boldsymbol{u}_1 - \boldsymbol{u}_2, \quad \phi = \phi_1 - \phi_2, \quad \mu = \mu_1 - \mu_2, \quad \theta = \theta_1 - \theta_2.\]

We first consider the equation satisfied by \(\boldsymbol{u}\). For every \(\boldsymbol{v}\in \boldsymbol{V}_{\sigma}\), it holds \[\begin{align} &\left\langle \partial_t \boldsymbol{u}, \boldsymbol{v}\right\rangle_{\boldsymbol{V}_{\sigma}} + \int_{\Omega} (\boldsymbol{u}_1 \cdot \nabla \boldsymbol{u}+ \boldsymbol{u}\cdot \nabla \boldsymbol{u}_2) \cdot \boldsymbol{v}\: \mathrm{d}x + \int_{\Omega} 2 \nu(\phi_1, \theta_1) D \boldsymbol{u}: D \boldsymbol{v}\: \mathrm{d}x\notag \\ &\qquad + \int_{\Omega} 2 (\nu(\phi_1, \theta_1) - \nu(\phi_2, \theta_2)) D \boldsymbol{u}_2: D \boldsymbol{v}\: \mathrm{d}x \notag \\ &\quad = \int_{\Omega} \big( \lambda(\theta_1) \nabla \phi_1 \otimes \nabla \phi_1 - \lambda(\theta_2) \nabla \phi_2 \otimes \nabla \phi_2 \big) : \nabla \boldsymbol{v}\: \mathrm{d}x + \int_{\Omega} (\boldsymbol{f}_{\mathrm{b}}(\phi_1, \theta_1) - \boldsymbol{f}_{\mathrm{b}}(\phi_2, \theta_2) ) \boldsymbol{v}\: \mathrm{d}x. \label{2D-dif-ua} \end{align}\tag{92}\] Taking \(\boldsymbol{v}= \boldsymbol{S}^{-1} \boldsymbol{u}\) in 92 yields \[\begin{align} \frac{1}{2} \frac{\: \mathrm{d}}{\: \mathrm{d}t} \| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|^2 + \int_{\Omega} 2 \nu(\phi_1, \theta_1) D \boldsymbol{u}: D (\boldsymbol{S}^{-1} \boldsymbol{u}) \: \mathrm{d}x = I_1 + I_2 + I_3 + I_4, \label{2D-dif-u} \end{align}\tag{93}\] where \[\begin{align} &I_1 = - \int_{\Omega} 2 (\nu(\phi_1, \theta_1) - \nu(\phi_2, \theta_2)) D \boldsymbol{u}_2:D (\boldsymbol{S}^{-1} \boldsymbol{u}) \: \mathrm{d}x, \\ &I_2 = - \big[ (\boldsymbol{u}_1 \otimes \boldsymbol{u}, \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) ) + (\boldsymbol{u}\otimes \boldsymbol{u}_2, \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) ) \big], \\ &I_3 = \int_{\Omega} \big( \lambda(\theta_1) \nabla \phi_1 \otimes \nabla \phi_1 - \lambda(\theta_2) \nabla \phi_2 \otimes \nabla \phi_2 \big) : \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \: \mathrm{d}x, \\ &I_4 = \int_{\Omega} (\boldsymbol{f}_{\mathrm{b}}(\phi_1, \theta_1) - \boldsymbol{f}_{\mathrm{b}}(\phi_2, \theta_2) ) \boldsymbol{S}^{-1} \boldsymbol{u}\: \mathrm{d}x. \end{align}\] By the definition of the Stokes operator, we have \(-\Delta(\boldsymbol{S}^{-1} \boldsymbol{u}) + \nabla \pi = \boldsymbol{u}\). Then it follows that \[\begin{align} &\int_{\Omega} 2 \nu(\phi_1, \theta_1) D \boldsymbol{u}: D (\boldsymbol{S}^{-1} \boldsymbol{u}) \: \mathrm{d}x\\ &\quad = -2 \int_{\Omega} \mathrm{div} ( \nu(\phi_1,\theta_1) D (\boldsymbol{S}^{-1} \boldsymbol{u}) ) \cdot \boldsymbol{u}\: \mathrm{d}x\\ &\quad = - \int_{\Omega} \nu(\phi_1,\theta_1) (\Delta \boldsymbol{S}^{-1} \boldsymbol{u}) \cdot \boldsymbol{u}\: \mathrm{d}x - 2 \int_{\Omega} \big[D (\boldsymbol{S}^{-1} \boldsymbol{u}) (\partial_1 \nu \nabla \phi_1 + \partial_2 \nu \nabla \theta_1)\big]\cdot \boldsymbol{u}\: \mathrm{d}x\\ &\quad = \int_{\Omega} \nu(\phi_1,\theta_1) | \boldsymbol{u}|^2 \: \mathrm{d}x - I_5 - I_6 - I_7, \end{align}\] where \[\begin{align} &I_5 = \int_{\Omega} \nu(\phi_1,\theta_1) \nabla \pi \cdot \boldsymbol{u}\: \mathrm{d}x, \\ &I_6 = 2 \int_{\Omega} \big[D (\boldsymbol{S}^{-1} \boldsymbol{u}) \partial_1 \nu \nabla \phi_1\big] \cdot \boldsymbol{u}\: \mathrm{d}x, \\ &I_7 = 2 \int_{\Omega} \big[ D (\boldsymbol{S}^{-1} \boldsymbol{u}) \partial_2 \nu \nabla \theta_1\big] \cdot \boldsymbol{u}\: \mathrm{d}x. \end{align}\] Therefore, we arrive at \[\begin{align} \frac{1}{2} \frac{\: \mathrm{d}}{\: \mathrm{d}t} \| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|^2 + \int_{\Omega} \nu(\phi_1,\theta_1) | \boldsymbol{u}|^2 \: \mathrm{d}x = \sum\limits_{k=1}^7 I_k. \end{align}\] We write \(I_1 = I_{1a} + I_{1b}\), where \[\begin{align} &I_{1a} = - \int_{\Omega} 2 (\nu(\phi_1, \theta_1) - \nu(\phi_1, \theta_2)) D \boldsymbol{u}_2:D (\boldsymbol{S}^{-1} \boldsymbol{u}) \: \mathrm{d}x, \\ &I_{1b} = - \int_{\Omega} 2 (\nu(\phi_1, \theta_2) - \nu(\phi_2, \theta_2)) D \boldsymbol{u}_2:D (\boldsymbol{S}^{-1} \boldsymbol{u}) \: \mathrm{d}x. \end{align}\] Using a similar reasoning as in [30], we have the following estimates \[\begin{align} &I_{1a} \leqslant \frac{\underline{\nu}}{20} \| \boldsymbol{u}\|^2 + \frac{\underline{\kappa}}{10} \| \nabla \theta \|^2 + C \| \nabla \boldsymbol{u}_2 \|^2 (\| \theta \|^2 + \| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|^2), \\ &I_2 \leqslant \frac{\underline{\nu}}{20} \| \boldsymbol{u}\|^2 + C ( \| \boldsymbol{u}_1 \|_{\boldsymbol{V}_{\sigma}}^2 + \| \boldsymbol{u}_2 \|_{\boldsymbol{V}_{\sigma}}^2 ) \| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|^2, \\ &I_3 \leqslant \frac{1}{40} \| \nabla \phi \|^2 + \frac{\underline{\kappa}}{10} \| \nabla \theta \|^2 + C ( \| \nabla \phi_1 \|_{L^{\infty}(\Omega)}^2 + \| \nabla \phi_2 \|_{L^{\infty}(\Omega)}^2 + \| \phi_1 \|_{H^2(\Omega)}^2 ) \| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|^2, \\ &I_4 \leqslant C \| \theta \|^2 + \| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|^2, \\ &I_7 \leqslant \frac{\underline{\nu}}{20} \| \boldsymbol{u}\|^2 + C \| \theta_1 \|_{H^2}^2 \| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|^2. \end{align}\] Let us recall the \(L^2\)-product estimate given in [66]: \[\| \phi D (\boldsymbol{S}^{-1} \boldsymbol{u}) \| \leqslant C \| \phi \|_{H^1(\Omega)} \| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \| \ln^{\frac{1}{2}}\left( e \frac{\| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|_{H^1(\Omega)}}{\| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|} \right).\] Using the Poincaré–Wirtinger inequality, Poincaré’s inequality and Young’s inequality, we can deduce that \[\begin{align} I_{1b} &\leqslant 2 \|\partial_1\nu\|_{L^\infty(\Omega)} \| D \boldsymbol{u}_2 \| \| \phi D (\boldsymbol{S}^{-1} \boldsymbol{u}) \| \\ &\leqslant C \| \nabla \boldsymbol{u}_2 \| \| \phi \|_{H^1(\Omega)} \| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \| \ln^{\frac{1}{2}}\left( e \frac{\| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|_{H^1(\Omega)}}{\| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u}) \|} \right) \\ &\leqslant \frac{1}{40} \| \nabla \phi \|^2 + C \| \nabla \boldsymbol{u}_2 \|^2 \mathcal{Y}(t) \ln \left( \frac{C}{\mathcal{Y}(t)} \right), \end{align}\] where \[\begin{align} \mathcal{Y}(t)= \| \nabla (\boldsymbol{S}^{-1} \boldsymbol{u})(t) \|^2 + \int_{\Omega} m(\phi_1(t)) | \nabla \mathcal{G}_{\phi_1} \phi (t)|^2 \: \mathrm{d}x + \| \theta(t)\|^2. \label{2D-dif-Y} \end{align}\tag{94}\] In the above estimate for \(I_{1b}\), we have used the fact that \(\mathcal{Y}(t) \in L^\infty(0,\infty)\) and the function \(s\ln\left(\frac{C}{s}\right)\) is increasing for \(s\in (0,L]\) if \(C\geqslant eL\). As far as the term \(I_5\) is concerned, using the Gagliardo–Nirenberg inequality and the \(L^4\)-estimate for the pressure \(\pi\) (see [67]), we have \[\begin{align} I_5 &= - \int_{\Omega} (\partial_1 \nu \nabla \phi_1 + \partial_2 \nu \nabla \theta_1 ) \cdot \boldsymbol{u}\pi \: \mathrm{d}x\\ &\leqslant ( \| \partial_1 \nu\|_{L^{\infty}(\Omega)} \| \nabla \phi_1 \|_{L^4(\Omega)} + \| \partial_2 \nu\|_{L^{\infty}(\Omega)} \| \nabla \theta_1 \|_{L^4(\Omega)} ) \| \boldsymbol{u}\| \| \pi \|_{L^4(\Omega)} \notag \\ &\leqslant C (\| \phi_1\|_{L^{\infty}(\Omega)}^{\frac{1}{2}} \| \phi_1 \|_{H^2(\Omega)}^{\frac{1}{2}} + \| \theta_1\|_{L^{\infty}(\Omega)}^{\frac{1}{2}} \| \theta_1 \|_{H^2(\Omega)}^{\frac{1}{2}} ) \| \boldsymbol{u}\|^{\frac{3}{2}} \| \nabla \boldsymbol{S}^{-1} \boldsymbol{u}\|^{\frac{1}{2}} \notag \\ &\leqslant \frac{\underline{\nu}}{20} \| \boldsymbol{u}\|^2 + C ( \| \phi_1\|_{H^2(\Omega)}^2 + \| \theta_1\|_{H^2(\Omega)}^2 ) \| \nabla \boldsymbol{S}^{-1} \boldsymbol{u}\|^2. \end{align}\] Finally, the term \(I_6\) can be estimated as in [66]: \[\begin{align} I_6 &\leqslant C \| \partial_1 \nu \|_{L^{\infty}} \| \nabla \phi_1 \|_{L^4(\Omega)} \| D \boldsymbol{S}^{-1} \boldsymbol{u}\|_{L^4(\Omega)} \| \boldsymbol{u}\| \\ &\leqslant C \| \phi_1 \|_{L^{\infty}(\Omega)}^{\frac{1}{2}} \| \phi_1 \|_{H^2(\Omega)}^{\frac{1}{2}}\| \nabla \boldsymbol{S}^{-1} \boldsymbol{u}\|^{\frac{1}{2}} \| \boldsymbol{u}\|^{\frac{3}{2}} \\ &\leqslant \frac{\underline{\nu}}{20} \| \boldsymbol{u}\|^2 + C \| \phi_1\|_{H^2(\Omega)}^2 \| \nabla \boldsymbol{S}^{-1} \boldsymbol{u}\|^2. \end{align}\]

Next, we consider the equations for \((\phi, \mu)\), which read as follows \[\begin{align} & \partial_t \phi + \boldsymbol{u}_1 \cdot \nabla \phi + \boldsymbol{u}\cdot \nabla \phi_2 = \mathrm{div} (m(\phi_1) \nabla \mu) + \mathrm{div} [(m(\phi_1) - m(\phi_2)) \nabla \mu_2], \tag{95} \\ &\mu = -\Delta \phi + W^{\prime}(\phi_1) - W^{\prime}(\phi_2). \tag{96} \end{align}\] Multiply 95 by \(\mathcal{G}_{\phi_1} \phi\) and then integrate over \(\Omega\). Recalling the definition of \(\mathcal{G}_{\phi_1}\) and the fact \(W^{\prime \prime} \geqslant -c_W\), we have \[\begin{align} -( \mathrm{div} (m(\phi_1) \nabla \mu), \mathcal{G}_{\phi_1} \phi) &= (\nabla \mu, m(\phi_1) \nabla \mathcal{G}_{\phi_1} \phi) \notag \\ &= (\mu, \phi) \notag \\ & \geqslant \| \nabla \phi \|^2 - c_W\|\phi\|^2\\ &\geqslant \frac{1}{2} \| \nabla \phi \|^2 - \frac{c_W^2}{2} \| \phi \|_{V_{(0)}^{\prime}}^2. \end{align}\] This implies that \[\begin{align} (\partial_t \phi , \mathcal{G}_{\phi_1} \phi) + \frac{1}{2} \| \nabla \phi \|^2 & \leqslant \frac{c_W^2}{2} \| \phi \|_{V_{(0)}^{\prime}}^2 + (\phi \boldsymbol{u}_1, \nabla \mathcal{G}_{\phi_1} \phi) + (\phi_2 \boldsymbol{u}, \nabla \mathcal{G}_{\phi_1} \phi) \notag \\ &\quad - \big((m(\phi_1) - m(\phi_2)) \nabla \mu_2 , \nabla \mathcal{G}_{\phi_1} \phi\big). \label{2D-dif-phi} \end{align}\tag{97}\] Using the Gagliardo–Nirenberg inequality, the Sobolev embedding theorem, Poincaré’s inequality and Young’s inequality, we get \[\begin{align} - ((m(\phi_1) - m(\phi_2)) \nabla \mu_2 , \nabla \mathcal{G}_{\phi_1} \phi) &\leqslant \|m'\|_{L^\infty(\Omega)} \| \phi \|_{L^4(\Omega)} \| \nabla \mu_2 \| \| \nabla \mathcal{G}_{\phi_1} \phi \|_{L^4(\Omega)} \\ &\leqslant C\| \nabla \mu_2 \| \| \nabla \mathcal{G}_{\phi_1} \phi \|^{\frac{1}{2}}\|\phi\| \| \nabla \phi \|^{\frac{1}{2}} \\ &\leqslant C\| \nabla \mu_2 \| \| \nabla \mathcal{G}_{\phi_1} \phi \| \| \nabla \phi \| \\ &\leqslant \frac{1}{40} \| \nabla \phi \|^2 + C\| \nabla \mu_2 \|^2 \| \nabla \mathcal{G}_{\phi_1} \phi \|^2, \end{align}\] and \[\begin{align} (\phi \boldsymbol{u}_1, \nabla \mathcal{G}_{\phi_1} \phi) &\leqslant \| \phi \|_{L^4(\Omega)} \| \boldsymbol{u}_1 \|_{L^4(\Omega)} \| \nabla \mathcal{G}_{\phi_1} \phi \| \\ &\leqslant \frac{1}{40} \| \nabla \phi \|^2 + C \| \boldsymbol{u}_1 \|_{L^4(\Omega)}^2 \| \nabla \mathcal{G}_{\phi_1} \phi \|^2, \\ (\phi_2 \boldsymbol{u}, \nabla \mathcal{G}_{\phi_1} \phi) &\leqslant \| \phi_2 \|_{L^{\infty}(\Omega)} \| \boldsymbol{u}\| \| \nabla \mathcal{G}_{\phi_1} \phi \| \\ &\leqslant \frac{\underline{\nu}}{20} \| \boldsymbol{u}\|^2 + C \| \nabla \mathcal{G}_{\phi_1} \phi \|^2. \end{align}\] The term \((\partial_t \phi , \mathcal{G}_{\phi_1} \phi)\) can be treated as in [50] such that \[\begin{align} (\partial_t \phi , \mathcal{G}_{\phi_1} \phi) &= \frac{\: \mathrm{d}}{\: \mathrm{d}t} \frac{1}{2} \int_{\Omega} m(\phi_1) | \nabla \mathcal{G}_{\phi_1} \phi |^2 \: \mathrm{d}x + \frac{1}{2} \int_{\Omega} \nabla \mathcal{G}\partial_t{\phi_1} \cdot m^{\prime \prime}(\phi_1) \nabla \phi_1 |\nabla \mathcal{G}_{\phi_1} \phi |^2 \: \mathrm{d}x\notag \\ &\quad + \int_{\Omega} \nabla \mathcal{G}\partial_t{\phi_1} \cdot m^{\prime}(\phi_1) (\nabla^2 \mathcal{G}_{\phi_1} \phi \nabla \mathcal{G}_{\phi_1} \phi) \: \mathrm{d}x, \label{transfer-phi-1} \end{align}\tag{98}\] with the following estimates (see [50]) \[\begin{align} \left|\frac{1}{2} \int_{\Omega} \nabla \mathcal{G}\partial_t{\phi_1} \cdot m^{\prime \prime}(\phi_1) \nabla \phi_1 |\nabla \mathcal{G}_{\phi_1} \phi |^2 \: \mathrm{d}x\right| \leqslant \frac{1}{40} \| \nabla \phi \|^2 + C \left( \| \nabla \mathcal{G}\partial_t{\phi_1} \|^2 + \| \phi_1 \|_{H^2(\Omega)}^4 \right) \| \nabla \mathcal{G}_{\phi_1} \phi \|^2, \end{align}\] \[\begin{align} \left|\int_{\Omega} \nabla \mathcal{G}\partial_t{\phi_1} \cdot m^{\prime}(\phi_1) (\nabla^2 \mathcal{G}_{\phi_1} \phi \nabla \mathcal{G}_{\phi_1} \phi) \: \mathrm{d}x\right| \leqslant \frac{1}{40} \| \nabla \phi \|^2 + C \left( \| \nabla \mathcal{G}\partial_t{\phi_1} \|^2 + \| \phi_1 \|_{H^2(\Omega)}^4 \right) \| \nabla \mathcal{G}_{\phi_1} \phi \|^2. \end{align}\]

Finally, we investigate the equation for \(\theta\), that is, \[\partial_t \theta + \boldsymbol{u}_1 \cdot \nabla \theta + \boldsymbol{u}\cdot \nabla \theta_2 = \mathrm{div} \big( ( \kappa(\theta_1) - \kappa(\theta_2) ) \nabla \theta_1 + \kappa(\theta_2) \nabla \theta \big).\label{2D-dif-the}\tag{99}\] Multiplying 99 by \(\theta\) and integrating over \(\Omega\), we have \[\begin{align} \frac{1}{2} \frac{\: \mathrm{d}}{\: \mathrm{d}t} \| \theta \|^2 + \int_{\Omega} \kappa(\theta_2) | \nabla \theta |^2 \: \mathrm{d}x = -\int_{\Omega} (\boldsymbol{u}\cdot \nabla \theta_2) \theta \: \mathrm{d}x -\int_{\Omega} ( \kappa(\theta_1) - \kappa(\theta_2) ) \nabla \theta_1 \cdot \nabla \theta \: \mathrm{d}x. \label{2D-dif-theb} \end{align}\tag{100}\] The second term on the right-hand side can be estimated in the same way as in 91 : \[\begin{align} -\int_{\Omega} ( \kappa(\theta_1) - \kappa(\theta_2) ) \nabla \theta_1 \cdot \nabla \theta \: \mathrm{d}x \leqslant \frac{\underline{\kappa}}{10} \| \nabla \theta \|^2 + C \| \theta_1 \|_{H^2(\Omega)}^2 \| \theta \|^2. \end{align}\] Concerning the first term on the right-hand side, we apply the Gagliardo–Nirenberg inequality, Hölder’s inequality and Young’s inequality and obtain \[\begin{align} -\int_{\Omega} (\boldsymbol{u}\cdot \nabla \theta_2) \theta \: \mathrm{d}x&\leqslant \| \boldsymbol{u}\| \| \nabla \theta_2 \|_{L^4(\Omega)} \| \theta \|_{L^4(\Omega)} \\ &\leqslant \frac{\underline{\nu}}{20} \| \boldsymbol{u}\|^2 + C \| \nabla \theta_2 \| \| \theta_2 \|_{H^2(\Omega)} \| \theta \| \| \nabla \theta \| \\ &\leqslant \frac{\underline{\nu}}{20} \| \boldsymbol{u}\|^2 + \frac{\underline{\kappa}}{10} \| \nabla \theta \|^2 + C \| \theta_2 \|_{H^2(\Omega)}^2 \| \theta \|^2. \end{align}\]

Combining all the above estimates, we deduce from 93 , 97 and 100 that \[\begin{align} &\frac{\: \mathrm{d}}{\: \mathrm{d}t} \mathcal{Y}(t) + \underline{\nu}\| \boldsymbol{u}\|^2 + \frac{1}{2} \| \nabla \phi \|^2 + \underline{\kappa} \| \nabla \theta \|^2 \leqslant C\mathcal{H}(t) \mathcal{Y}(t) \ln \left( \frac{C}{\mathcal{Y}(t)} \right), \label{2D-Diff} \end{align}\tag{101}\] where \(C>0\) is sufficiently large, \(\mathcal{Y}(t)\) is defined as in 94 and \[\begin{align} \mathcal{H}(t) & = 1 + \| \boldsymbol{u}_1(t) \|_{\boldsymbol{V}_{\sigma}}^2 + \| \boldsymbol{u}_2(t) \|_{\boldsymbol{V}_{\sigma}}^2 + \| \phi_1(t) \|_{W^{2,3}(\Omega)}^2 + \| \phi_2(t) \|_{W^{2,3}(\Omega)}^2 + \| \phi_1(t) \|_{H^2(\Omega)}^4 \\ &\quad + \|\nabla \mu_2(t) \|^2 + \| \nabla \mathcal{G}\partial_t{\phi_1}(t) \|^2 + \| \theta_1(t) \|_{H^2(\Omega)}^2 + \| \theta_2(t) \|_{H^2(\Omega)}^2. \end{align}\] Noticing that \(\mathcal{H}(t) \in L^1(0,T)\) for any \(T>0\), applying Osgood’s Lemma (see [68]), we can deduce from 101 that, if \(\mathcal{Y}(0)=0\), then \(\mathcal{Y}(t) =0\) for all \(t\in [0,T]\). Since \(T>0\) is arbitrary, we thus prove the uniqueness of global weak solutions on \([0,\infty)\).

The proof of Theorem 2 is complete. 0◻

Declarations↩︎

Conflict of interest. The authors have no competing interests to declare that are relevant to the content of this article.
Funding. The research of H. Wu was partially supported by Natural Science Foundation of Shanghai (Grant number 25ZR1401023).
Acknowledgments. The authors are grateful to the reviewer for several valuable comments. Part of the work was done during H. Wu’s participation in the Thematic Program on Free Boundary Problems at the Erwin Schrödinger International Institute for Mathematics and Physics (ESI), whose hospitality is gratefully acknowledged. H. Wu is a member of the Key Laboratory of Mathematics for Nonlinear Sciences (Fudan University), Ministry of Education of China.

References↩︎

---

1. School of Mathematical Sciences, Fudan University, Handan Road 220, Shanghai 200433, P. R. China. Email: chenlingxi@gmail.com↩︎

2. Corresponding author. School of Mathematical Sciences, Fudan University, Handan Road 220, Shanghai 200433, P. R. China. Email: haowufd@fudan.edu.cn↩︎