Well-posedness of a bulk-surface convective Cahn–Hilliard system with dynamic boundary conditions

We consider a general class of bulk-surface convective Cahn–Hilliard systems with dynamic boundary conditions. In contrast to classical Neumann boundary conditions, the dynamic boundary conditions of Cahn–Hilliard type allow for dynamic changes of the contact angle between the diffuse interface and the boundary, a convection-induced motion of the contact line as well as absorption of material by the boundary. The coupling conditions for bulk and surface quantities involve parameters K, L ∈ [0 , ∞ ], whose choice declares whether these conditions are of Dirichlet, Robin or Neumann type. We ﬁrst prove the existence of a weak solution to our model in the case K, L ∈ (0 , ∞ ) by means of a Faedo–Galerkin approach. For all other cases, the existence of a weak solution is then shown by means of the asymptotic limits, where K and L are sent to zero or to inﬁnity, respectively. Eventually, we establish higher regularity for the phase-ﬁelds, and we prove the uniqueness of weak solutions given that the mobility functions are constant.


Introduction
In recent times, the description of immiscible two-phase flows in a bounded domain by diffuse-interface models has become very popular.This is especially because those systems can usually be handled more easily in terms of mathematical analysis than their sharpinterface counterparts.In such models, the location of the two fluids inside the container is represented by an order parameter, the so-called phase-field.The time evolution of this phase-field function is often described by a convective Cahn-Hilliard equation.There, the velocity field of the mixture enters via the material derivative of the phase-field.The time evolution of the velocity field is usually described by a fluid equation, e.g., the Navier-Stokes equation.Very frequently used Navier-Stokes-Cahn-Hilliard models for two-phase flows are the Model H (see [25,27]), which covers the case where both fluids have the same individual density, and the Abels-Garcke-Grün model (AGG model) (see [1]), which is even capable of describing the situation where both fluids have different (i.e., unmatched ) densities.
For the Model H and the AGG model, the standard choice are homogeneous Neumann boundary conditions for the Cahn-Hilliard quantities (i.e., the phase-field and the chemical potential) as well as a no-slip boundary condition for the velocity field.However, these choices lead to some crucial limitations of these models: • The contact angle between the diffuse interface and the boundary of the domain is fixed at ninety degrees at all times, which is unrealistic for many real-world applications.• The motion of the contact line, where the diffuse interface intersects the boundary of the domain, is driven only by diffusion.This means that a motion of the contact line caused by convection is not taken into account.
• No transfer of material between the bulk and the boundary of the domain is allowed.
Therefore, any absorption of material by the boundary cannot be described.
For a more detailed discussion of these issues, we refer the reader to [23].
To overcome these limitations, a new Navier-Stokes-Cahn-Hilliard model with dynamic boundary conditions was derived in [23].It can be regarded as an extension of the AGG model and is capable of describing two-phase flows with unmatched densities.In the case of matched densities, some first analytic results (namely the existence of weak solution and their uniqueness under certain additional assumptions) were presented in [23].The case of unmatched densities, however, is much more involved.This is mainly because in the Navier-Stokes equation, the density function then depends on the phase-field, and a further flux term related to the interfacial motion will appear.In the case of unmatched densities, at least for certain parameter choices in the dynamic boundary conditions, the existence of a weak solution to the model introduced in [23] was shown in [17] by an implicit time discretization scheme.A different strategy to prove the existence of weak solutions to the model proposed in [23], in a unified framework for all parameter choices in the dynamic boundary conditions, is to first analyze the underlying bulk-surface convective Cahn-Hilliard subsystem separately.Using this information, a weak solution to the full Navier-Stokes-Cahn-Hilliard system is then to be constructed by means of a suitable fixed point argument.
The function φ : Q → R is the bulk phase-field.It is an order parameter, which represents the distribution of the two immiscible materials within the domain Ω.Similarly, the surface phase-field ψ : Σ → R represents the distribution of two materials on the surface.The functions F and G are double-well potentials, which lead to the effect that φ(t) and ψ(t) will attain values close to ±1 in most parts of Ω and Γ (which correspond to the pure phases of the materials).In the remaining intermediate regions, which are called diffuse interfaces, φ and ψ are expected to exhibit a continuous transition between the values −1 and 1.The functions µ : Q → R and θ : Σ → R represent the bulk chemical potential and the surface chemical potential, respectively.The non-negative, scalar functions m Ω (φ) and m Γ (ψ) are called mobilities.They depend on the phase-fields and describe where diffusion processes occur.Moreover, the functions v : Q → R d and w : Σ → R d are prescribed velocity fields that correspond to the flow of the two materials in the bulk and on the surface, respectively.In this paper, we always assume that v•n = 0 and w•n = 0 on Σ.We point out that in many cases (e.g., in the model derived in [23]), v and w are related by the condition w = v| Σ on Σ.However, as this relation will not have any impact on our mathematical analysis, we consider v and w as independent functions to keep our model as general as possible.Furthermore, ǫ > 0 is a parameter related to the thickness of the diffuse interface in the bulk, whereas ǫ Γ > 0 corresponds to the width of the diffuse interface on the boundary.The parameter κ > 0 acts as a weight for the surface Dirichlet energy ψ → Γ |∇ Γ ψ| 2 dΓ, which has a smoothing effect on the phase separation at the boundary.
The time evolution of (φ, µ) is determined by the bulk convective Cahn-Hilliard subsystem (1.1a), (1.1b) , whereas the evolution of (ψ, θ) is described by the surface convective Cahn-Hilliard subsystem (1.1c), (1.1d) , which is coupled to the bulk by expressions involving the normal derivatives ∂ n φ and ∂ n µ.The bulk and surface quantities are further coupled by the boundary conditions (1.1e) and (1.1f), which depend on parameters K, L ∈ [0, ∞] and α, β ∈ R.
In (1.1e) and (1.1f), the parameters K, L ∈ [0, ∞] are used to distinguish different cases, each corresponding to a certain solution behaviour related to a physical phenomenon.The case that has been studied most extensively in the literature is K = 0.In this case, (1.1e) is to be interpreted as the Dirichlet type boundary condition φ = αψ on Σ.This choice makes particular sense along with α = 1, if the materials on the boundary are simply considered as an extension of those in the bulk.Originally, the boundary condition (1.1f) was introduced in the non-convective case (i.e., v ≡ 0 and w ≡ 0) in the following literature: • The choice L = 0 was proposed in [14] and [24].Then, (1.1f) can be restated as the Dirichlet condition µ = βθ on Σ, which means that the chemical potentials µ and θ are always in a chemical equilibrium.In this case, a rapid transfer of material between bulk and boundary can be expected (see, e.g., [29]).• The choice L = ∞ was introduced in [32].In this case, (1.1f) is a homogeneous Neumann boundary condition on µ, which means that the mass flux between bulk and surface is zero.Consequently, no transfer of material between bulk and surface will occur.• The choice L ∈ (0, ∞) was first used in [29] to interpolate between the extreme cases L = 0 and L = ∞.In this case, a transfer of material between bulk and surface will occur, and the number L −1 is related to the rate of absorption of bulk material by the boundary (cf.[29]).
In the convective case (i.e., with non-trivial velocity fields v and w), the boundary condition (1.1f) was also used in the Navier-Stokes-Cahn-Hilliard model derived in [23] and in the Cahn-Hilliard-Brinkman model studied in [12].Furthermore, the model derivation in [23] shows that the parameter β acts as a weight for the mass flux between bulk and surface.This means that β can also be negative and therefore, we simply assume β ∈ R.
However, we also want to cover the case, where the materials on the boundary are not the same as those in the bulk.For instance, this is the case if the bulk materials are transformed at the boundary by a chemical reaction.In this case, the surface phase-field might not be proportional to the trace of the bulk phase-field, so the choice K = 0 would not be appropriate.In our model, this is taken into account by the choice K ∈ (0, ∞].In the case K = ∞, (1.1f) degenerates to a homogeneous Neumann boundary condition for the phase-field.For most applications, this might not be the preferred choice as it is known that such a Neumann boundary condition enforces the diffuse interface to always intersect the boundary at a perfect angle of ninety degrees.In fact, this is one of the aforementioned limitations we usually want to overcome by the usage of dynamic boundary conditions.However, we include the case K = ∞ for the sake of completeness.In the case K ∈ (0, ∞), (1.1f) can be regarded as a Robin type boundary condition.It is suitable to describe a scenario, where ψ and the trace of φ are not proportional, but it still allows for dynamical changes of the contact angle between the diffuse interface and the boundary.In the context of bulk-surface Cahn-Hilliard equations, condition (1.1f) with K ∈ (0, ∞) was used in [28] (in the non-convective case) and in [12] (in the convective case).In [28] and [12], it was further shown that the limit K → 0 can be used to recover the boundary condition (1.1f) with K = 0. Especially in the case that the materials on the boundary differ from those in the bulk, the parameter α could be any real number (even with a negative sign).Therefore, to keep the model as general as possible, we allow for α ∈ R.However, for our mathematical analysis, we will need the additional relation αβ |Ω| + |Γ| = 0. Of course, this is trivially satisfied if α and β have the same sign.
We further want to point out that including the case K ∈ (0, ∞) also helps our mathematical analysis.This is because the existence of a weak solution to (1.1) can be shown by a suitable Faedo-Galerkin scheme only in the cases, where the boundary conditions are of the same type (i.e., K = L = 0, K = L = ∞ or 0 < K, L < ∞), as otherwise, the spaces of admissible test functions in the weak formulation do not match.In this paper, we will first construct a weak solution of (1.1) in the case 0 < K, L < ∞.Then, we prove the existence of a weak solution in the remaining cases by sending the parameters K and L to zero or to infinity, respectively.
Let us now discuss some important properties of our model.The system (1.1) is associated with the total energy where the function is used to distinguish the different cases corresponding to the choice of K. Sufficiently regular solutions of the system (1.1) satisfy the mass conversation law for all t ∈ [0, T ].As mentioned above, this means that a transfer of material between bulk and surface is allowed only in the cases L ∈ [0, ∞).Moreover, sufficiently regular solutions satisfy the energy identity for all t ∈ [0, T ].We point out that in the non-convective case (i.e., v = 0 and w = 0), the right-hand side of (1.4) is clearly non-positive.This means that the energy dissipates over the course of time and the terms on the right-hand side of (1.4) can be interpreted as the dissipation rate.In this case, the system (1.1) can be derived as an H −1 type gradient flow of the free energy E K subject to the inner product (•, •) L,β, * that will be introduced in Section 2, (P5) (cf.[29, Remark 2.2] and [28,Section 3]).Moreover, if L = ∞, the system (1.1) can also be derived by an energetic variational approach that combines the least action principle and Onsager's principle of maximum energy dissipation (cf.[32, Section 2] and [28,Appendix]).
In the convective case (i.e., v and w are non-trivial), the energy identity (1.4) does, in general, not imply dissipation of the energy.However, if the velocity field is not just prescribed but determined by a Navier-Stokes equation (cf.[23]) or a Brinkman/Stokes equation (cf.[12]), an energy dissipation law for the corresponding total energy can be obtained.
The convective Cahn-Hilliard equation with dynamic boundary conditions was analyzed, for instance, in [10,11,22], and also in [15][16][17]23] as part of a system in which the velocity field is described by an additional fluid equation.
For more information on the Cahn-Hilliard equation with classical homogeneous Neumann boundary conditions or dynamic boundary conditions, we refer to the recent review paper [41] as well as the book [36].
Structure of this paper.In Section 2 we collect some notation, assumptions, preliminaries and important tools.After introducing the notion of weak solutions of (1.1), our main results are stated in Section 3. In Section 4, we construct a weak solution in the case (K, L) ∈ (0, ∞) 2 via a Faedo-Galerkin approach.Afterwards, in Section 5, we investigate the asymptotic limits K → 0, K → ∞, L → 0 and L → ∞, which prove the existence of weak solutions of (1.1) in the limit cases.In Section 6, under suitable additional assumptions, we establish higher spatial regularity for the phase-fields in the context of weak solutions.In the case of constant mobilities, we eventually prove the continuous dependence of weak solutions on the velocity fields and the initial data.In particular, this continuous dependence result directly entails the uniqueness of the weak solution.
2 Notation, assumptions, and preliminaries In this section, we introduce some notation, assumptions, and preliminaries that are supposed to hold throughout the remainder of this paper.

Notation
Let us first introduce some basic notation.
(N2) Let Ω ⊂ R d with d ∈ {2, 3} be a bounded Lipschitz domain in R d , and let Γ := ∂Ω denote its boundary.For any s ≥ 0 and p ∈ [1, ∞], the Lebesgue and Sobolev spaces for functions mapping from Ω to R are denoted as L p (Ω) and W s,p (Ω).We write • L p (Ω) and • W s,p (Ω) to denote the standard norms on these spaces.In the case p = 2, we use the notation H s (Ω) = W s,2 (Ω).In particular, H 0 (Ω) can be identified with L 2 (Ω).The Lebesgue and Sobolev spaces on Γ are denoted by L p (Γ) and W s,p (Γ) along with the corresponding norms • L p (Γ) and • W s,p (Γ) , respectively.For vectorvalued functions mapping from Ω into R d , we use the notation L p (Ω), W s,p (Ω) and H s (Ω).The spaces L p (Γ), W s,p (Γ) and H s (Γ) are defined analogously.For any real numbers s ≥ 0 and p ∈ [1, ∞] and any Banach space X, the Bochner spaces of functions mapping from an interval I into X are denoted by L p (I; X) and W s,p (I; X).Furthermore, for any interval I and any Banach space X, the space C(I; X) denotes the set of continuous functions mapping from I to X.
(N3) For any Banach space X, its dual space is denoted by X ′ .The corresponding duality pairing of elements φ ∈ X ′ and ζ ∈ X is denoted by φ, ζ X .If X is a Hilbert space, we write (•, •) X to denote its inner product.
to denote the generalized means of u and v, respectively.Here, |Ω| denotes the ddimensional Lebesgue measure of Ω, whereas |Γ| denotes the (d − 1)-dimensional Hausdorff measure of Γ.If u ∈ L 1 (Ω) or v ∈ L 1 (Γ), the generalized mean can be expressed as respectively.
(N5) For any bounded domain Ω ⊂ R d (d ∈ N) with Lipschitz boundary Γ := ∂Ω, we introduce the spaces We point out that in the definition of L 3 div (Ω), the relation div v = 0 in Ω has to be understood in the sense of distributions.This already implies that v • n ∈ H −1/2 (Γ), and therefore, the relation v • n = 0 on Γ is well-defined.

Assumptions
We make the following general assumptions.
(A3) The mobility functions m Ω : R → R and m Γ : R → R are bounded, continuous and uniformly positive.This means that there exist constants (A4) The potentials F : R → [0, ∞) and G : R → [0, ∞) are continuously differentiable and there exist exponents p and q satisfying as well as constants c F ′ , c G ′ ≥ 0 such that the first-order derivatives satisfy the growth conditions ) for all s ∈ R.These assumptions already imply the existence of constants c F , c G ≥ 0 such that F and G satisfy the growth conditions for all s ∈ R.
(A5) There exist constants a F , a G > 0 and b F , b G ≥ 0 such that the potentials F and G introduced in (A4) satisfy for all s ∈ R.
(A6) The potentials F and G introduced in (A4) are twice continuously differentiable, and there exist constants c F ′′ , c G ′′ ≥ 0 such that the second-order derivatives of F and G satisfy the growth conditions ) for all s ∈ R. We point out that these assumptions already imply the growth conditions (2.1)-(2.4)stated in (A4).
Remark 2.1.A standard choice for F and G is the polynomial double-well potential It satisfies the assumptions (A4)-(A6) with p = q = 4.However, singular potentials such as the logarithmic Flory-Huggins potential or the double-obstacle potential are not admissible as they do not satisfy any polynomial growth condition.Nevertheless, the investigation of system (1.1) with singular potentials is an interesting topic for future research.

Preliminaries
(P1) For any real numbers s ≥ 0 and p ∈ [1, ∞], we set L p := L p (Ω) × L p (Γ), and H s := H s (Ω) × H s (Γ), provided that the boundary is sufficiently regular.As usual, we identify L 2 with H 0 .Note that H s is a Hilbert space with respect to the inner product H s .We recall that the duality pairing can be expressed as (P2) Let L ∈ [0, ∞] and β ∈ R. We introduce the closed linear subspace and define Endowed with the inner product (•, •) H 1 and its induced norm, the space . By means of the Riesz representation theorem, this product can be extended to a duality pairing on (H 1 L,β ) ′ × H 1 L,β , which will also be denoted as Note that V 1 L,β is a Hilbert space with respect to the inner product (•, •) H 1 and its induced norm.
and we define a bilinear form on for all (φ, ψ) ∈ H 1 .The bilinear form (•, •) L,β defines an inner product on V 1 L,β , and (P5) For any L ∈ [0, ∞] and β ∈ R, we define the space Using the Lax-Milgram theorem, one can show that for any (φ, ψ) ∈ V −1 L,β , there exists a unique weak solution to the following elliptic problem with bulk-surface coupling: ) This means that S L,β (φ, ψ) satisfies the weak formulation for all (φ, ψ) ∈ V −1 L,β , for a constant C ≥ 0 depending only on Ω, L and β.We can thus define the solution operator as well as an inner product and its induced norm on For the case L ∈ (0, ∞), we refer the reader to [30, Theorem 3.3 and Corollary 3.5] for a proof of these statements.In the other cases, the results can be proven analogously.
(P6) We further recall the following bulk-surface Poincaré inequalitiy, which has been established in [30, Lemma A.1]: Let K ∈ [0, ∞) and α, β ∈ R with αβ |Ω| + |Γ| = 0 be arbitrary.Then there exists a constant C P > 0 depending only on K, α, β and Ω such that We conclude this section by presenting a simple inequality that will be frequently used in our mathematical analysis.
Hence, the claim directly follows.

Main results
As mentioned in (A2), we set ǫ = ǫ Γ = κ = 1.This does not mean any loss of generality as the exact values of ǫ, ǫ Γ and κ do not have any impact on the mathematical analysis (as long as they are positive).By this choice, the system (1.1) can be restated as follows: ) ) The total energy associated with this system reads as We now introduce the notion of a weak solution to system (3.1).
(iii) The functions φ, ψ, µ and θ satisfy the weak formulation The functions φ and ψ satisfy the mass conservation law The functions φ, ψ, µ and θ satisfy the energy inequality Now, we are ready to formulate the main results of this paper.The first main result provides the existence of a weak solution to system (3.1) in all cases (K, Theorem 3.2.(Existence of weak solutions of (3.1)) Suppose that the assumptions ) be given velocity fields.In the case K ∈ {0, ∞}, we further assume that (A5) holds.Then, there exists a weak solution (φ, ψ, µ, θ) of the system (3.1) in the sense of Definition 3.1, which has the additional regularity Let us now assume that (A6) holds, that the domain Ω is of class C k for k ∈ {2, 3}, and in the case d = 3, we further assume that (A4) holds with p ≤ 4.Then, we additionally have and the equations are fulfilled in the strong sense.
Remark 3.3.The continuity property (3.10) also holds in the case K = L = 0 provided that α = 0, β = 0 and the potentials F and G satisfy the compatibility condition In that case, we have (F ′ (φ), G ′ (ψ)) ∈ H 1 0,β = D β , which allows us to use Proposition A.1 similarly as in Section 6.1.
Proof of Theorem 3.2.The existence of a weak solution in the sense of Definition 3.1 is established in Moreover, in the case K ∈ (0, ∞], the additional regularity (3.7) follows from the corresponding aforementioned theorems.All remaining results are shown in Theorem 6.1.
Remark 3.4.We point out that Theorem 5.2, Theorem 5.4, Theorem 5.6 and Theorem 5.8 are not only useful to prove Theorem 3.2, but also provide further valuable insights about the asymptotic limits K → 0, K → ∞, L → 0 and L → ∞, respectively, on the level of weak solutions.
Our second main result shows continuous dependence of weak solutions on the velocity fields and the initial data in the case of constant mobilities.As a direct consequence, this entails uniqueness of the weak solution.
The proof of Theorem 3.5 will be presented in Subsection 6.2.
Proof.Step 1: Discretization via a Faedo-Galerkin scheme.It is well known that the problems have countable many eigenvalues, which can be written as increasing sequences {λ j Ω } j∈N and {λ j Γ } j∈N , respectively.The associated eigenfunctions {ζ j } j∈N ⊂ H 1 (Ω) and {ξ j } j∈N ⊂ H 1 (Γ) form an orthonormal Schauder basis of L 2 (Ω) and L 2 (Γ), respectively.In particular, we fix the eigenfunctions associated to the first eigenvalues . Moreover, the eigenfunctions {ζ j } j∈N and {ξ j } j∈N also form an orthogonal Schauder basis of H 1 (Ω) and H 1 (Γ), respectively.
For any m ∈ N, we introduce the finite-dimensional subspaces along with the orthogonal L 2 (Ω)-projection P Am , and the orthogonal L 2 (Γ)-projection P Bm .
In particular, there exist constants C Ω , C Γ > 0 depending only on Ω and Γ, respectively, such that for all ζ ∈ H 1 (Ω) and ξ ∈ H 1 (Γ), For any m ∈ N, and t ∈ [0, T ], we make the ansatz Here, the scalar, time-dependent coefficients a m j , b m j , c m j , d m j , j = 1, . . ., m are assumed to be continuously differentiable functions that are not yet determined.They need to be designed in a way such that the discretized weak formulation m > 0 to the corresponding initial value problem.Without loss of generality, we assume that T * m ≤ T and that (a m , b m ) ⊤ is non-extendable, i.e., T * m is chosen as large as possible.Now, we can reconstruct (c m , d m ) ⊤ : [0, T * m ) → R 2m by the aforementioned 2m-dimensional system of algebraic equations.In view of (4.1), we thus have shown the existence of functions solving (4.2) on the time interval [0, T * m ) subject to the initial conditions (4.3).
Step 2: Uniform estimates.We establish suitable estimates for each approximate solution (φ m , ψ m , µ m , θ m ), which are uniform with respect to the index m.In particular, let T m < T * m be arbitrary.In the following, let C denote a generic non-negative constant depending only on the initial data and the constants introduced in (A3)-(A4) including the final time T , but not on m or T m .Testing (4.2a) with (µ m , θ m ), (4.2b) with − (∂ t φ m , ∂ t ψ m ), adding the resulting equations, and integrating with respect to time from 0 to t, we derive the discrete energy identity for all t ∈ [0, T m ].Now, due to the growth conditions on F, G (see (A4)), the Sobolev embeddings H 1 (Ω) ֒→ L 6 (Ω) and H 1 (Γ) ֒→ L q (Γ), the trace embedding and the properties of the projections P Am and P Bm , we use the initial conditions to infer To obtain a suitable bound on the energy and the integral terms on the left-hand side of (4.4), we have to deal with the convective terms appearing on the right-hand side.To this end, we start by deriving for all t ∈ [0, T m ] the estimate using again the embeddings H 1 (Ω) ֒→ L 6 (Ω) and H 1 (Γ) ֒→ L 6 (Γ).Recalling that F, G ≥ 0, and that the mobility functions m Ω and m Γ are uniformly positive (see (A3)), we infer from (4.5) and (4.6) that Here, the last inequality follows from the bulk-surface Poincaré inequality (P6), which usage is justified by the following reasoning.First, testing (4.2a) with (β, 1) yields the discrete mass conservation law for all t ∈ [0, T m ], where we used that Hence, considering and invoking (4.8), we infer This allows us to use the bulk-surface Poincaré inequality (P6), which yields Thus, using Gronwall's lemma, we conclude from (4.7) that ds for all t ∈ [0, T m ].We thus readily deduce that In view of (4.7), this additionally entails Now, due to Lemma 2.2, we have on [0, T m ].Hence, after taking the supremum over all (ζ * , ξ * ) ∈ H 1 satisfying (ζ * , ξ * ) H 1 ≤ 1, taking the square of the resulting inequality and integrating in time over [0, T m ], we deduce that Exploiting the uniform estimates (4.10) and (4.16) we conclude Step 3: Extension of the approximate solution onto [0, T ].Using the definition of the approximate solution (4.1) and the uniform estimate (4.10), we obtain for any Since the coefficients (c m , d m ) ⊤ can be reconstructed from (a m , b m ) ⊤ , we further conclude that the coefficients (c m , d m ) ⊤ also exist on [0, T ].This automatically entails that the approximate solution (φ m , ψ m , µ m , θ m ) exists on [0, T ] and satisfies the discretized weak formulation (4.2) on [0, T ] for all m ∈ N. Additionally, as the particular choice of T m did not play any role in the derivation of the uniform estimates established in Step 2, we infer that the estimates (4.10), (4.11), (4.15), (4.16) and (4.17) remain true after replacing T m with T .In summary, we conclude that for each m ∈ N, the approximate solution (φ m , ψ m , µ m , θ m ) satisfies the uniform estimate
Since the approximate solutions (φ m , ψ m ) satisfy the discrete mass conservation law (see (4.8)) for all m ∈ N, we can use the convergence (4.20) to deduce that (φ, ψ) satisfies the mass conservation law (3.5) as well.Thus, Definition 3.1(iv) is fulfilled.
To verify the energy inequality (3.6), we consider an arbitrary non-negative test function σ ∈ C ∞ c (0, T ).Multiplying (4.4) with σ and integrating in time from 0 to T yields As m → ∞, we infer from (4.28) that Additionally, using (4.24) and (A4), Fatou's lemma implies Moreover, from the strong convergence in (4.26) we obtain as m → ∞.Lastly, in view of (4.20) and the weak lower semicontinuity of norms, we have As a consequence of (4.32), (4.33) and the weak lower semicontinuity of norms with respect to the weak convergence, another application of Fatou's lemma entails that Further, due to the convergences (4.36) and (4.37), we find that Lastly, recalling the growth conditions on F and G (see (A4)), we use the initial conditions (4.3) as well as the convergence properties of the projections P Am and P Bm together with Lebesgue's general convergence theorem to infer that for almost every t ∈ [0, T ].To prove that (4.49) holds for every t ∈ [0, T ], first note that all integral terms in (4.49) depend continuously on time.Furthermore, due to (φ, ψ) ∈ C([0, T ]; L 2 ), we deduce that the following functions are lower semicontinuous on [0, T ].For the first two functions, this is a consequence of the lower semicontinuity of the respective norms, while for the last two functions, it follows from Fatou's lemma.This already entails that the weak solution (φ, ψ, µ, θ) satisfies the energy inequality for all times t ∈ [0, T ].Thus, Definition 3.1(v) is fulfilled.
We have therefore shown that the quadruplet (φ, ψ, µ, θ) is a weak solution of system (3.1) in the sense of Definition 3.1.
Remark 4.2.This proof should work similarly in the cases K = L = 0 and K = L = ∞.
In the first case, one uses a Faedo-Galerkin scheme as in [23] based on eigenfunctions of a suitable bulk-surface elliptic problem with Dirichlet type coupling condition (cf.[30,Theorem 3.3]).In the second case, the bulk-surface Cahn-Hilliard equation (3.1) reduces to two, uncoupled Cahn-Hilliard equations, one in Ω and one on Γ.Therefore, a Faedo-Galerkin basis as in the above proof can be used.Also note that, while we do not have the bulk-surface Poincaré inequality at our disposal anymore, we can use the standard Poincaré inequality for functions with vanishing mean, since the approximate solutions satisfy the mass conservation law (3.5)for L = ∞.
5 Asymptotic limits and existence of weak solutions to the limit models In this section, we investigate the asymptotic limits K → 0 and K → ∞, and L → 0 and L → ∞ of the system (3.1).
Remark 5.1.In this section, we will need the additional assumption (A5) when approaching the limit cases K ∈ {0, ∞}.This is because the constant in the bulk-surface Poincaré inequality (P6) depends on K in some way, but we do not know this dependence explicitly.In particular, it is unclear how this constant behaves if we send K to zero or infinity, respectively.Therefore, we cannot rely on this Poincaré inequality to obtain suitable uniform bounds, but instead, we use (A5) to directly obtain uniform bounds from the energy functional.
The mass conservation law follows by passing to the limit m → ∞ in the mass conservation law for (φ Km , ψ Km ).Alternatively, one can test the weak formulation (3.4a) with (β, 1) ∈ H 1 , integrate with respect to time from 0 to t, and employ the fundamental theorem of calculus to infer (3.5).We conclude that Definition 3.1(iv) is satisfied.
As the boundary term involving K m in the energy is non-negative, we note that for all t ∈ [0, T ], and thus for all non-negative test functions σ ∈ C ∞ c (0, T ).Here, the first inequality follows by proceeding similarly as in the proof of Theorem 3.2.We now use the energy inequality (3.6) written for (φ Km , ψ Km , µ Km , θ Km ) to further bound the right-hand side of (5.15).Using lower semicontinuity arguments (similar to those in the proof of Theorem 3.2), and recalling (5.2), passing to the limit m → ∞ leads to the corresponding energy inequality for (φ * , ψ * , µ * , θ * ).This proves that the quadruplet (φ * , ψ * , µ * , θ * ) is a weak solution of (3.1) in the sense of Definition 3.1 with K = 0.
Remark 5.5.Suppose that for any K ∈ (0, ∞), the phase-field φ K is sufficiently regular such that the boundary condition (3.1e) holds in the strong sense, that is This is actually fulfilled under additional assumptions on the regularity of Γ and the parameter p (if d = 3), see Theorem 6.1.Then, estimate (5.16) can be reformulated as As the right-hand side of this inequality tends to zero as K → ∞, this explains why the homogeneous Neumann boundary condition ∂ n φ * = 0 a.e. on Σ appears in the limit model corresponding to K = ∞.
Proof of Theorem 5.4.We consider an arbitrary sequence (K m ) m∈N ⊂ (0, ∞) such that K m → ∞ as m → ∞.Without loss of generality, we assume that K m ∈ [1, ∞) for all m ∈ N.For any m ∈ N, let (φ Km , ψ Km , µ Km , θ Km ) denote a weak solution of the system (3.1) in the sense of Definition 3.1 with initial data (φ 0 , ψ 0 ) corresponding to the parameter K m .In this proof, we use the letter C to denote generic positive constants independent of K m and m.Let now m ∈ N be arbitrary.
As we have seen in the proof of Theorem 3.2 (see 4.5), the initial energy satisfies since K m ≥ 1 for all m ∈ N.This allows us to use the same argumentation as in the proof of Theorem 5.2, and we infer This already implies (5.16).We now test the weak formulation (3.4b) with (ζ, ξ) ∈ H 1 .Using again that K m ≥ 1, we find that This estimate allows us to apply Lemma 2.2, and using (5.18) and (5.19), we conclude In combination with (5.18) we thus deduce Lastly, exploiting (5.18), we conclude from the weak formulation (3.4a) by a comparison argument that In summary, combining (5.18) and (5.20)-(5.22),we thus have shown that In view of the uniform estimate (5.23), the Banach-Alaoglu theorem and the Aubin-Lions-Simon lemma imply the existence of functions (φ * , ψ * , µ * , θ * ) such that φ Km , ψ Km → (φ * , ψ * ) weakly-star in L ∞ (0, T ; H 1 ), strongly in C([0, T ]; H s ) for all s ∈ [0, 1), (5.25) as m → ∞, along a non-relabeled subsequence.Additionally, thanks to (5.18) and the trace theorem, we infer as m → ∞.We readily deduce that Definition 3.1(i) is satisfied.
Using the Sobolev embedding H s ֒→ L 2 for s ∈ (0, 1) along with (5.25), we find that as m → ∞.In view of the initial conditions satisfied by (φ Km , ψ Km ), we deduce that Definition 3.1(ii) is fulfilled.
As in the proof of Theorem 5.2, the convergences (5.24)-(5.27)are sufficient to pass to the limit in the weak formulations (3.4) associated with K m to conclude that the limit (φ * , ψ * , µ * , θ * ) fulfils the weak formulation (3.4) associated with K = ∞.Thus, we infer that Definition 3.1(iii) is satisfied.
The mass conservation law as stated in Definition 3.1(iv) also follows with the same reasoning as in the proof of Theorem 5.2.
For the energy inequality, note that for all non-negative σ ∈ C ∞ c (0, T ), we have ( The first inequality follows the same argumentation that we already have seen before, whereas the second inequality is due to (5.27).We finish the proof by mimicking the proof of the energy inequality in Theorem 3.2, which is based on the convergence results (5.24)-(5.26)and the assumptions (A3)-(A4).
We thus have shown that the quadruplet (φ * , ψ * , µ * , θ * ) is a weak solution to the system (3.1) in the sense of Definition 3.1 for K = ∞.

5.3
The limit L → 0 and the existence of a weak solution if In this section, we investigate the asymptotic limits L → 0 and L → ∞ of the system (3.1) for fixed K ∈ [0, ∞].
Proof of Theorem 5.6.We consider an arbitrary sequence (L m ) m∈N ⊂ (0, ∞) such that L m → ∞ as m → ∞ and a corresponding weak solution (φ Lm , ψ Lm , θ Lm , µ Lm ) to the system (3.1) in the sense of Definition 3.1 to the initial data (φ 0 , ψ 0 ).In this proof, we denote by C arbitrary positive constants independent of L m and m, which may change their value from line to line.Let now m ∈ N be arbitrary.
If K ∈ [0, ∞), we conclude from the energy inequality (3.6) that If K = ∞, we do not have the bulk-surface Poincaré inequality at our disposal.Instead, we have to argue as in the proof of Theorem 5.2 and make use of the additional assumption (A5) to obtain the estimate (5.31).In particular, (5.31) yields (5.30).Arguing as in the proof of Theorem 3.2, we additionally infer that If K ∈ (0, ∞], we may use Lemma 2.2 and further obtain (µ Lm , θ Lm ) L 4 (0,T ;L 2 ) ≤ C, (5.33) which in combination with (5.31) yields (µ Lm , θ Lm ) L 2 (0,T ;H 1 ) ≤ C. (5.34) In the case K = 0, we argue analogously as in the proof of Theorem 5.2 to deduce that from which we obtain (5.34) as well.
For the time derivatives, we again proceed similarly as in the proof of Theorem 3.2 but choose the test function space D β instead of H 1 .We then obtain due to the weak formulation (3.4a) and the uniform bound (5.31) that βθ Lm − µ Lm → βθ * − µ * weakly in L 2 (0, T ; L 2 (Γ)), (5.40) as m → ∞, along a non-relabeled subsequence.Furthermore, we conclude from (5.31) that as m → ∞.In combination with (5.40) we infer that µ * = βθ * a.e. on Σ due to the uniqueness of the limit.Proceeding similarly as in the case K → 0 (see the proof of Theorem 5.2), we eventually show that the quadruplet (φ * , ψ * , µ * , θ * ) is a weak solution of the system (3.1) in the sense of Definition 3.1 for L = 0.
Proof of Theorem 5.8.We consider an arbitrary sequence (L m ) ⊂ (0, ∞) such that L m → ∞ as m → ∞ and a corresponding weak solution (φ Lm , ψ Lm , θ Lm , µ Lm ) to (3.1) in the sense of Definition 3.1 to the initial data (φ 0 , ψ 0 ).Since L m → ∞ as m → ∞, we can assume without loss of generality that L m ∈ [1, ∞) for all m ∈ N. In this proof, we use the letter C to denote generic positive constants independent of L m and m, which may change their value from line to line.Let now m ∈ N be arbitrary.
Using the same argumentation as in the proof of Theorem 5.6, we infer For the time derivatives, we proceed similarly as in the proof of Theorem 3.2 and obtain due to the weak formulation (3.4a) and the uniform bounds (5.43) that (5.47) Here we additionally used L m ≥ 1 for all m ∈ N. In view of the uniform estimates (5.43), (5.44)-(5.46)and (5.47), the Banach-Alaoglu theorem and the Aubin-Lions-Simon lemma imply the existence of functions φ * , ψ * , µ * and θ * such that φ Lm , ψ Lm → (φ * , ψ * ) weakly-star in L ∞ (0, T ; H 1 K,α ), strongly in C([0, T ]; H s ) for all s ∈ [0, 1), (5.49) as m → ∞, along a non-relabeled subsequence.Due to (5.43), we additionally have as m → ∞.Arguing further as in the case K → ∞ (see the proof of Theorem 5.4), we eventually show that the quadruplet (φ * , ψ * , µ * , θ * ) satisfies Definition 3.1(i)-3.1(iii)and 3.1(v).For the mass conservation law, simply note that the weak formulations for φ * and ψ * , as well as the test functions, are not coupled, which allows us to use ζ ≡ 1 and ξ ≡ 1 as test functions, respectively.Integrating the resulting equations in time from 0 to t, and employing the fundamental theorem of calculus, we infer We conclude this section with the following remark.
Remark 5.10.If the uniqueness of the respective limiting models is known (e.g. if the mobility functions are constant, cf.Theorem 3.5), we even obtain convergence of the full sequence (not only a subsequence) in Theorem 5.2, Theorem 5.4, Theorem 5.6 and Theorem 5.8, respectively, by means of a standard subsequence argument.
6 Higher regularity, continuous dependence and uniqueness In this section, we present the results on higher regularity for the phase-fields as well as the continuous dependence and uniqueness of weak solutions to the system (3.1) as stated in Theorem 3.2 and Theorem 3.5, respectively.
Remark 6.2.In the other cases, where the boundary conditions involving K and L are of the same type (i.e., K = L = 0 or K = L = ∞), it should also be possible to obtain the regularities (6.1) and (6.2) even for d = 3 and p ∈ (4,6).This is because in these cases, the system (1.1) can be discretized by a Faedo-Galerkin scheme (cf.Remark 4.2) and therefore, the higher order regularity estimates can be performed on the level of the approximate solutions.In the case K = L = 0, we refer the reader to the proof of [23,Theorem 3.3], where this line of argument is carried out in detail.

Continuous dependence and uniqueness
Proof of Theorem 3.5.As the functions m Ω and m Γ are constant, we assume, without loss of generality, that m Ω ≡ 1 and m Γ ≡ 1.In the following, we use the letter C to denote generic positive constants depending only on Ω, the parameters of the system (3.1), and the initial data and the prescribed velocity field.

(4. 12 )
on [0, T m ].Taking the supremum over all (ζ * , ξ * ) ∈ H 1 with (ζ * , ξ * ) H 1 ≤ 1, and exploiting (4.10), we conclude that This shows that the solution (a m , b m ) ⊤ is bounded on the time interval [0, T * m ) by a constant that is independent of m and T * m .Hence, by classical ODE theory, this allows us to extend the solution beyond the time T * m .However, as (a m , b m ) ⊤ was assumed to be non-extendable, this is a contradiction.We thus infer that the solution (a m , b m ) ⊤ actually exists on [0, T ].

. 1 )
and the limit (φ * , ψ * , µ * , θ * ) is a weak solution to the system (3.1) in the sense of Definition 3.1 with K = 0. Remark 5.3.(a) As the right-hand side in (5.1) tends to zero as K → 0, this explains why the Dirichlet type boundary condition φ * = αψ * a.e. on Σ appears in the limit model corresponding to K = 0.