On a non-isothermal diffuse interface model for two-phase flows of incompressible fluids

We introduce a diffuse interface model describing the evolution of a mixture of two different viscous incompressible fluids of equal density. The main novelty of the present contribution consists in the fact that the effects of temperature on the flow are taken into account. In the mathematical model, the evolution of the macroscopic velocity is ruled by the Navier-Stokes system with temperature-dependent viscosity, while the order parameter representing the concentration of one of the components of the fluid is assumed to satisfy a convective Cahn-Hilliard equation. The effects of the temperature are prescribed by a suitable form of the heat equation. However, due to quadratic forcing terms, this equation is replaced, in the weak formulation, by an equality representing energy conservation complemented with a differential inequality describing production of entropy. The main advantage of introducing this notion of solution is that, while the thermodynamical consistency is preserved, at the same time the energy-entropy formulation is more tractable mathematically. Indeed, global-in-time existence for the initial-boundary value problem associated to the weak formulation of the model is proved by deriving suitable a-priori estimates and showing weak sequential stability of families of approximating solutions.


Introduction
We study a non-isothermal diffuse interface model for the flow of a mixture of two viscous incompressible Newtonian fluids of equal density in a bounded domain Ω ⊂ R 3 . In classical models the interface between the two fluids is assumed to be a 2−dimensional sufficiently smooth surface; in this case capillarity phenomena are related to contact angle conditions and a jump condition for the stress tensor across the interface. This classical description fails when some parts of the interface merge or reconnect (developing singularities) due to droplet formation or coalescence of several droplets. Indeed, despite the large amount of mathematical literature on free boundary problems related to fluids with a classical sharp interface, most papers are confined to the case of flows without singularities in the interface and so far there is no satisfactory existence theory of weak solutions for a two-phase flow of two viscous, incompressible, immiscible fluids with a classical sharp interface. Thus, in order to avoid analytical problems related to interface singularities, an alternative approach, based on diffuse interface models, can be used. In this setting, the classical sharp interface, represented by a lowerdimensional surface, is replaced by a thin interfacial region, whose "thickness" is described by a small parameter ε > 0. Therefore a partial mixing of the macroscopically immiscible fluids is allowed; in order to describe this phenomenon, a new variable ϕ is introduced. This quantity may represent the concentration difference or the concentration of one component of the fluid.
The original idea of diffuse interface model for fluids goes back to Hohenberg and Halperin [21] and it is usually referred with the name "H-model". Later, Gurtin et al. [20] gave a continuum mechanical derivation based on the concept of microforces. For a review of the development of diffuse interface models and their applications we refer to [1,3] and the references therein.
The present contribution is aimed at extending the H-model to a non-isothermal setting, developing a thermodynamically consistent theory. The system of partial differential equations resulting from this approach couples the incompressible Navier-Stokes system for the velocity u with a Cahn-Hilliard system with convection, where we also account for the effects of the (absolute) temperature θ. Namely, we consider the following equations: div u = 0, (1.1) u t + u · ∇ x u + ∇ x p = div S − εdiv(∇ x ϕ ⊗ ∇ x ϕ), S = ν(θ)Du, (1.2) ϕ t + u · ∇ x ϕ = ∆µ, (1.3)

4)
c V (θ)θ t + c V (θ)u · ∇ x θ + θ(ϕ t + u∇ϕ) − div(κ(θ)∇ x θ) = ν(θ)|Du| 2 + |∇ x µ| 2 , (1.5) in Ω × (0, T ), being Ω a bounded and sufficiently regular subset of R 3 and T > 0 a given final time, which may be arbitrarily large. Here p is the pressure, S = ν(θ)Du represents the dissipative part of the stress tensor, where Du = (∇ x u+ ∇ t x u)/2, and ν(θ) > 0 is the viscosity of the mixture. Moreover, c V (θ) stands for the specific heat, κ(θ) indicates the heat conductivity, ε > 0 is a (small) parameter related to the "thickness" of the interfacial region, and µ is an auxiliary variable (usually named chemical potential) which helps particularly for the statement of the weak formulation of the model. Finally, F (ϕ) is some suitable energy density whose expression is specified below. The capillarity forces due to surface tension are modeled by an extra-contribution ε∇ x ϕ ⊗ ∇ x ϕ in the global stress tensor appearing in the right-hand side of (1.1). System (1.1)-(1.5) will be closed by adding the initial conditions and suitable boundary conditions. Namely, the Cahn-Hilliard and temperature equations will be complemented by no-flux conditions, while the velocity u will be assumed to satisfy the socalled complete slip conditions. As we will see (cf. Subsection 2.6 for more details), these choices are crucial as we formulate the weak version of the model. Nevertheless, some other choices could be considered as well (for instance, the case of periodic boundary conditions can be treated similarly). It is also worth noting that, integrating (1.3) in space and using the no-flux condition together with (1.1) and the complete slip boundary condition for u, one gets back the mass conservation property (cf. (2.4)

below).
Isothermal versions of our model have been studied by several authors (see, e.g., [2,5,18,39,42] and references therein) and a mathematical theory can be now considered to be well-established. On the other hand, at least up to our knowledge, a non-isothermal model for two-phase fluids has been analyzed only in the reference [40], where a linearization of the internal energy balance is used in order to describe the evolution of the temperature. This permits the authors to get rid of the quadratic terms in the right hand side of (1.5) and of the coupling beween (1.4) and (1.5). On the other hand, the resulting model turns out to be thermodynamically consistent only in a neighbourhood of the equilibrium temperature.
In this contribution, we describe the non-isothermal evolution of the fluid by means of a model that keeps its thermodynamical consistency in a wide temperature range. Moreover, we can prove global-in-time existence for a suitable weak formulation of the associated initial-boundary value problem in three dimensions of space and without any magnitude restriction on the data. The analysis of non-isothermal problems in mathematical modelling of advanced materials is gaining more and more importance in the recent years. We refer for instance, without aiming at completeness, to [4,8,9,10,11,13,15,17,23,24,25,26,32,33,34,35,36,37,38], where temperature-dependent models are presented for describing the evolution of several types of substances, like elastic media, plastic materials (possibly with hysteresis, fatigue and damage), shape memory alloys, water-ice mixtures, and liquid crystals. The idea of replacing the heat equation with the energy and entropy balances in the weak formulation has been originally developed in [6,12] in the framework of heat conduction phenomena in fluids and in [14] in the case of solid-liquid phase transitions. It is worth observing that the related notion of weak solution, which will be introduced in full detail in Subsection 2.3, is consistent with the standard (strong) one. Actually, it is not difficult to prove that, at least for sufficiently smooth weak solutions, the total energy balance together with the entropy inequality imply the original form of the heat (or, more precisely, internal energy balance) equation (1.5). On the other hand, since this regularity in our case is not at all known (for instance due to the occurrence of the 3D Navier-Stokes system), this notion of solution turns out to be particularly useful because it allows us to prove a global in time existence result in 3D and at the same time it guarantees the thermodynamical consistency of the model. Better regularity properties are expected to hold for weak solutions in the 2D case. This will be the subject of a forthcoming paper, where we will also analyze the long-time dynamics of the model.
Plan of the paper. In the next Section 2 we provide a physical derivation of the model by following a variant of the general approach devised by Frémond in [16]. Namely, the equations of the system are obtained by imposing the balances of energy and entropy in terms of the free energy functional Ψ and of the pseudopotential of dissipation Φ and assuming standard consitutive relations. In the subsequent Section 3 we introduce the main assumptions on data, which permit us to state the weak formulation of the problem and the main existence theorem. The proof of this result occupies the remainder of the paper and is split into two steps: a-priori estimates, which are described in Section 4, and weak-sequential stability, which is proved in the last Section 5.

Derivation of the model
We suppose that a two-component fluid occupies a bounded spatial domain Ω ⊂ R 3 , with a sufficiently regular boundary Γ. We let n denote the outer normal unit vector to Γ. We denote by u = u(t, x) the associated velocity field in the Eulerian reference system. Moreover, we introduce the absolute temperature θ(t, x) and the order parameter ϕ(t, x), representing the concentration difference, or the concentration of one component, of the fluid. Furthermore, we denote as the material derivative of a generic function w, while w t (or also ∂ t w) stands for the partial derivative with respect to t. We set H := L 2 (Ω) and V := H 1 (Ω). We will often write H in place of H 3 , or V in place of V 3 , when vector-valued functions are considered. In particular, (·, ·) will stand for the usual scalar product both in H and in H 3 . For every f ∈ V ′ we indicate by f the spatial mean of f over Ω, i.e.
where ·, · denotes the duality pairing between V ′ and V and |Ω| stands for the Lebesgue measure of Ω. We note as H 0 , V 0 and V ′ 0 the closed subspaces of functions (or functionals) having zero mean value in H, V , and, respectively, in V ′ . Then, by the Poincaré-Wirtinger inequality, represents a norm on V 0 which is equivalent to the norm inherited from V . In particular · V0 is a Hilbert norm and we can introduce the associated Riesz isomorphism mapping J : Moreover, if u is as above, then for all v ∈ V 0 . Finally, we can identify H 0 with H ′ 0 by means of the scalar product of H so to obtain the Hilbert triplet V 0 ⊂ H 0 ⊂ V ′ 0 , where inclusions are continuous and dense. In particular, if z ∈ V and v ∈ V 0 , it is easy to see that (2.2)

Free-energy and pseudopotential of dissipation
We would like to apply here the general approach proposed in the monograph [16,Chapters 2,3] in order to build our diffuse interface model for incompressible fluids with conserved order parameter. We start by introducing, in agreement with the basic principles of classical Thermodynamics, the free-energy and the pseudopotential of dissipation. To this aim, we specify the set of the state variables, describing the actual configuration of the material: Correspondingly, the set of the dissipative variables, whose evolution describes the way along which the system tends to dissipate energy, is given by where Du := ∇ x u + ∇ t x u 2 denotes the symmetric gradient of u.
Motivated by the Ginzburg-Landau theory for phase transitions, we choose the free energy density ψ and the free energy functional Ψ in the form where ε is a positive constant related to the interface thickness. The function F in (2.3) penalizes the deviation of the length |ϕ| from its natural value 1; generally, F is assumed to be a sum of a dominating convex (and possibly non smooth) part and a smooth non-convex perturbation of controlled growth. Typical example are the standard double-well potential F (ϕ) = (|ϕ| 2 −1) 2 and the so-called logarithmic potential Notice that in our case we will assume F to be a double-well potential and this hypothesis turns out to be essential in our analysis (cf. Sec. 4). Moreover, f represents the purely caloric part of the free energy and is linked to the specific heat c V (θ) = Q ′ (θ) by the relation . In what follows we will assume c V (θ) ∼ θ δ with δ ∈ (1/2, 1). While the condition δ < 1 is quite natural (in particular it ensures the concavity of the entropy as function of θ), the bound from below δ > 1/2 is crucial for the purpose of analyzing our PDE system (cf. Sec. 4 below). The evolution of the system is characterized by a second functional Φ, called pseudopotential of dissipation, assumed to be nonnegative and convex with respect to the dissipative variables. In order to present the explicit expression of Φ, however, we have to impose in some way the mass conservation constraint. To this aim, we first decompose ϕ = ϕ 0 + m 0 , where m 0 represents the mean value of the initial datum ϕ 0 . Then, the conservation of mass corresponds to prescribe ϕ 0 to take its values in H 0 during the whole evolution of the system.
Let us also note as V n the subspace of H 1 (Ω; R 3 ) consisting of the functions u such that u · n = 0 on Γ. Then, given a time-dependent family of divergence-free vector fields u(t, ·) ∈ V n and a scalar field ϕ = ϕ(t, x) satisfying the mass conservation constraint ϕ(t, ·) = ϕ(0, ·) for a.e. t ∈ (0, T ), we have In other words, Dϕ Dt has zero spatial mean. Therefore, if ϕ is so smooth to satisfy Dϕ Dt ∈ V ′ 0 a.e. in time, we can define and, consequently, µ 0 ∈ V 0 . Hence, we can set where the "local component" φ of the "dissipation density" is given by Here, ·, · denotes the duality between V ′ 0 and V 0 , ν = ν(θ) > 0 is the viscosity coefficient and κ = κ(θ) > 0 represents the heat conductivity. The incompressibility of the fluid is formally enforced by I 0 , i.e., the indicator function of {0} (given by I 0 = 0 if div u = 0 and +∞ otherwise). The last term in (2.6), which accounts for mass conservation, is nonstandard and requires some words of explanation. Basically, it corresponds to a (squared) V ′ 0 -norm of the material derivative of ϕ and, hence, depends on the dissipative variable Dϕ Dt in a nonlocal way. As will be seen below, this gives rise to a convective Cahn-Hilliard dynamics, in contrast with the Allen-Cahn dynamics that would result from the choice of the H-norm (cf. [16] for more details).
The functionals Φ and Ψ are assumed to be defined for all sets of variables E and δE for which they make sense. In other words, the finiteness, say, of Ψ determines the class of admissible state variables. So, in particular, if a time-dependent set of variables is given such that, a.e. in (0, T ), Ψ and Φ are finite, ϕ 0 (t, ·) ∈ V 0 , and u(t, ·) ∈ V n , then u is divergence-free and, by (2.4), Dϕ Dt has zero mean value.
We also note that, whenever Dϕ Dt ∈ H 0 , then, by elliptic regularity, it follows Moreover, we can equivalently rewrite the pseudopotential Φ as follows: Indeed, integrating by parts in space and using the definition of J (2.2), it turns out that: An alternative strategy to derive the (isothermal) Cahn-Hilliard equation starting from a balance of the so-called "microscopic motions" can be found in [19] and could be extended to the present case of binary fluids. However, while the analysis of the isothermal Cahn-Hilliard equation has been deeply investigated starting from the pioneering paper [7] and up to the most recent contributions and Navier-Stokes systems (cf., e.g., [2,18,20]), at least up to our knowledge a rigorous derivation of a thermodynamically consistent model including also the internal energy balance equation is still lacking (cf., e.g. [31] for a thermodynamically consistent model coupling the Cahn-Hilliard system with a singular heat equation). Our method is intended to fill this gap in the more intricated case where the effects of the macroscopic velocity u are also taken into account.

Constitutive relations
We start by introducing the energy density B and the energy flux vector H, both assumed to be the sum of their non-dissipative and dissipative components, namely, Moreover, we set H d ≡ 0. Here and in the sequel for simplicity we will use the symbol ∂ not only for partial derivatives, but also for variational derivatives of possibly non-convex funtionals. Relation (2.11) defines B d as the first variation (more precisely, subdifferential) of Φ with respect to Dϕ Dt in the space H 0 . To see that it coincides, indeed, with J −1 Dϕ Dt , we observe that, for any v ∈ H 0 , there holds The dissipative components of the heat and entropy fluxes (denoted respectively by q d and Q d ) are whereas the nondissipative components q nd and Q nd will be determined later on in such a way that the second law of Thermodynamics is satisfied (cf. (2.31) below). Of course, we assume that In what follows we also ask that κ(θ) . This choice is mainly motivated by mathematical reasons (actually, it guarantees some additional integrability of θ). A physical justification for it is provided, e.g., in [41]. Also the stress tensor σ is decomposed into the dissipative component and the non dissipative part σ nd to be determined below. The entropy of the system is given by and, finally, the internal energy e reads where Q(θ) = f (θ) − θf ′ (θ) represents, physically speaking, the antiderivative of the specific heat c V .

Field equations
In accordance with Newton's second law, the balance of momentum reads where g is a given external force. The balance of internal energy takes the form where we recall that the internal energy flux is decomposed as q = q d + q nd . Moreover, we notice that, on the right hand side of (2. 19), there appears a new (with respect to the standard theory of [16]) term N . This contribution is aimed at balancing the nonlocal dependence of the last term in the pseudopotential of dissipation Φ (cf. (2.6)) with respect to the dissipative variable Dϕ Dt , which is one of the peculiarities of Cahn-Hilliard systems (cf. also [26] for similar techniques applied to different nonlocal contributions). The expression of N (in terms of the dissipative variables δE) will be obtained below in such a way to comply with the second law of Thermodynamics. Since the main role of N is to model the nonlocal interactions between points inside Ω, it will result from our computations that Ω N (x) dx = 0, in agreement with natural expectations.
Finally, the equation ruling the evolution of the order parameter ϕ can be derived from the principle of virtual powers. Indeed, following the general theory developed in [16,Chap. 2 However, in the present setting the above relation is not completely rigorous, since it does not properly incorporate the boundary conditions and the mass conservation constraint. So, to be more precise, we first have to rewrite the expression of Ψ in the form where we used the decomposition ϕ = ϕ 0 + m 0 already introduced in Section 2.1.
In order to impose this constraint mathematically, we restate (2.20) as a generalized gradient flow problem in the space H 0 as follows: Let us note that asking for ϕ 0 to lie in the domain of the differential δ H0,ϕ 0 Ψ means that there exists a (unique) function z ∈ H 0 such that δ H0,ϕ 0 Ψ(ϕ 0 ) can be represented by z in the scalar product of H 0 (i.e., of H). In this way, (2.22) incorporates both the homogeneous Neumann boundary conditions for ϕ and the mass conservation property. Moreover, it is immediately seen that such a function z must have the expression Combining (2.22) and (2.23) with (2.11), we then get Applying the distributional Laplace operator to both hand sides and noting that −∆J −1 v = v for any v ∈ H 0 (cf. (2.2)), we then arrive at the system where the auxiliary variable µ is introduced mainly for mathematical convenience and takes the name of chemical potential. In fact, we may note that µ 0 = −J −1 ( Dϕ Dt ) = µ − µ (cf. (2.24), (2.5)). The non-dissipative components of the stress σ nd and of the flux q nd , as well as the "nonlocality compensation" term N , are determined by means of (2.19) and the constitutive relations derived above, in order for the second law of Thermodynamics to be satisfied. Indeed, computing De Dt from (2.17) by means of the standard Helmholtz relations Moreover, by (2.23)-(2.26), we have (pointwise) To deduce the expressions for the non-dissipative components of the stress σ nd and of the flux q nd as well as that of N , we impose validity of the Clausius-Duhem inequality in the form where Q denotes the entropy flux and it is linked to the flux q by the relation θQ = q. Recalling (2.10)-(2.14) and developing the left hand side, we get Then, in order to obtain the non-negativity of the right hand side (cf. (2.31)), we can assume, e.g., the following constitutive relations With these choices, we get Ω N (x) dx = 0, as expected. Moreover, a straighforward computation shows that the internal energy balance (2.19) can be rewritten as Notice that the dissipation terms on the right hand side are in perfect agreement with the expression (2.8) of the pseudopotential of dissipation Φ. Indeed, as already mentioned, one has µ 0 = µ − µ due to (2.25)-(2.26).

Strong formulation
On account of the derivation sketched above, we can now write the PDE system representing the strong formulation of our model. Firstly, collecting (2.15), (2.18) and (2.32), we obtain that the evolution of the velocity u is ruled by the Navier-Stokes system, given by Incompressibility: div u = 0; (2.34) Conservation of momentum: where p is the pressure, and Regarding the evolution of ϕ, (2.25)-(2.26) give rise to the following Cahn-Hilliard equation: (2.37) Finally, we rewrite (2.33) as the Internal energy balance: It is worth observing that, neglecting the temperature, equations (2.34-2.38) reduce to the model derived in [20] by means of a different method.

Balances for total energy and entropy
A key point in the statement of the weak formulation of our model consists in replacing the "heat" equation (2.38) with the balances of total energy and of entropy. These relations are, indeed, mathematically more tractable, but they keep all the main features of the problem (in particular, the first and second laws of Thermodynamics are still respected). We give here a formal derivation of these relations, beginning with the total energy balance. To obtain it, we start multiplying (2.35) by u, which gives (2.39) Next, we multiply the first (2.37) by µ, the second by ϕ t , and take the difference. After standard manipulations, we obtain We now substitute into the last term the expression of µ given by the second (2.37): Then, we can take the sum of (2.38), (2.39) and (2.42) to obtain the Total energy balance: and the heat flux Next, multiplying (2.38) by 1/θ, we obtain the Entropy equation: We anticipate that the function Λ is well defined in view of the assumptions on c V stated in (3.2) below. Moreover, it is worth noting that (2.46) is an equality at this level, but it will turn to an inequality (cf. (3.20) below) in the framework of the rigorous definition of weak solution that will be introduced later on (cf. Def. 1 and Theorem 1). Of course, this phenomenon is due to the quadratic terms on the right hand side, which do not behave well with respect to weak limits. However, it is worth noticing that, in case we could prove that there exist a smooth solution to the weak formulation of the model, then for that solution it would be possible to recover the entropy equality (2.46). Moreover, (2.46) is equivalent to (2.19) in that setting. In this sense, the weak formulation analyzed in the next section turns out to be compatible with the strong formulation (1.1-1.5) given at the beginning, at least when sufficient regularity holds.

Initial and boundary conditions
In order to get a well-posed problem, we have to specify suitable initial and boundary conditions. Compatibly with the physical derivation, we will essentially assume that the system is insulated from the exterior. This leads to taking the following no mass flux (through the boundary) condition: where we recall that n is the external normal. Next, we assume that This position prescribes a "contact angle" of π/2 between the diffuse interface and the boundary of the domain. Moreover, we take no-flux boundary conditions for the temperature: q · n |Γ = 0. The first condition states that the normal component of the boundary velocity is zero (so, the fluid cannot exit from Ω, but it can move tangentially to the boundary). The second position prescribes that there is no external contribution to the viscous stress. This, in a sense, excludes friction effects with the boundary. The above choice has also mathematical implications. Indeed, following the lines of [6], we will use it in order estimate the pressure term appearing in the total energy balance (2.43). It is worth noting that an estimate of the pressure can be reached also in the case when Ω is the unit torus and periodic boundary conditions are taken for all unknowns. In particular, our results could be extended to that situation with trivial modifications. For other types of boundary conditions, integrability of the pressure is, instead, an open issue. We remark once more that, thanks to (2.48) and the first (2.51), integrating (2.37) in space and time, we get mass conservation:

52)
This feature is characteristic of Cahn-Hilliard-type models and is in agreement with the underlying physics. Finally, the system is complemented by the initial conditions

Main result 3.1 Assumptions on coefficients and data
Before formulating the main result of the paper, we list here the hypotheses imposed on the constitutive functions. First of all, just for the sake of simplicity, we take ε = 1 and g = 0. The case of a nonconstant forcing term could be treated, indeed, with trivial modifications. Next, we assume that F (ϕ) is the classical double-well potential, namely More general expressions of F ′ having cubic growth at ∞ may be admissible as well, but we prefer to keep right from the beginning the expression (3.1) in order not to overburden the presentation. We assume that the thermal conductivity, the specific heat and the viscosity of the mixture depend on θ in the following way: for all θ ≥ 0, some 0 < ν < ν, and some β > 0, δ > 0 complying with the following restrictions: In view of (3.2), we can compute The ansatz δ < 1 comes from the assumption of physical consistency (actually, it implies that the thermal component Λ of the entropy is concave, as prescribed by Thermodynamics); the other limitations mainly have a mathematical motivation and are needed in order to obtain the necessary regularity to pass to the limit. Notice however that a power-like behavior for the heat-conductivity is typical of several types of fluids (cf., e.g., [41]). For brevity we also set p β,δ := β + 2 3 (δ + 1). (3.6) This exponent will be needed in the a-priori estimates derived below (cf. for instance (4.26) and (5.23)). We conclude by specifying our hypotheses on the initial data: Here and below, L 2 div indicates the space of divergence-free L 2 functions.

Weak formulation
First of all, we rewrite the momentum equation (2.18), with g = 0, in the more explicit form: This permits us to introduce the notion of weak solution to our model problem:

Definition 1.
A weak solution to the non-isothermal diffuse interface model for two-phase flows of fluids is a quadruplet (u, ϕ, µ, θ) satisfying the incompressibility condition div u = 0 a.e. in (0, T ) × Ω, the weak momentum balance
It is worth noting that (3.9) incorporates both the incompressibility constraint and the initial condition u(0, ·) = u 0 ; moreover, it accounts for the complete-slip conditions (2.51). The first equation (3.10) of the Cahn-Hilliard system is in weak form and also accounts for the no-flux condition (2.48), while we will be able to prove sufficient regularity on ϕ in order for (3.11) to hold pointwise (with the no-flux condition (2.49) in the sense of traces). To get (3.12), we tested (2.43) by ξ, integrated by parts in time and used the Cahn-Hilliard system (2.37). More precisely, we wrote

Main existence theorem
Our main result reads as follows: Under the assumptions stated in Subsection 3.1, the non-isothermal diffuse interface model for two-phase flows of fluids admits at least a weak solution, in the sense of Definition 1, in the following regularity class: θ > 0 a.e. in (0, T ) × Ω, log θ ∈ L 2 (0, T ; H 1 (Ω)), (3.19) δ and β being specified in (3.3). Moreover, this solution complies with the following weak form of the entropy production inequality: (3.20) holding for any ξ ∈ C ∞ 0 ([0, T ) × Ω), ξ ≥ 0, and where we have set Remark 1. Let us notice that in case we could prove existence of a sufficiently smooth weak solution (in particular, regular enough in order to integrate back by parts the terms in (3.12)), then it would be possible to show that such a solution also satisfies the "standard" form of the heat equation (1.5).
Equivalently, the entropy inequality (2.46) would hold as an equality in that case. Hence, the current notion of weak solution turns out to be compatible both with Thermodynamics and also with the "strong" one.

A priori bounds
The remainder of the paper is devoted to the proof of Theorem 1. We start by briefly sketching our strategy. In this section, we will prove some formal a-priori estimates holding for a hypothetical quadruple (u, ϕ, µ, θ) solving the "strong" formulation of the model stated in Subsection 2.3. Actually, these estimates will mainly follow as direct consequences of the Total energy balance (2.43) and Entropy inequality (2.46). Of course, to make this procedure fully rigorous, one should rather work on a proper regularization or approximation of the strong system and prove that it admits at least one solution being sufficiently smooth in order to comply with the estimates. However, the system stated in Subsection 2.3 is rather complex and the related approximation, on the one hand, would be particularly long and technical, and, on the other hand, should not present particular difficulties, or novelties, from the analytical viewpoint. Indeed, the nonisothermal Navier-Stokes system given by (2.35) and (2.38) can be treated along the lines developed in the monograph [12], while the regularization of the Cahn-Hilliard equation (2.37) is completely standard since no singular or nonsmooth terms are involved. For all these reasons, we decided to skip the details of this argument and rather proceed formally. We just quote the papers [2,13,14,15,27,28,29] where more details of possible approximations of related systems are given. In Section 5, having the a-priori estimates at disposal, we will then prove that any sequence (u n , ϕ n , µ n , θ n ) complying with the bounds uniformly in n admits at least one limit point (u, ϕ, µ, θ) which solves the weak formulation of the system (i.e., satisfies the conditions stated in Definition 1). This procedure, which will be referred to as "weak sequential stability" of families of solutions, can be seen as a simplified version of the compactness argument that one should use to remove some form of regularization or approximation.

Energy estimates
Integrating the total energy balance (2.43), we deduce the following a priori estimates:

Temperature estimates
We now integrate the temperature equation (2.38) in space and time, with the aim of getting a L 2bound for the two quadratic terms on the right hand side. Note that, at this level (or, to be more precise, in the approximation we decided to skip), it is crucial to have the strong relation (2.38) (or an approximated version of it "containing" the same information) at disposal. In other words, this procedure cannot be reproduced by workly directly on weak solutions since the total energy balance (3.12) alone is not sufficient. That said, we need to control the terms on the left hand side of (2.38). First of all, (4.1) brings Moreover, by boundary conditions and incompressibility, Finally, using the first equation of (2.37), we are able to deduce and we need to control the last two terms. Actually, the first one is absorbed by the last term on the right hand side of (2.38), while the second one is controlled thanks to (4.9). Hence, from the right hand side of (2.38) we can "read" the a priori estimates Let us note that here it is essential to have a cubic growth in 3D for the potential F ′ . Now we choose ξ = ϕ in (3.10), we test (3.11) by ∆ϕ, and take the sum. Some terms cancel out due to (2.34) and our choice of boundary conditions. Hence, using in particular (4.5), (4.15) and (4.16) for treating the other terms, it is not difficult to get the additional regularity on ϕ: ||∆ϕ|| L 2 ((0,T )×Ω) ≤ c ⇒ ||ϕ|| L 2 (0,T ;H 2 (Ω)) ≤ c, (4.17) where the classical regularity theorems for elliptic equations have also been used. Now, we would like to show that To this aim, we use classical interpolation inequalities. Indeed, ϕ ∈ L ∞ (0, T ; H 1 (Ω)) ∩ L 2 (0, T ; H 2 (Ω)) ֒→ L 2/ϑ (0, T ; H 1+ϑ (Ω)) ֒→ L 2/ϑ (0, T ; L 6 1−2ϑ (Ω)), provided that ϑ ∈ (0, 1/2), whence On the other hand At this point, combining (4.19) and (4.20) we need to find ϑ < 1/2 such that ϑ + ϑ 2 = 1 2 and 1 − 2ϑ 3 + 3 − 2ϑ 6 = 1 2 and this leads to ϑ = 1/3. For this value of ϑ we then have F ′′ (ϕ) ∈ L 3 (0, T ; L 9 (Ω)), ∇ x ϕ ∈ L 6 (0, T ; L 18/7 (Ω; R 3 )), whence (4.18).

Weak sequential stability
In this section, we assume to have a sequence (u n , ϕ n , µ n , θ n ) of solutions satisfying (a proper approximation of) the strong system of Subsection 2.3. Then, by virtue of the argument developed in the previous part, we can assume that this family complies with the proved a priori bounds uniformly with respect to n. Our aim is showing, by weak compactness arguments, that at least a subsequence converges in a suitable way to a weak solution to our problem (i.e. to a limit quadruple (u, ϕ, µ, θ) satisfying the statement given in Definition 1). Actually, to further simplify the notation, we intend that all the convergence relations appearing in the following are to be considered up to the extraction of (not relabelled) subsequences. That said, collecting the bounds proved before, we have where we recall that β and δ fufill assumption (3.3). Next, we need an estimate for the pressure p n . To achieve it, we will follow the approach devised in [6] for dealing with the Navier-Stokes-Fourier system. Referring to that paper for more details, a formal way to get integrability of p n consists in computing the divergence of (3.8), which gives rise to ∆p n = div div S n − u n ⊗ u n − ∇ x ϕ n ⊗ ∇ x ϕ n . (5.6) In other words, p n solves, at least formally, some kind of elliptic problem. Hence, an estimate for it can be proved by relying on suitable regularity theorems. To be more precise, (5.6) has to be interpreted in a "very weak" sense. Namely, p n turns out to satisfy the integral identity for any test function ξ ∈ H 2 (Ω) with ∇ x ξ · n| Γ = 0. The above formulation incorporates the boundary conditions, though in a way which is not completely obvious. In particular, one uses in an essential way the complete slip conditions (2.51) for u and the no-flux conditions for the other variables. Indeed, the Neumann condition (2.49) allows us to deal with the extra stress −∇ x ϕ n ⊗ ∇ x ϕ n . To see this at least formally, one multiplies (2.35) by ∇ x ξ. Then, thanks to (2.49) and (2.51), it is possible to integrate by parts and get (5.7) without the occurrence of additional boundary terms. This procedure can be made rigorous by following closely the lines of the approximation argument described in [6,Sec. 4], to which we refer the reader for more details.
That said, we want to apply suitable elliptic regularity theorems to (5.6) (or, more precisely, to its weak formulation (5.7)). To this aim, we first establish some bounds on the terms appearing on the right hand side.
On the other hand, thanks to (4.23), we have Then, collecting (5.9)-(5.12) and using elliptic regularity in (5.7) (again, the precise details are given in [6]), we deduce p n → p weakly in L Consequently, a comparison of terms in (3.8) gives where X is a suitable Sobolev space of negative order. By the Aubin-Lions lemma, (5.1) and (5.9), we then get u n → u strongly in L 3 + ((0, T ) × Ω; R 3 ). (5.14) Next, we observe that u n ·∇ x ϕ n is bounded in L 2 ((0, T )× Ω), due to (4.2) and (4.22). Hence, recalling (5.3) and comparing terms in the first (2.37) we infer Combining this with (5.2) and applying once more the Aubin-Lions lemma, we get for all σ > 0. Now, let us assume that, for every n ∈ N, the approximate solution (u n , ϕ n , µ n , θ n ) fulfills the strong system stated in Subsec. 2.4 (actually, its hypothetical approximation). Then, multiplying by suitable test functions and integrating by parts, it readily follows that (u n , ϕ n , µ n , θ n ) also complies with Definition 1. In particular, starting from the "strong" formulation (2.43) of the total energy balance (assumed to hold at the n-level), accounting for (3.14) and performing suitable integrations by parts, it is not difficult to deduce its weak counterpart ∂ t 1 2 |u n | 2 + e n + div 1 2 u n |u n | 2 + e n u n + div(p n u n ) − ∆κ(θ n ) − div(S n u n ) − ∆µ 2 n + div (u n · ∇ x ϕ n )∇ x ϕ n + div(∇ x ∇ x ϕ n ∇ x µ n ) − div div(∇ x µ n ⊗ ∇ x ϕ n ) = 0, (5.16) which holds at least in the sense of distributions, and where e n is defined as Let us note that (5.16) reduces to (3.12) after testing it by ξ ∈ C ∞ 0 ([0, T ) × Ω). It is worth remarking that, in the approximation, (u n , ϕ n , µ n , θ n ) also satisfies the strong form (2.46) of the entropy relation.
Then, we can see what happens as we let n ր ∞. Actually, on account of the properties proved above, it is a standard matter to see that both the Cahn-Hilliard system (3.10)-(3.11) and the momentum equation (3.9) pass to the desired limits. Hence, to conclude the proof, we have to prove (3.12) and (3.20) in the limit. This issue is a bit more delicate and is dealt with in the next two subsections.
Actually, the limit of |u n | 2 is identified as |u| 2 thanks to the strong convergence (5.14), while the limit e of e n (cf. (5.17)) still needs to be identified in terms of θ and ϕ. To this aim, we need strong convergence of θ n in L p ((0, T ) × Ω) for some p and, to achieve it, we use a monotonicity argument like in the paper [15]. Actually, from (5.15) we infer , strongly, say, in L 1 + ((0, T ) × Ω).
In view of the above discussion, all terms in the first row of (5.16) pass to the desired limits. On the other hand, the last three terms in the second row can be managed thanks to (5.20-5.22). Finally, combining (5.15) with (5.23) and comparing terms in (3.11), we obtain that µ n → µ = −∆ϕ + F ′ (ϕ) − θ strongly in L 2 ((0, T ) × Ω). (5.24) Hence, ∆µ 2 n → ∆µ 2 at least in the sense of distributions. This allows us to take the limit of (5.16), which immediately reduces to (3.12) after testing by ξ ∈ C ∞ 0 ([0, T ) × Ω) and integrating by parts.

Proof of the entropy inequality
To conclude the proof of Theorem 1, we need to prove the entropy production inequality (3.20). As noted above, we can assume that a stronger relation, i.e. (2.46), holds at the n-level and we aim at taking its (supremum) limit as n → ∞. In particular, this will give rise to the ≥ sign in (3.20). For the reader's convenience, we start by reporting the statement of a useful lower semicontinuity result due to A.D. Ioffe [22]: f (x, ·, ·) is lower semicontinuous on R n × R m for every x ∈ O, (5.25) f (x, u, ·) is convex on R m for every (x, u) ∈ O × R n . (5.26) Let also (u k , v k ), (u, v) : O → R n × R m be measurable functions such that Then, To begin, we test (2.46) at the n-level by a nonnegative test function ξ ∈ C ∞ 0 ([0, T ) × Ω), as specified in the statement of Theorem 1, and we integrate by parts. What we get is exactly the n-version of (3.20) (with the equal sign). In particular, integration by parts gives rise to the function h(θ n ) (cf. (3.21)). Thanks to (4.10) and (5.23) (cf. also Remark 2 below), it is not difficult to see that h(θ n ) → h(θ), say, strongly in L 1 + ((0, T ) × Ω), (5.27) whence in particular To deal with the other terms on the left hand side of (3.20) n , we observe that, due to (5.23) and (3.5), Λ(θ n ) → Λ(θ) strongly in L p ((0, T ) × Ω) for all p ∈ 1, p β,δ δ , whence, recalling (5.14), Λ(θ n )u n → Λ(θ)u, say, strongly in L 1 ((0, T ) × Ω; R 3 ).

Remark 2.
In order for our estimates to make sense, we implicitly assumed in the course of the proof the temperature θ to be (almost everywhere) positive. This fact is used in several estimates which, otherwise, would not make sense. Positivity of θ n should be shown, indeed, at the n-level, i.e., for the hypothetical regularized problem which we decided not to detail here. Actually, for regularized solutions, which are usually smoother, it is often possible to prove some stronger property (like θ n (x) ≥ c n > 0 a.e.), by applying a suitable maximum principle. We cannot give here a proof of this fact, since this would require to provide the details of the regularization. However, we can at least show that, if θ n is almost everywhere positive, and satisfies the estimates given in Section 4, then positivity is preserved in the limit. To see this, we first notice that, by (4.10), ∇ x log θ n L 2 ((0,T )×Ω;R 3 ) ≤ c. (5.29) Then, using that the mean value of ϕ n is conserved, i.e., ϕ n (t) = ϕ n (0) for every t ∈ (0, T ) and integrating (2.46) both in space and in time, we readily obtain Ω Λ(θ n (t)) ≥ Ω Λ(θ 0 ) for a.e. t ∈ (0, T ) or, equivalently, In particular, θ > 0 almost everywhere also in the limit.