Two-Phase Flows with Bulk–Surface Interaction: Thermodynamically Consistent Navier–Stokes–Cahn–Hilliard Models with Dynamic Boundary Conditions

We derive a novel thermodynamically consistent Navier–Stokes–Cahn–Hilliard system with dynamic boundary conditions. This model describes the motion of viscous incompressible binary fluids with different densities. In contrast to previous models in the literature, our new model allows for surface diffusion, a variable contact angle between the diffuse interface and the boundary, and mass transfer between bulk and surface. In particular, this transfer of material is subject to a mass conservation law including both a bulk and a surface contribution. The derivation is carried out by means of local energy dissipation laws and the Lagrange multiplier approach. Next, in the case of fluids with matched densities, we show the existence of global weak solutions in two and three dimensions as well as the uniqueness of weak solutions in two dimensions.

. Difference between the Sharp Interface (left) and Diffuse Interface (right) methods the formulation of free boundary problems. In the latter, also called phase-field method, the interface is represented by a layer with finite thickness, corresponding to the level sets of the concentration function whose evolution is governed by a macroscopic equation. The main advantage of this approach is the Eulerian formulation of the phase-field equation, which avoids tracking the interface as in free boundary problems. Although conceptually opposites, these approaches are strictly related since a proper scaling of DI systems approximates SI models in the so-called sharp interface limit (see, e.g., [5]).
In the context of the DI theory, the cornerstone system for the motion of two viscous and incompressible fluids with matched (constant) densities is the model H, which was rigorously derived in [35]. The model H consists of the following Navier-Stokes-Cahn-Hilliard system Here, Q := Ω × (0, T ), where Ω ⊂ R d with d = 2, 3 denotes a bounded, smooth domain with boundary Γ and T > 0 is a given final time. In (1.1), ρ is the constant density of the mixture, v is the velocity of the mixture, p is the pressure, φ is the difference of the fluid concentrations, μ is the chemical potential. In addition, ν is the concentration depending viscosity, D is the symmetric gradient operator, ε is a positive parameter related to the thickness of the interface, m Ω is the (bulk) mobility function and F is a double-well potential. The fundamental assumption in deriving the model H is that the density of the mixture is constant, meaning that the two fluids have the same constant density ρ (matched densities). In the past two decades, several works investigated suitable Navier-Stokes-Cahn-Hilliard generalizations of the model H aiming to describe both incompressible mixtures with unmatched densities and compressible two-phase flows. Without claiming completeness, we refer the reader to [6,12,18,21,30,36,43,49,50]. Among them is the thermodynamically consistent model for viscous incompressible mixtures with unmatched densities proposed by Abels, Garcke and Grün (AGG model) in the seminal work [6], which reads as follows:  (1.4) where the total energy E is given by (1.6) denote the kinetic energy and the free energy in the bulk, respectively. We point out that in this model, the fluid concentrations on the surface do not affect the dynamics in the bulk. Concerning the mathematical analysis, the Cauchy problems corresponding to (1.1) and (1.2), each endowed with (1. 3), have been studied in [1,11,22,33] and in [3,4,7,31,32], respectively. We also refer the interested reader to [5] (see also [6]) for the connection between the above DI systems and the two-phase Navier-Stokes free boundary problem.
Even though the AGG model ((1.2),(1.3)) is already very advanced and is capable of describing the complex case where the two fluids have different specific densities, it still inherits some limitations from the underlying (convective) Cahn-Hilliard equation with homogeneous Neumann boundary conditions. The main limitations are: (L1) The boundary condition (1.3b) enforces the diffuse interface separating both fluids to always intersect the boundary Γ at a perfect angle of ninety degrees. Of course, this restrictive condition will not always be satisfied in concrete applications as the contact angle of the interface might deviate from ninety degrees and even change dynamically over the course of time. Moreover, as discussed in [47], the no-slip boundary condition (1.3a) is not well-suited for describing general moving contact line phenomena. This is because situations where the velocity field actually contributes to the motion of the contact line of the diffuse interface in the bulk at the boundary of the domain cannot be described. (L2) The boundary condition (1.3c) can be regarded as a no-flux boundary condition as it already implies J · n = 0 on Σ. Therefore, the system ((1.2),(1.3)) can merely describe the situation where the mass of both fluids in the bulk Ω is conserved but it is not capable of describing a transfer of material between bulk and surface which could, for instance, be caused by absorption or adsorption processes as well as chemical reactions taking place on the boundary (see, e.g., [38]).
Due to these well-known limitations of the Cahn-Hilliard model, and in order to describe short-range interactions between bulk and surface more precisely, physicists proposed that the total free energy should contain an additional contribution on the surface being also of Ginzburg-Landau type (see, e.g., [10]): and E bulk is defined as in (1.6). Here, ψ is an additional phase-field variable describing the material distribution on the surface, the parameter δ > 0 corresponds to the width of the diffuse interface on the surface and the constant κ ≥ 0 acts as a weight for surface diffusion effects. In most cases, ψ is just assumed to be the trace of the phase-field φ on the boundary, i.e., φ| Σ = ψ on Σ. (1.8) However, also more general relations between φ and ψ (so-called transmission conditions) have been investigated in the literature (see, e.g., [16,39]). The function G stands for an additional potential on the surface. If phase separating processes on the boundary are to be described, G is also chosen to be double-well shaped. Based on this free energy, several dynamic boundary conditions for the Cahn-Hilliard (CH) equation have been introduced in the literature. For instance, the Allen-Cahn type dynamic boundary condition (1.10) to replace (1.3b) was proposed in [37] and further analyzed in [15]. Here, the symbol Δ Γ denotes the Laplace-Beltrami operator (see, e.g., [20] for differential operators on surfaces). In recent times, dynamic boundary conditions which also exhibit a Cahn-Hilliard type structure have become very popular, especially since they allow for a better description of the transfer of material between bulk and surface (see, e.g., [38]). In these models, the boundary condition (1.3b) is replaced by a Cahn-Hilliard type equation on the surface, that is Here, θ denotes the chemical potential on the boundary, the function m Γ describes the mobility and the parameter κ ≥ 0 acts as a weight for surface diffusion effects. Furthermore, the boundary condition (1.3c) needs to be replaced by a suitable coupling condition for the bulk chemical potential μ and the surface chemical potential θ. Recently, for a parameter β = 0, the following choices have been considered: (1.11c) • If L = 0, (1.11c) is to be interpreted as the Dirichlet type coupling condition βθ = μ on Σ. This means that the chemical potentials μ and θ are assumed to always remain in chemical equilibrium. In this case, the dynamic boundary conditions (1.11) were proposed in [24] (for κ = 0) and in [34] (for general κ ≥ 0). Sufficiently regular solutions of the system ((1.8),(1.9),(1.11)) satisfy the mass conservation law for all t ∈ [0, T ]. In the case β > 0, the well-posedness of the system ((1.8),(1.9),(1.11)) was studied in [34] and its long-time behavior was investigated in [29,34]. • For L = ∞, the dynamic boundary conditions (1.11) were derived in [42] by means of an energetic variational approach based on Onsager's law and the least action principle. Here, due to (1.11c), the chemical potentials are not directly coupled. However, mechanical interactions of the materials between bulk and surface are still taken into account by the trace relation (1.8). Since (1.11c) with L = ∞ implies that the mass flux between bulk and surface is zero, the bulk mass and the surface mass are conserved separately. To be precise, sufficiently regular solutions of the system ((1.8),(1.9),(1.11)) satisfy the mass conservation law and the energy dissipation law (1.13) for all t ∈ [0, T ]. The well-posedness of the system ((1.8), (1.9), (1.11)) was established in [28,42] and its long-time behavior was studied in [42,45]. • In the case L ∈ (0, ∞), the boundary conditions (1.11) were proposed and analyzed in [38]. Here, the chemical potentials μ and θ are coupled by a Robin type boundary condition. The constant 1/L is related to the kinetic rate associated with adsorption/desorption processes or chemical reactions on the boundary. Sufficiently regular solutions of the system ((1.8),(1.9),(1.11)) satisfy the mass conservation law (1.12) and the energy dissipation law for all t ∈ [0, T ]. The weak and strong well-posedness of the system ((1.8),(1.9),(1.11)) was established in [38]. It was further shown that the case L ∈ (0, ∞) can be understood as an interpolation between L = 0 and L = ∞ as these cases are obtained as asymptotic limits on the level of strong solutions to the system ((1.8),(1.9),(1.11)) as the parameter L is sent to zero or to infinity, respectively. For the investigation of long-time dynamics, we refer to [29]. We further point out that a nonlocal variant of the system ((1.8),(1.9),(1.11)) was proposed and investigated in [41]. Recent reviews of results concerning the Cahn-Hilliard equation with dynamic boundary conditions can be found in [44,54].
In the context of Navier-Stokes-Cahn-Hilliard models for two-phase flows, in order to overcome the aforementioned limitation (L1), the authors in [47] proposed a variational derivation through Onsager's principle of maximal energy dissipation of a new class of generalized Navier slip boundary conditions for two-phase flows where ψ = φ| Σ as in (1.8), and Here, the subscript τ denotes the tangential component of a vector w, i.e., w τ = w − (w · n)n, ∇ Γ is the tangential gradient on Γ and Δ Γ denotes the Laplace-Beltrami operator on Γ (see, e.g., [20] for differential operators on surfaces). The phenomenological parameters β, κ and γ are positive, and G represents a surface double-well potential. The surface dynamic of the concentration ψ is described by (1.16b) that is a convective Allen-Cahn type equation on the surface. It can be understood as an adaptation of the dynamic boundary condition (1.10) (with δ = 1) to the situation of an additional volume-averaged velocity field. Notice that the boundary condition (1.16c) (together with either (1.1) or (1.2)) entails that Ω φ(t) dx is invariant (as for (1.3c)), whereas Γ ψ(t) dΓ is not a conserved quantity. The model H (1.1) endowed with (1.16) was first studied in [25]. The authors proved the existence of global weak solutions assuming that both F and G are polynomial-like functions or F is a singular potential and G is polynomial function. In the former case, they also showed the convergence of any weak solution to a stationary state. Later on, the global existence of weak solutions for (1.2) endowed with (1.16) has been achieved in [26]. More recently, the compressible Navier-Stokes-Cahn-Hilliard system supplemented with (1.16) has been studied in [14], and the Navier-Stokes-Allen-Cahn and the Navier-Stokes-Voigt-Allen-Cahn systems with (1.16) have been analyzed in [23]. For the stochastic model H with dynamic boundary conditions (1.16), we mention the recent results in [27]. It is worth pointing out that none of the aforementioned works established any uniqueness result for weak solutions in either two or three dimensions. In fact, the only uniqueness result established so far concerns quasi-strong solutions for the Navier-Stokes-Voigt-Allen-Cahn systems with (1.16) in [23].
In this paper, we derive and study a new class of thermodynamically consistent Navier-Stokes-Cahn-Hilliard systems with dynamic boundary conditions. Our model derivation is based on local mass balance laws in the bulk and on the surface in which the mass fluxes are still to be determined. Arguing similarly as in [6], after considering local energy dissipation laws, we apply the Lagrange multiplier approach to complete our model derivation by identifying the unknown mass fluxes on the macroscopic level. The resulting Navier-Stokes-Cahn-Hilliard (NSCH) model reads as follows: In this model, m Ω and m Γ are nonnegative functions respresenting the mobilities. The positive parameters ε and δ are related to the thickness of the diffuse interfaces in the bulk and on the boundary, respectively. The constant κ ≥ 0 acts as a weight for surface diffusion effects. Moreover, β ∈ R satisfies β |Ω| + |Γ| = 0, and L ∈ [0, ∞]. In contrast to the dynamic boundary conditions (1.16) which are of Allen-Cahn type, the dynamic boundary conditions in system (1.17) are of Cahn-Hilliard type. More precisely, they can be understood as a convective variant of the dynamic boundary conditions ((1.8), (1.11)) presented above. The sharp interface limit corresponding to the system (1.17) is to be studied in a subsequent work.
The system (1.17) is associated with the total energy where E kin is the kinetic energy introduced in (1.6) and E free = E bulk + E surf is the free energy defined in (1.7). Sufficiently regular solutions of system (1.17) satisfy the mass conservation law (1.19) and the energy dissipation law (1.21) as long as the solution exists. After presenting the model derivation of the system (1.17) in Sect. 2, we will study the case of fluids having matched densities, namely ρ ≡ρ 1 =ρ 2 . This means that the flux J vanishes. We rewrite the system (1.17) as follows The model (1.22) can be regarded as a generalization of the model H with dynamic boundary conditions of Cahn-Hilliard type. We will focus on the mathematical analysis of system (1.22) in the case L = 0 and β > 0 (see system (3.12) in Sect. 3). In Theorem 3.3, we first prove the existence of global weak solutions to system (3.12) with regular (polynomial-like) potentials F and G in bounded Lipschitz domains (cf. (A2)-(A5) below for the specific assumptions). Weak solutions satisfy the variational formulation of the problem (3.12) and a weak energy dissipation law (see (3.18)). A weak solution is obtained as the limit of smooth solutions to suitably regularized problems based on a Faedo-Galerkin scheme, where appropriate compactness properties are derived on the basis of uniform and global estimates gained from the energy dissipation law (1.20). However, an ad hoc construction of the solution ansatz is required since both (φ, ψ) and (μ, θ) need to be approximated in the same finite-dimensional subspace in order to show the validity of (1.20) for such regularized solutions (cf. the change of variables (3.22)). In comparison with the works in literature [23,25,26], our argument is rather simple and relies only on a one-level Faedo-Galerkin approximation. Furthermore, it is shown that (μ, θ) ∈ L 4 (0, T ; L 2 ), which is a novel regularity for NSCH systems (and even CH equations) with regular potentials (cf., for instance, [11,22]). Next, in Theorem 3.4, we then establish the uniqueness of the weak solutions in the two dimensional case. This is the first uniqueness result in the literature for a NSCH system with dynamic boundary conditions. It is important to underline that the nonlinear coupling in the NSCH system with dynamic boundary conditions (1.22) is much stronger than in the model H (1.1). Roughly speaking, this is due to the elastic (surface) term [ψ∇θ] τ in (1.22h), in addition to the classical capillarity (bulk) term −εdiv(∇φ ⊗ ∇φ) (also referred to as Korteweg force). More precisely, gradient terms arising from equivalent formulations of the Korteweg force (cf. (3.11)) disappear in bulk integrals in view of the incompressibility constraint div v = 0. Instead, surface integrals involving the product between a gradient function and v does not vanish since div τ v = 0 on Σ. For this reason, the term [ψ∇θ] τ cannot be handled as in [33] through dual estimates. Therefore, the appropriate functional to control the difference of two solutions corresponds to the total energy. This functional is capable of reproducing the same cancellation of the nonlinear terms as in the derivation of (1.20). In order to rigorously justify this fact, two additional assumptions need to be imposed: the domain Ω has a C 3 -boundary, which is needed to gain further spatial regularity of (φ, ψ) and to recover the relations of the bulk/surface chemical potentials (μ, θ) almost everywhere; the relation between bulk and surface potentials F (s) = βG(s) for all s ∈ R, which enables us to apply a particular chain formula (see Appendix A).

Model Derivation
In this section we will derive the novel Navier-Stokes-Cahn-Hilliard system with dynamic boundary conditions (1.17).

Considerations Based on Local Mass Balance Laws
We consider the time evolution of two fluids (indexed with i = 1, 2) in a (sufficiently smooth) domain where Γ = ∂Ω. Let ρ i : Q ∪ Σ → R, i = 1, 2, denote the mass densities of these two fluids. The local mass balance equations in the bulk read as (2.1) Here, J i are the mass fluxes. As we want to allow for a transfer of material between bulk and surface, the local mass balance equation on the boundary is given by In this relation, K i stand for the mass fluxes on the boundary. The right-hand side β J i describes the transfer of mass between bulk and surface. Here, the constant β ∈ R acts as a weight of this mass transfer.
the mass balances (2.1) and (2.2) can be rewritten as Let now φ i : Q → R and ψ i : Σ → R, i = 1, 2, denote the volume fractions of the two fluids in the bulk and on the boundary, respectively. Provided that each fluid has a constant densityρ i , i = 1, 2, the volume fractions can be identified as In particular, this naturally entails the trace relation Under the assumption that the excess volume is zero, we have We further set Furthermore, let v : Q ∪ Σ → R d denote the volume averaged velocity field associated with the two fluids. It can be expressed as We assume that both fluids cannot permeate the boundary Γ, which leads to the condition v n = v · n = 0 on Σ. (2.11) We further write to denote the mass fluxes relative to the volume averaged velocity. The mass balance equations (2.14) and (2.15) can thus be expressed as In this context, J i and K i , i = 1, 2, can be regarded as diffusive flow rates. We now define the total mass density as This leads to the relations Next, assuming a conservation of linear momentum with respect to the velocity field v, we have Here,T denotes the stress tensor that needs to be specified by constitutive assumptions. In the following, we write This allows us to define the objective tensor which is frame indifferent. Hence, (2.21) can be rewritten as We now multiply (2.14) by 1/ρ i , i = 1, 2, and we add the resulting equations. Recalling (2.7) and (2.9), this yields (2.24) which means that the volume averaged velocity field is divergence-free. Multiplying (2.15) by 1/ρ i , i = 1, 2, and adding the resulting equations, we infer We now definẽ We next multiply (2.14) by 1/ρ i , i = 1, 2. Subtracting the resulting equations, we deduce that Moreover, multiplying (2.14) by 1/ρ i , adding the resulting equations, and recalling (2.7) and (2.24), we conclude that Proceeding similarly with (2.15), and recalling (2.7) and (2.25), we obtain Using (2.26) and (2.29), we infer that This means that J and 1 2 (ρ 2 −ρ 1 )J φ can merely differ by an additive divergence-free function. As in [6], we thus assume that In the same fashion, we derive the identity From (2.7), (2.8) and (2.16), assuming that φ and ψ only attain values in [−1, 1], we further infer that the density ρ = ρ(φ) is given as and due to (2.6), the trace ρ(φ)| Σ can be expressed as

Local Energy Dissipation Laws
Proceeding as in [6], we introduce the following energy density in the bulk: The first summand on the right-hand side stands for the kinetic energy density whereas the second summand denotes the free energy density in the bulk. On the boundary, we introduce the additional energy density where g stands for the free energy density on the surface. We now consider an arbitrary test volume V (t) ⊂ Ω, t ∈ [0, T ], that is transported by the flow associated with v. Let ν denote the outer unit normal vector field of V (t). In an isothermal situation, the second law of thermodynamics leads to the dissipation inequality Here J e and K e are energy fluxes that will be specified later. As we consider a closed system, there is no transfer of energy over the boundary Γ and thus, the domain of the first integral in the second line is just ∂V (t) ∩ Ω instead of ∂V (t). We further point out that ∂ Γ (∂V (t) ∩ Γ) ⊆ Γ is to be understood as the relative boundary of the set ∂V (t) ∩ Γ within the submanifold Γ, and ν Γ stands for the corresponding outer unit normal vector field. Applying Gauß's divergence theorem on both integrals in the second line, and recalling that ν = n on ∂V (t) ∩ Γ, we reformulate (2.41) as Applying the Reynolds transport theorem, we find that Similarly, using the transport theorem for evolving hypersurfaces (see, e.g., [9, Theorem 32]), we obtain In particular, the above inequality holds true for all test volumes V (t) ⊂ Ω with ∂V (t) ∩ Γ = ∅. We thus infer the local dissipation law in the bulk, which reads as Let now α > 0 be arbitrary and suppose that V (t) ⊂ Ω is a test volume with |V (t)| being sufficiently small such that Then, we use (2.45) to infer that Since α > 0 and the test volume V (t) were arbitrary (except for the above restriction on V (t)), we conclude the following local dissipation law on the boundary:

Completion of the Model Derivation by the Lagrange Multiplier Approach
We now complete the model derivation by means of the Lagrange multiplier approach. We introduce the functions μ and θ that will be fixed in the subsequent approach. In view of the identities (2.37) and (2.38), the local energy dissipation laws (2.46) and (2.47) can be written as Here, the functions μ and θ can be understood as Lagrange multipliers. In the following, for brevity, we will just write ρ, f and g instead of ρ(φ), f (φ, ∇φ) and g(ψ, ∇ Γ ψ). By the definition of the energy densities e Ω and e Γ (see (2.39) and (2.40)), we rewrite (2.48) and (2.49) as By the chain rule, the derivatives ∂ t f and ∂ t g can be expressed as In the following, for any function h : Q → R, we use the notation to denote its material derivative. In the bulk, we proceed exactly as in [6] to reformulate (2.50) as To ensure (2.54), we now choose the chemical potential μ and the energy flux J e as By these choices, the first two lines of the right-hand side in (2.54) vanish. We further assume the mass flux J φ to be of Fick's type, that is is a nonnegative function representing the mobility. Hence, (2.54) reduces to We now define the tensor Here, the variable p denotes the pressure, and I denotes the identity matrix. S is the viscous stress tensor that corresponds to irreversible changes of the energy due to friction. For Newtonian fluids, S is usually chosen as is a nonnegative function representing the viscosity of the fluids, and Dv is the symmetric gradient of v. By this choice, (2.58) is satisfied since we obtain by means of the identity pI : ∇v = p div v = 0 in Q. We now assume that the energy density f is of Ginzburg-Landau type, that is Here, the parameter ε > 0 is related to the thickness of the diffuse interface that separates the two fluids, and F is a potential that usually exhibits a double-well structure. Hence, the chemical potential μ reads as and the total stress tensor T is given by Plugging (2.57) into (2.37), we further obtain the equation Eventually, recalling (2.22) and (2.63), we use (2.19) to derive the equation We now consider the local energy dissipation law (2.51) on the boundary. Recalling the formulas for J e (see (2.56)) and T (see (2.63)), we infer that due to the definition of the surface gradient and the boundary condition (2.11). In order to ensure that the inequality (2.67) is fulfilled, we choose the chemical potential θ, the mass flux K ψ , and the energy flux K e as follows: is a nonnegative function representing the mobility. This means that the first two lines of the right-hand side in (2.67) vanish. Moreover, since v n = 0 (see (2.11)), we have v = v τ and thus, Therefore, in order to ensure that (2.67) is satisfied, we make the constitutive assumption where γ = γ(ψ) is a nonnegative function. This equation can be regarded as an inhomogeneous Navier slip boundary condition. We now assume that the energy density g is also of Ginzburg-Landau type, that is Here, δ > 0 is related to the thickness of the diffuse interface on the boundary, κ ≥ 0 acts as a weight for surface diffusion effects, and G is a potential that usually exhibits a double-well structure. Hence, recalling the definition of the energy density f (see (2.61)), the chemical potential θ is given as Thanks to (2.68), we further have Hence, (2.67) reduces to The first two terms on the right-hand side are clearly nonpositive as γ(ψ) and m Γ (ψ) are nonnegative. This means that (2.76) is fulfilled if one of the following boundary conditions holds: Finally, by substituting (2.57) and (2.70) into (2.38), we obtain for prescribed initial data v 0 , φ 0 and ψ 0 , we have thus derived the system (1.17).

Notation and Preliminaries
We fix some notation and assumptions that are supposed to hold throughout the remainder of this paper. Notation.
(N1) The set of natural numbers excluding zero is denoted by N, and we write . We write Γ := ∂Ω. For any real numbers k ≥ 0 and p ∈ [1, ∞], the Lebesgue and Sobolev spaces for functions defined on Ω with values in R are denoted as L p (Ω) and W k,p (Ω). We write · L p (Ω) and · W k,p (Ω) to denote the standard norms on these spaces. If p = 2, we use the notation H k (Ω) = W k,2 (Ω). We point out that H 0 (Ω) can be identified with L 2 (Ω). Analogously, the Lebesgue and Sobolev spaces on Γ are denoted by L p (Γ) and W k,p (Γ) with corresponding norms · L p (Γ) and · W k,p (Γ) , respectively. In the case of vector fields on Ω with values in R n for some n ∈ N with n > 1, we use the notation L p (Ω), W k,p (Ω) and H k (Ω). For simplicity, their norms are denoted as in the scalar case by · L p (Ω) , · W k,p (Ω) and · H k (Ω) , respectively. Let I be a closed interval in R and X be a Banach space. The space C(I; X) denotes the set of continuous functions from I to X and, for k ∈ N, C k (I; X) denotes the space of k-times continuously differentiable functions from I to X. In particular, we simply write C(I) and C k (I) if X = R. Moreover, C w (I; X) denotes the space of functions mapping from I to X which are continuous on I with respect to the weak topology on X. This means that for any function f ∈ C w (I; X) and every Furthermore, for any real numbers k ≥ 0 and p ∈ [1, ∞], the Bochner spaces of functions defined on an interval I in R with values in X are denoted by L p (I; X) and W k,p (I; X). (N3) For any Banach space X, we write X to denote its dual space. The associated duality pairing of elements φ ∈ X and ζ ∈ X is denoted as φ, ζ X . If X is a Hilbert space, we write (·, ·) X to denote its inner product. (N4) For any bounded Lipschitz domain Ω, we define to denote the (generalized) spatial mean of u. Here, |Ω| denotes the d-dimensional Lebesgue measure of Ω. The spatial mean of a function v ∈ H 1 (Γ) (or v ∈ L 1 (Γ), respectively) is denoted as v Γ and defined analogously. (N5) For any bounded Lipschitz domain Ω and d ∈ {2, 3}, we introduce the space To prove the existence of weak solutions in the case of matched densities, we make the following assumptions: with Lipschitz boundary Γ = ∂Ω and a final time T > 0. We further use the notation (A2) The constants occurring in the system (1.22) satisfy ε, δ, β, κ > 0 and L = 0. Since the choice of δ, ε and κ has no impact on the mathematical analysis, we will simply set δ = ε = κ = 1 without loss of generality. (A3) The potentials F : R → [0, ∞) and G : R → [0, ∞) are twice continuously differentiable and there exist two exponents p and q satisfying as well as constants c F , c G ≥ 0 such that the second-order derivatives satisfy the growth conditions for all s ∈ R.
These assumptions already entail that there exist constants c F , c G , c F , c G ≥ 0 such that F , G , F and G satisfy the growth conditions (A5) The viscosity of the mixture ν : R → R and the friction parameter γ : R → R are continuous, bounded and uniformly positive. This means that there exist positive constants ν 0 , ν 1 , γ 0 , γ 1 such that for all s ∈ R, Remark 3.1. We point out that the polynomial double-well potential is a suitable choice for F and G as it satisfies (A3) with p = 4 and q = 4. However, singular potentials like the logarithmic Flory-Huggins potential or the double-obstacle potential are not admissible in this setting as they do not even satisfy (3.1)-(3.2).
Preliminaries. We next introduce several function spaces, inner products, norms and operators that will be used throughout this paper. (P1) For any real numbers k ≥ 0 and p ∈ [1, ∞], we set and we identify L 2 with H 0 . Note that H k is a Hilbert space with respect to the inner product and its induced norm · H k := (·, ·) endowed with the inner product (·, ·) D β := (·, ·) H 1 and its induced norm. The space D β is a Hilbert space. Moreover, we define the product for all (φ, ψ), (ζ, ξ) ∈ L 2 . By means of the Riesz representation theorem, this product can be extended to a duality pairing on D β × D β , which will also be denoted as ·, · D β . In particular, the spaces D β , L 2 , D β form a Gelfand triplet, and the operator norm on D β is given by for all (φ, ψ) ∈ D β . (P3) We further recall the following bulk-surface Poincaré inequalities: (P3.1) For any α > 0, there exists a constant C P depending only on α and Ω such that for all (φ, ψ) ∈ D α with α |Ω| φ Ω + |Γ| ψ Γ = 0. This bulk-surface Poincaré inequality is a special case of the one established in [40,Lemma A.1] (with the parameters therein being chosen as K = 0 and α, β > 0). (P3.2) There exists a constant C P depending only on Ω such that

Existence and Uniqueness of Weak Solutions
In this section, we consider the case of matched densities, that is ρ ≡ρ 1 =ρ 2 . This means we consider the system (1.22) supplemented with the boundary condition (1.22g) with L = 0 and β > 0. For the analysis, without loss of generality, we set ρ = κ = ε = δ = 1. We further employ the usual decomposition Hence, replacing the pressure p by Here, we recall that in the case of matched densities, the term J vanishes. Therefore, in the Navier-Stokes equation (3.12a) and in the corresponding boundary condition (3.12h), the terms related to J do not appear any more. The total energy associated with this system is The notion of weak solutions to system (3.12) is defined as follows: for all test functions w ∈ H 1 div (Ω), (ζ, ξ) ∈ D β , (η, ϑ) ∈ D 1 . (iv) The functions φ and ψ satisfy the mass conservation law (3.17) (v) The functions v, φ, ψ, μ and θ satisfy the weak energy dissipation law The existence of such a weak solution is ensured by the following theorem.

Theorem 3.3.
Suppose that the assumptions (A1)-(A5) hold and let v 0 ∈ L 2 div (Ω) and (φ 0 , ψ 0 ) ∈ D 1 be arbitrary. Then, there exists a weak solution of system (3.12) in the sense of Definition 3.2. If we additionally suppose that the domain Ω is of class C k with k ∈ {2, 3}, and in the case d = 3, we further assume that assumption (A3) holds with p < 6, then (φ, ψ) ∈ L 2 (0, T ; H k ) as well as a.e. in Q, (3.19a) In two dimensions, we show the uniqueness of the weak solutions under additional assumptions on the domain Ω, the viscosity ν, the mobilities m Ω and m Γ and the potentials F and G. (3.20) Then, the weak solution of system (1.22) given by Theorem 3.3 is unique. In addition, given two weak solu-

21)
holds for all t ∈ [0, T ], where , and the constant C > 0 depends only on Ω, the parameters of the system and the norms of the initial data.

Proofs
Proof of Theorem 3.3. We construct a weak solution to system (3.12) by discretizing the weak formulation (3.16) through a Faedo-Galerkin scheme.
In order to obtain suitable uniform estimates for the Galerkin approximation, it is necessary that the phase-fields (φ, ψ) and the chemical potentials (μ, θ) are approximated in the same finite-dimensional subspace. In turn this is only possible if the phase-fields and the chemical potentials are coupled by the same Dirichlet type boundary condition (cf. (3.12g)), which is not the case for general β > 0. However, by a change of variables, we can reformulate system (3.12) in such a way that the Dirichlet type boundary A. Giorgini, P. Knopf JMFM conditions for the redefined variables are exactly the same. To this end, we set α := √ β > 0, and we introduce the functions ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Ψ := α −1 ψ, i.e., ψ = αΨ, By this change of variables, system (3.12) is equivalent to In view of Definition 3.2, a weak solution (v, φ, Ψ, μ, Θ) of system (3.22) needs to have the regularity and satisfy the weak formulation which reads as for all test functions w ∈ H 1 div (Ω) and (ζ, ξ), (η, ϑ) ∈ D α . We now prove the existence of a weak solution of system (3.23) (in the aforementioned sense) and afterwards, we use the change of variables (3.22) to infer the existence of a weak solution to the original model (3.12). The weak formulation (3.25) can indeed be discretized by a suitable Faedo-Galerkin scheme.
Step 1: Discretization via a bulk-surface Faedo-Galerkin scheme. It has already been established in [40,Theorem 4.4] that the second-order elliptic eigenvalue problem has countably many eigenvalues {λ i } i∈N with λ i ≥ 0 for all i ∈ N, and the corresponding eigenfunctions {(ζ k , ξ k )} k∈N ⊂ D α can be chosen such that they form an orthonormal basis of L 2 . Here, we fix the eigenfunction associated with the first eigenvalue λ 1 = 0 as (ζ 1 , ξ 1 ) = 1 Moreover, in view of their construction (see [40]), the eigenfunctions (ζ k , ξ k ) with k ≥ 2 satisfy the mean value constraint For any k ∈ N, we introduce the finite dimensional subspace and we write Next, we consider the Stokes operator with Navier boundary conditions (see, e.g., [2, Appendix A]) where P is the Leray projection operator, with domain Since A is a positive and self-adjoint operator on L 2 div (Ω) with compact inverse, and there exists a sequence of countably many positive eigenvalues λ i i∈N and the corresponding eigenfunctions {w i } i∈N ⊂ H 1 div (Ω), which are determined by can be chosen in such a way that they form an orthonormal basis of L 2 div (Ω). For any k ∈ N, we define the finite dimensional subspace U k = span {w 1 , . . . , w k } ⊂ H 1 div (Ω), and the L 2 div (Ω)-orthogonal projection onto U k is denoted by P U k . Next, for any k ∈ N and t ∈ [0, T ], we make the ansatz a.e. in Ω × Γ.
(3.33c) A. Giorgini, P. Knopf JMFM Here, the scalar, time-dependent coefficients a k i , b k i and c k i , i = 1, . . . , k are assumed to be continuously differentiable functions that are still to be determined. They need to be designed in such a way that the discretized weak formulation holds for all w ∈ U k and all (ζ, ξ) ∈ Z k , supplemented with the initial conditions  ξ 1 ), . . . , (ζ k , ξ k ), we conclude that (a k , b k ) is determined by a system of 2k ordinary differential equations whose right-hand side depends continuously on a k , b k and c k . Due to (3.35) and (3.36), this ODE system is subject to the initial condition Moreover, testing (3.34c) with (ζ 1 , ξ 1 ), . . . , (ζ k , ξ k ), we infer that c k is explicitly given by a system of k algebraic equations whose right-hand side depends continuously on b k . We now replace the variable c k appearing in the right-hand side of the aforementioned ODE system to obtain a closed 2k-dimensional ODE system describing (a k , b k ) whose right-hand side depends continuously on (a k , b k ) . We can thus apply the Cauchy-Peano theorem to infer the existence of at least one local solution (a k , b k ) : to the corresponding initial value problem. Without loss of generality, we assume that T * k ≤ T and that (a k , b k ) is right-maximal, meaning that T * k is chosen as large as possible. Now, we can reconstruct c k : [0, T * k ) → R k by the aforementioned system of k algebraic equations. In view of (3.33), we thus obtain functions ; D α , which fulfill the discretized weak formulation (3.34) on the time interval [0, T * k ).
Step 2: Uniform estimates. We establish suitable estimates for each approximate solution (v k , φ k , ψ k , μ k , θ k ) which are uniform with respect to the index k. In particular, let T k < T * k be arbitrary. In the following, let C denote generic positive constants depending only on the initial data and the constants introduced in (A2)-(A5) including the final time T , but is independent of k and T k .
Furthermore, it follows from Strauss' lemma (see [52, In the two dimensional case, in view of the enhanced time integrability of ∂ t v, we even obtain v ∈ C([0, T ]; L 2 div (Ω)). In summary, we conclude that the regularity condition (3.24) is satisfied and thus, by the change of variables (3.22), Definition 3.2(i) is fulfilled.
To verify the weak energy dissipation law, let σ ∈ C ∞ 0 (0, T ) be an arbitrary non-negative test function. Multiplying (3.37) by σ and integrating on [0, T ], we obtain (3.78) By the strong convergence in (3.55), it easily follows that as k → ∞. Similarly, by (3.61), we have as k → ∞. By means of (3.57), we also infer (3.81) In addition, in light of (A3) and (3.59), Fatou's lemma implies Combining (3.79)-(3.82), we obtain Thanks to (3.68)-(3.71) and the weak lower-semicontinuity of norms with respect to weak convergence, another application of Fatou's lemma entails that Recalling the growth conditions on F and G (see (A3)), we use the initial conditions (3.35) and (3.36) as well as the convergence properties of the projections P U k and P Z k along with Lebesgue's general convergence theorem to infer

E v(t), φ(t), αΨ(t) σ(t) dt
which, in turn, implies that for almost all t ∈ [0, T ]. It remains to show that the energy inequality (3.87) actually holds true for all t ∈ [0, T ]. First, we notice that the both integral terms in (3.87) depend continuously on time. Recalling the growth assumptions A. Giorgini, P. Knopf JMFM on G , we further derive the estimate for all s, t ∈ [0, T ]. This implies that the function t → Γ G(αΨ(t)) dS is continuous on [0, T ]. Furthermore, recalling φ ∈ C([0, T ]; L 2 (Ω)) as well as (3.73), we deduce that the functions are lower semicontinuous on [0, T ]. Here, for the first three functions, we used the weak lower semicontinuity of the respective norms. For the fourth function, we applied Fatou's lemma. Combining all these properties, it follows that the weak solution (v, φ, Ψ, μ, Θ) satisfies the weak energy dissipation law (3.87) for all times t ∈ [0, T ]. Hence, reversing the change of variables (3.22), we conclude that Definition 3.2(v) is fulfilled.
We have thus shown that the quintuplet (v, φ, ψ, μ, θ) is a weak solution of system (3.12) in the sense of Definition 3.2.
Step 5: Higher regularity. We now assume that Ω is of class C 2 and, if d = 3, we further assume p < 6. Thanks to the regularity theory for second order eigenvalue problems in [40,Proposition 4.1], we deduce from the ansatz (3.33b) that (φ k , Ψ k ) ∈ L 2 (0, T ; D α ∩ H 2 ) for all k ∈ N. In view of (3.34c), the pair φ k (t), Ψ k (t) is a weak solution of the elliptic problem for almost all t ∈ [0, T ]. We now aim to prove uniform estimates that are independent of k by exploiting (3.52a) and (3.52b). For this reason, the symbol C will stand for a generic constant independent of k. We first observe that Recalling (A3) and using Sobolev's embedding theorem, we derive the estimate for almost all t ∈ [0, T ]. Let now > 0 be arbitrary. Without loss of generality, we assume p ∈ (4, 6). Using the Gagliardo-Nirenberg inequality and Young's inequality, we observe that (3.90) In particular, this yields f k (t), g k (t) ∈ L 2 for almost all t ∈ [0, T ]. Introducing the functions ψ k := αΨ k and θ k := α −1 Θ k for any k ∈ N, we infer that for almost all t ∈ [0, T ], the pair φ k (t), ψ k (t) ∈ L 2 (0, T ; D 1 ∩ H 2 ) is a weak solution of the elliptic system Hence, by regularity theory for elliptic problems with bulk-surface coupling (see [40,Theorem 3.3 for almost all t ∈ [0, T ]. Choosing sufficiently small, we can absorb the -dependent term on the righthand side by the left-hand side. Integrating the resulting estimate with respect to t from 0 to T , we derive that Recalling the convergence properties established in Step 4, we now pass to the limit k → ∞. By means of the Banach-Alaoglu theorem, we infer from (3.92) that φ, ψ ∈ L 2 (0, T ; H 2 ). We further obtain as k → ∞. Hence, by passing to the limit in (3.91), we eventually conclude that (3.19) is fulfilled in the strong sense. We now suppose that Ω is even of class C 3 and, if d = 3, we further assume p < 6. In the following, we merely proceed formally by directly considering the solution (φ, ψ) of (3.19). Nevertheless, the reasoning below can be rigorously justified by a suitable approximation argument.
Let > 0 be arbitrary and without loss of generality, we merely consider the case p ∈ [5,6). Applying the Gagliardo-Nirenberg inequality and the Agmon inequality, we find for almost all t ∈ [0, T ]. Proceeding similarly, we derive the estimate

95)
Proof of Theorem 3.4. We recall that the existence of a weak solution already follows from Theorem 3.3.
As the functions ν, m Ω and m Γ are assumed to be constant, we assume, without loss of generality, that ν ≡ 1, m Ω ≡ 1 and m Γ ≡ 1. In the following, we write C to denote generic positive constants depending only on Ω, the parameters of the system (3.12) and the norms of the initial data.

Conflict of interests
The authors do not have any financial or non-financial interests that are directly or indirectly related to the work submitted for publication.
Open Access. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

A. Appendix: Chain Rule Formula
Let κ and σ be two positive constants. We define the function Φ σ κ : Proof. It is easy to see that Φ σ κ is convex and proper. The lower semicontinuity of Φ σ κ immediately follows from the weak lower semicontinuity of the norm with respect to the weak convergence.