{'fraction_non_alphanumeric': 0.07664195738945623, 'fraction_numerical': 0.026385361049193926, 'mean_word_length': 3.5888439581351093, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On a Cahn-Hilliard system with source term and thermal memory

A nonisothermal phase field system of Cahn-Hilliard type is introduced and analyzed mathematically. The system constitutes an extension of the classical Caginalp model for nonisothermal phase transitions with a conserved order parameter. It couples a Cahn-Hilliard type equation with source term for the order parameter with the universal balance law of internal energy. In place of the standard Fourier form, the constitutive law of the heat flux is assumed in the form given by the theory developed by Green and Naghdi, which accounts for a possible thermal memory of the evolution. This has the consequence that the balance law of internal energy becomes a second-order in time equation for the thermal displacement or freezing index, that is, a primitive with respect to time of the temperature. Another particular feature of our system is the presence of the source term in the equation for the order parameter, which entails additional mathematical difficulties because the mass conservation of the order parameter is lost. We provide several mathematical results under general assumptions on the source term and the double-well nonlinearity governing the evolution: existence and continuous dependence results are shown for weak and strong solutions to the corresponding initial-boundary value problem.


Introduction
A common assumption in phase segregation processes of binary mixtures is to postulate that the mixture under investigation undergoes the phase separation at a constant temperature. However, in numerous applications the evolution does not take place under isothermal conditions. The first contribution aiming at including temperature effects in the theory of phase separation is due to Caginalp [7][8][9]. It was motivated by the Stefan problem for the evolution of the interface in a solid-liquid phase transition and in a Hele-Shaw type flow between two fluids with different viscosities.
Another typical assumption in the context of the Cahn-Hilliard equation is the mass conservation property that arises as a direct consequence of the standard no-flux boundary condition prescribed for the chemical potential associated with the phase field variable. While this condition is very natural for the engineering applications that Cahn and Hilliard had in mind originally (see [10]), the recent employment of the Cahn-Hilliard equation to describe other phenomena driven by phase segregation demands the incorporation of an external source term S in the model that reflects the fact that the system may not be isolated and the loss or production of mass is possible. Without claiming to be exhaustive, let us mention that numerous liquid-liquid phase segregation problems arise in cell biology [16] and in tumor growth models [20]. For this reason, we also included the presence of a source term in our investigation.
The standard isothermal Cahn-Hilliard system has been extensively studied in the past decades: see, e.g., [28] and the references therein. On the other hand, the mathematical understanding of nonisothermal Cahn-Hilliard systems is, thirty years after the seminal works by Alt and Pawlow (see [1,3] and, in particular, [2]) and twenty years after the groundbreaking work [17] by Gajewski for the nonlocal case, still far from being complete. Before presenting our system, let us discuss some recent literature. Concerning some analytic results of the aforementioned system by Caginalp, we mention the related contibutions [11,12,27]. Next, employing micro-force balance theory, Miranville and Schimperna proposed a further derivation in [29], and the well-posedness of a related system has been addressed in [26]. Moreover, we point out the recent contribution [15] by De Anna et al., where two new thermodynamically consistent models related to nonisothermal Cahn-Hilliard systems have been derived. Finally, we refer to [18,19] for some mathematical results on a relaxed version of the above systems endowed with dynamic boundary conditions. Motivated by the aforementioned remarks, we aim at analyzing a nonisothermal Cahn-Hilliard type system with source term in this paper. To this end, let Ω ⊂ R 3 be the spatial domain where the evolution takes place, and T > 0 a given final time. We then consider the following initial-boundary value problem: 3) ∂ n ϕ = ∂ n µ = ∂ n (κ 1 ∂ t w + κ 2 w) = 0 on Σ = ∂Ω × (0, T ), (1.4) ϕ(0) = ϕ 0 , w(0) = w 0 , ∂ t w(0) = w 1 in Ω. (1.5) In the above system, the unknowns have the following physical meaning: ϕ is a normalized difference between the volume fractions of pure phases in the binary mixture (the dimensionless order parameter of the phase transformation, which should attain its values in the interval [−1, 1]), µ is the associated chemical potential, and w is the so-called thermal displacement (or freezing index), which is directly connected to the temperature ϑ (which in the case of the Caginalp model is actually a temperature difference) through the relation w(·, t) = w 0 + F log (r) := (1 + r) ln(1 + r) + (1 − r) ln(1 − r) − c 1 r 2 if |r| ≤ 1, +∞ otherwise, (1.8) with the convention 0 ln(0) := lim rց0 r ln(r) = 0. In (1.8), c 1 > 1 so that F log is nonconvex. Another example is the double obstacle potential , where, with c 2 > 0, Singular potentials like (1.8) and (1.9) are difficult to handle from the mathematical viewpoint, but have the great advantage that if a solution exists, then it automatically inherits the property of being physically meaningful, that is, ϕ ∈ [−1, 1]. In general, this cannot be guaranteed for regular potentials like the quartic (1.7), which, in this sense, provides just an approximation of the more physical choices. In cases like (1.9), one has to split F into a nondifferentiable convex part β (the indicator function of [−1, 1] in the present example) and a smooth (usually quadratic) perturbation π. Accordingly, the second equation (1.2) has then to be understood as the differential inclusion or, equivalently, with the help of a selection ξ, as the identity µ = −∆ϕ + ξ + π(ϕ) + a − b∂ t w with ξ ∈ ∂ β(ϕ).
The above system is a formal extension of the Cahn-Hilliard system introduced by Caginalp in [8] (see also the derivation in [6,Ex. 4.4.2,(4.44), (4.46)]); it corresponds to the Allen-Cahn counterpart analyzed recently in [14]. The main differences between our system and the one originally introduced in [8] are the following: • In [8], we have a = λ (the specific latent heat).
• In [8], the heat flux is assumed in the standard Fourier form q = −κ 1 ∇ϑ, while we we follow the works by Green and Naghdi [23][24][25] and Podio-Guidugli [32] and postulate that (1.10) Note that this assumption accounts for a possible previous thermal history of the phenomenon. We also observe that the no-flux condition q · n = 0 then gives rise to the third boundary condition in (1.4).
• The third -and main -difference is that (1.1)-(1.2) comprises a Cahn-Hilliard system with a source term S := f − γϕ, which is independent of temperature.
The presence of S radically changes the behavior of the Cahn-Hilliard equation since the mass conservation property is no longer fulfilled. In fact, due to the no-flux boundary condition for µ in (1.4), a formal consideration, i.e., testing (1.1) by 1/|Ω|, readily reveals that the mass balance law of the order parameter ϕ is ruled by In this direction, we highlight that, especially when working with singular potential like (1.9), the control of the mean value of ϕ plays a crucial role in the mathematical analysis of Cahn-Hilliard-type systems. Besides, in the case when f is a positive constant such that f ∈ (−γ, γ), the equations (1.1)-(1.2) correspond to the so-called Cahn-Hilliard-Oono system, see, e.g., [22] and [13].
Finally, let us mention that the differential structure with respect to w in equation (1.3) is sometimes also referred to as the strongly damped wave equation (see, e.g., [31] and the references therein).
Let us conclude this section by presenting an outline of the paper. In the following section, we state the main results and list the corresponding assumptions. Then, from Section 3 onward, we start proving the mentioned results. In particular, Section 3 is devoted to showing some continuous dependence results enjoyed by the system (1.1)-(1.5). In Section 4, we then introduce and solve a preparatory approximating problem that will allow us to prove in Section 5 the existence of weak solutions, as well as some regularity results.

Statement of the problem and main results
Throughout the paper, Ω indicates a bounded and connected open subset of R 3 (the lower dimensional cases can be treated in the same way) with smooth boundary Γ := ∂Ω. In the following, |Ω| and n denote the Lebesgue measure of Ω and the outward unit normal vector field on Γ, respectively. Given a final time T > 0, we set Q t := Ω × (0, t) for t ∈ (0, T ] and Q := Q T . (2.1) Given a Banach space X, we denote its norm by · X , with the exceptions of L p spaces on Ω and Q, whose norms are denoted by · p for 1 ≤ p ≤ ∞ (if no confusion can arise), and of the space H introduced below. For brevity, we use the same symbol for the norm in a space and in any power thereof. Furthermore, for Banach spaces X and Y , we notice that the linear space X ∩ Y becomes a Banach space when equipped with its natural graph norm Then, we introduce the shorthands and endow these spaces with their natural norms. For simplicity, we write · instead of · H . Moreover, we denote by ( · , · ) and · , · the standard inner product of H and the duality pairing between the dual space V * of V and V itself, respectively. We identify H with a subspace of V * in the usual way, i.e., in order that This makes (V, H, V * ) a Hilbert triplet. We notice that all of the embeddings W ֒→ V ֒→ H ֒→ V * are dense and compact. Next, we define the generalized mean value v of a generic element v ∈ V * by setting v : where we have written 1 for the constant function that takes the value 1 in Ω. It is clear that v reduces to the usual mean value if v ∈ H. The same notation v is employed also if v is a time-dependent function.
Let us come to the structural assumptions we make for our analysis. First, γ, a, b, κ 1 , κ 2 and λ are positive constants. (2.4) Next, in order to allow for general double-well potentials in (1.2), we assume that F : R → (0, +∞] admits the decomposition F = β + π, where (2.5) β : R → [0, +∞] is convex, l.s.c., and fulfills β(0) = 0, (2.6) π ∈ C 1 (R), and its derivative is Lipschitz continuous. (2.7) Moreover, we set β := ∂ β and π := π ′ , (2.8) where ∂ denotes the subdifferential operator, and notice that β : R → 2 R is maximal monotone with corresponding domain D(β) and that 0 ∈ β(0). We observe that all of the examples (1.7)-(1.9) of potentials introduced before do satisfy the conditions required above. Of course, in the case of nonregular potentials like the double obstacle (1.9), the second equation (1.2) has to be intended as the differential inclusion or, equivalently, with the help of a selection ξ ∈ β(ϕ) a.e. in Q, as the identity As for the data, we assume that f ∈ L ∞ (Q) and g ∈ L 2 (0, T ; H), (2.9) ϕ 0 ∈ W, w 0 ∈ V and w 1 ∈ H . (2.10) However, we also need some compatibility conditions between the data f and ϕ 0 and the domain of β. These are in fact already expected, as we are dealing with possible singular potentials and a mass source. In particular, let us repeat that the contribution S := f − γϕ in (1.1) plays the role of a (phase-dependent) mass source/sink in the model. Indeed, by formally testing (1.1) by 1/|Ω|, and using (1.4), we infer that the mass balance law of the system reads d dt Therefore, it is natural to expect to have some compatibility conditions between the structure of the source term S, thus on the constant γ and the function f , and the possibly singular potential β. Namely, setting ρ := f ∞ γ , (2.11) and noting that ϕ 0 ∈ C 0 (Ω) by (2.10), we require that all of the quatities where ( · ) − and ( · ) + denote the negative and positive part functions, respectively.
Remark 2.1. The assumptions on f and ϕ 0 can be weakened slightly. However, in doing so, we would have to replace (2.11)-(2.12) by more complicated compatibility conditions. Moreover, when regularizing our problem as we are going to do in the forthcoming Section 4, we would have to regularize ϕ 0 as well. This would lead to estimates depending on the regularization parameter, so that further uniform estimates had to be performed.
At this point, we can rigorously state our notion of (weak) solution to the aforementioned problem under study. A weak solution to the system (1.1)-(1.5) is a quadruplet (ϕ, µ, ξ, w) enjoying the regularity properties (2.16) and satisfying for every v ∈ V and a.e. in (0, T ) , (2.17) µ = −∆ϕ + ξ + π(ϕ) + a − b∂ t w and ξ ∈ β(ϕ) a.e. in Q , (2.18) for every v ∈ V and a.e. in (0, T ) , (2.20) The present paper is devoted to the study of the well-posedness of the above problem and of the regularity of its solutions. Our first result is an existence theorem. ϕ ∈ L 2 (0, T ; W 2,6 (Ω)) and ξ ∈ L 2 (0, T ; L 6 (Ω)) , (2.21) as well as the estimate with a positive constant K 1 that depends only on the structure of the system, Ω, T , and upper bounds for the norms of the data and the quantities related to assumptions (2.9)-(2.12).
Uniqueness cannot be expected, in general, as it usually occurs in Cahn-Hilliard type problems with nonregular potentials. However, we have the result stated below, which ensures continuous dependence on f and g for the components ϕ and w of every solution with fixed initial data. In the statement, we use the following notation for convolution products with 1: (2.23) holds true with a positive constant K 2 that depends only on the structure of the system, Ω, T , and an upper bound for the norms of ξ 1 and ξ 2 in L 1 (Q).
Partial uniqueness in general and full uniqueness if β is single valued trivially follow, as stated below.

27)
and satisfies the estimate In particular, if (ϕ, µ, ξ, w) is a solution in the sense of Theorem 2.5, then h belongs to L 2 (0, T ; H), and the elliptic regularity theory yields that where C Ω depends only on Ω. Thus, the same norm can be estimated by a constant which is proportional to K 3 . By solving κ 1 ∂ t w + κ 2 w = u for w, we obtain that w ∈ H 1 (0, T ; H 2 (Ω)) or w ∈ H 1 (0, T ; W ), provided that w 0 ∈ H 2 (Ω) or w 0 ∈ W , respectively.
In the next sections, when proving our results, we widely use Hölder's inequality, as well as the Young, Poincaré, Sobolev and compactness inequalities recalled below: Here, C Ω is a constant that depends only on Ω, while C Ω,p,δ depends on p and δ, in addition. Moreover, the symbol · * denotes the norm in V * defined by the forthcoming formula (2.35). The last two inequalities are related to the (three-dimensional) embedding V ֒→ L p (Ω) which holds for p ∈ [1, 6] and is compact if p < 6.
Finally, we take advantage of a tool that is commonly used in the study of problems related to the Cahn-Hilliard type equations: consider, for ψ ∈ V * , the problem of finding (2.33) Obviously, if ψ ∈ H, then this problem is just the usual homogeneous Neumann problem for the Poisson equation −∆u = ψ. Now, since Ω is connected, for ψ ∈ V * , (2.33) is solvable if and only if ψ has zero mean value. Moreover, if this condition is satisfied, then there exists a unique solution possessing zero mean value. This entails that the operator is well defined and yields an isomorphism between the above spaces. Besides, it follows that the formula defines a Hilbert norm in V * that is equivalent to the standard dual norm. From the above definitions one trivially derives that Ω ∇Nψ · ∇v = ψ, v for every ψ ∈ dom(N) and v ∈ V , We conclude this section by stating a general rule concerning the constants that appear in the estimates to be performed in the following. The small-case symbol c stands for a generic constant whose actual value may change from line to line, and even within the same line, and depends only on Ω, the shape of the nonlinearities, and the constants and the norms of the functions involved in the assumptions of the statements. In particular, the values of c do not depend on the parameters ε > 0 and n ∈ N that will be introduced in the next sections. A small-case symbol with a subscript like c δ indicates that the constant may depend on the parameter δ, in addition. On the contrary, we mark precise constants that we can refer to by using different symbols (see, e.g., (2.31)).

Continuous dependence
This section is devoted to the proof of Theorem 2.3. Let f i and g i , i = 1, 2, satisfy (2.9), and let (ϕ i , µ i , ξ i , w i ) be any two corresponding solutions as in the statement. We set, for convenience, ϕ := ϕ 1 − ϕ 2 , and define µ, ξ, w, f and g analogously. We first make some preliminary observations. Recalling (2.13) (see (2.2) for the definition of W ) and testing (2.17) by 1/|Ω|, we find that Then, on the one hand, by multiplying this equality by Ω v, we deduce that On the other hand, by (formally) multiplying (3.1) by sign(ϕ), where sign : R → R is the sign function defined by sign(r) := r/|r| if r = 0 and sign(0) = 0, we infer that We now start the proof of the theorem. We use the properties (2.34)-(2.39) of the operator N and recall the notation (2.23) for convolution products with 1. We write equation (2.17) for both solutions and take the difference, obtaining an equality from which we subtract (3.2) to arrive at the identity for every v ∈ V and a.e. in (0, T ). Since (ϕ−ϕ)(t) has zero mean value for every t ∈ [0, T ], we are allowed to test the above equation by N(ϕ − ϕ). Integration with respect to time then leads to, for every t ∈ [0, T ], Next, we write (2.18) for both solutions, multiply the difference by −(ϕ−ϕ), and integrate over Q t , finding that Finally, we write (2.19) for both solutions and take the convolution with 1. Then, we test the difference of the corresponding equalities by At this point, we add (3.4)-(3.6) to each other and notice that some cancellations occur. Moreover, β is monotone, and thus all of the remaining terms on the left-hand side are nonnegative. We treat those on the right-hand side individually. First, we have that where we have used the trivial inequalities v * ≤ c |v| ≤ c v * , which hold for every v ∈ V * .
Next, we fix a constant M such that ξ i L 1 (Q) ≤ M for i = 1, 2. Then, recalling (3.3), we have that Also, in view of the Lipschitz continuity of π, the obvious inequality v ≤ v for v ∈ H, and the compactness inequality (2.32), we find that Moreover, by Young's inequality and arguing as above, we have that on account of (3.3). We deal with the next integral using integration by parts to infer that In addition, it is clear that Finally, we note that Upon collecting (3.4)-(3.6) and the inequalities shown above, we obtain that where c has the dependence required for the constant K 2 in the statement of the theorem. On the other hand, (3.3) implies that By combining this with the previous inequality, we are in a position to apply Gronwall's lemma and obtain the desired estimate (2.24), which concludes the proof.

Approximation
In this section, we introduce and solve a proper approximating problem depending on the parameter ε ∈ (0, 1). First of all, we replace the functional β and the maximal monotone graph β by their Moreau-Yosida regularizations β ε and β ε , respectively (see, e.g., [5, pp. 28 and 39]). We recall that 0 ≤ β ε (r) = r 0 β ε (s) ds ≤ β(r) for every r ∈ R , (4.1) β ε is monotone and Lipschitz continuous with β ε (0) = 0 , where β • (r) denotes the element of the section β(r) having minimum modulus. The approximating problem to be considered consists in finding a triplet (ϕ ε , µ ε , w ε ) satisfying the regularity properties and solving the following system of variational identities or equations and initial conditions: for every v ∈ V and a.e. in (0, T ) , (4.9) ϕ ε (0) = ϕ 0 , w ε (0) = w 0 and ∂ t w ε (0) = w 1 . (4.10) We remark that here we obviously do not need to consider any selection ξ as β ε is regular and single valued. Here is our basic result. The rest of this section is devoted to the proof of the above theorem. Clearly, uniqueness is a consequence of Theorem 2.3, since β ε satisfies all the assumptions postulated for β in (2.5)-(2.7), and it is single valued, in addition (cf. Corollary 2.4).
To prove the existence of a solution, we start from a Faedo-Galerkin scheme. To this end, we introduce the nondecreasing (ordered) sequence {λ j } of eigenvalues and the corresponding complete orthonormal sequence {e j } of eigenfunctions of the eigenvalue problem for the Laplace operator with homogeneous Neumann boundary conditions. Namely, we have that −∆e j = λ j e j in Ω and ∂ n e j = 0 on Γ for j = 1, 2, . . . , (4.11) Ω e i e j = δ ij for every i and j, (4.12) with the standard Kronecker symbols δ ij . Moreover, we set V n := span{e j : 1 ≤ j ≤ n} for n = 1, 2, . . . (4.13) and recall that the union of these spaces is dense in both V and H. Notice that all of the eigenfunctions are smooth since Ω is smooth. Furthermore, as Ω is connected, we have that λ 1 = 0 < λ 2 , and V 1 is the subspace of constant functions.
The discrete problem consists then in finding a triplet (ϕ n , µ n , w n ) of functions satisfying ϕ n ∈ H 1 (0, T ; V n ) , µ n ∈ L 2 (0, T ; V n ) and w n ∈ H 2 (0, T ; V n ) , (4.14) and solving the discrete problem for every v ∈ V n and a.e. in (0, T ) , for every v ∈ V n and a.e. in (0, T ) , The strategy of the proof can be schematized as follows. First, we show that the above problem has a unique solution. Then, we perform a number of a priori estimates that allow us to pass to the limit as n tends to infinity. In this way, we identify a limit triple (ϕ ε , µ ε , w ε ), which then is shown to be a solution to the problem (4.7)-(4.10) enjoying the desired regularity properties.
In terms of these true unknowns the equations (4.15)-(4.17) take the form where the matrix A = (A ij ) n i,j=1 and the vectors f = (f i ) n i=1 and g = (g i ) n i=1 are given by while F ε : R n → R n is the function whose i-th component (i = 1, . . . , n) is given by Clearly, f and g are L 2 functions and F ε is Lipschitz continuous. Moreover, the initial conditions (4.18) provide initial conditions for the vectors ϕ n , w n and w ′ n . We first eliminate ϕ ′ n from (4.21) by exploiting (4.19) and then eliminate every occurrence of µ n by means of (4.20). In this way, we obtain a well-posed Cauchy problem for the pair ( ϕ n , w n ) coupled with the chemical potential equation (4.20), and it is clear that the new problem is equivalent to the previous one. Hence, we find a unique solution with the regularity ϕ n ∈ H 1 (0, T ; R n ), w n ∈ H 2 (0, T ; R n ), and µ n ∈ L 2 (0, T ; R n ), so that the discrete problem has a unique solution, as claimed.
A preliminary estimate. We recall that V n ⊃ V 1 and that V 1 is the subspace of constant functions. Hence, we can test (4.15) by 1/|Ω| to obtain that Before continuing, it is worth making some observations on projections which are collected in the following remark.
Remark 4.2. Let P n : H → V n be the H-orthogonal projection operator. We list some inequalities that hold true for every n ∈ N, as well as convergence properties as n tends to infinity. For every v ∈ H, we clearly have that P n v ≤ v , and P n v → v strongly in H.
Assume now that v ∈ V . Then it is easy to see that also For a detailed proof, see, e.g., [13,Rem. 4.2]. In particular, we deduce that Next, assume that v ∈ W . Then, we have that (v, e j )λ j e j .
We deduce that ∆P n v = P n ∆v, and we can apply the above inequalities and convergence properties to ∆v as well in order to recover further information on P n v. We obtain, with a constant C Ω that depends only on Ω, that Notice that all this can be applied to the initial values of the discrete solution as they are projections on V n . Now, we consider time-dependent functions. A simple combination of the above properties with the Lebesgue dominated convergence theorem shows the following: if we assume that v ∈ L 2 (0, T ; H) or v ∈ L 2 (0, T ; V ) and define v n by setting v n (t) := P n (v(t)) for a.a. t ∈ (0, T ), then v n → v strongly in L 2 (0, T ; H) or L 2 (0, T ; V ), respectively.
At this point, we can start estimating, and we recall that the symbol c stands for possibly different constants independent of ε and n according to our general rule regarding constants stated at the end of Section 2. We repeatedly owe to the properties (2.34)-(2.39) related to the operator N without further reference.
First uniform estimate. We first observe that Nv ∈ V n for every v ∈ V n satisfying v = 0. Indeed, both v and w := Nv can be expressed in terms of the eigenfunctions e j , and we have that (v, e j ) e j .
Hence, (w, e j ) = 0 for every j > n (since λ j > 0 for j > 1), i.e., w ∈ V n . Once this is established, we take the difference between (4.15) written for a generic v ∈ V n and (4.23) multiplied by Ω v, write the resulting equality at the time s ∈ (0, T ) and test it by N(ϕ n − ϕ n )(s). Then, by integrating over (0, t) ⊂ (0, T ), we obtain that At the same time, we test (4.16), written at the time s, by −(ϕ n (s) − ϕ n (s)) and integrate over (0, t). It results that Finally, we take the convolution between (4.17) and 1 (see (2.23)) and test the resulting equality by ∂ t w n . After time integration, we obtain that At this point, we add (4.25)-(4.27) to each other and notice that a cancellation occurs. After rearranging, we deduce that The integral involving β ε can be estimated from below by combining (4.22) and (4.24) as follows: All of the other terms on the left-hand side are nonnegative. For those on the right-hand side, we perform separate estimates.
Since the embedding H ֒→ V * is continuous, the first term I 1 is uniformly bounded by the assumption (2.10) on ϕ 0 , Remark 4.2, and estimate (4.24). Next, we have that Owing to Young's inequality, the Lipschitz continuity of π, and (4.24) once more, we have, for every δ > 0, where in the second line we also used the compactness inequality (2.32). Next, arguing similarly, we obtain that thanks to (4.24) and to our assumptions on the initial data w 1 and ϕ 0 (by applying Remark 4.2) and on g. The last term to be estimated is first treated by an integration by parts. Finally, by also using Young's inequality and the estimate for ∇w n (0) obtained by applying Remark 4.2, we have, for every δ > 0, that On the other hand, we also have that At this point, we recall (4.28) and all the above estimates, choose δ small enough, and apply Gronwall's lemma. We obtain that ϕ n − ϕ n L ∞ (0,T ;V * )∩L 2 (0,T ;V ) + β ε (ϕ n ) L 1 (Q) + ∂ t w n L 2 (0,T ;H) + ∇(w n − w n (0)) L ∞ (0,T ;H) ≤ c , Consequence. By testing (4.16) by 1/|Ω|, and owing to (4.29), we infer that Second uniform estimate. We test the equations (4.15), (4.16), and (4.17), by µ n , −∂ t ϕ n , and (b/λ)∂ t w n , respectively, sum up and notice that the terms involving the products ∂ t ϕ n µ n and ∂ t w n ∂ t ϕ n cancel each other. Then, we integrate in time and rearrange to obtain that where all of the terms on the left-hand side are nonnegative. Moreover, as before, we can recall Remark 4.2 in order to estimate the terms involving the initial data, and just the one containing β ε needs further comments. Since ϕ 0 belongs to W by (2.10), ϕ n (0) converges to ϕ 0 strongly in W , hence uniformly. On the other hand, by the quoted assumption, min ϕ 0 and max ϕ 0 belong to the interior of D(β). Thus, for some n 0 and every n ≥ n 0 , all of the values of ϕ n (0) belong to a compact interval I contained in the interior of D(β). By also recalling (4.1), we thus may conclude that Ω β ε (ϕ n (0)) ≤ Ω β(ϕ n (0)) ≤ max r∈I β(r) = c .
It is understood that n ≥ n 0 from now on, which is no restriction since we aim at letting n tend to infinity eventually. Let us come to the other terms on the right-hand side. The last one can be dealt with employing Young's inequality (and then Gronwall's lemma), and the integral that precedes it has already been estimated, since it is a multiple of the mean value. Moreover, since π grows at most quadratically by condition (2.7), we can infer from the compactness inequality (2.32) and (4.29) that The other integrals that need some treatment are those containing µ n . We have that where we have used the Poincaré inequality (2.30), our assumptions on f (see (2.9)), (4.29), and (4.30). By coming back to (4.31), collecting the above estimates and observations, and applying the Gronwall lemma, we conclude that ∇µ n L 2 (0,T ;H) + ϕ n L ∞ (0,T ;V ) + β ε (ϕ n ) L ∞ (0,T ;L 1 (Ω)) + w n W 1,∞ (0,T ;H)∩H 1 (0,T ;V ) ≤ c . (4.32) Third uniform estimate. Next, recalling that every constant is allowed as a test function, we test (4.16) by ϕ n (t) − ϕ n (t) and rearrange. Omitting the time variable for brevity, we have a.e. in (0, T ) that In view of (4.24) and of our assumption (2.12), we can bound the integral involving β ε from below using (4.22): As for the right-hand side, the first term needs some treatment. Thanks to Poincaré's inequality (2.30) and to (4.32), we have that Ω µ n (ϕ n − ϕ n ) = Ω (µ n − µ n )(ϕ n − ϕ n ) ≤ c ∇µ n ϕ n − ϕ n ≤ c ∇µ n .
The sum of the other terms is bounded from above by c ϕ n 2 + |ϕ n | 2 + ∂ t w n 2 + c .
Combining (4.33), the inequalities just derived, and the previous estimates, we see that the function t → Ω |β ε (ϕ n (t))| is bounded from above by an L 2 (0, T ) function independently of both n and ε, that is, it holds whence we trivially derive an estimate in L 2 (0, T ) for the mean value of β ε (ϕ n ). Then, from (4.16), we can estimate the L 2 (0, T ) norm of µ n . This, (4.32), and the use of the Poincaré inequality once more, imply that µ n L 2 (0,T ;V ) ≤ c . Fourth uniform estimate. We recall Remark 4.2 and use the notations introduced there. We fix some v ∈ L 2 (0, T ; V ) and define v n ∈ L 2 (0, T ; V n ) by setting v n (t) = P n (v(t)) for a.a. t ∈ (0, T ). Then, we test (4.15) by v n , and integrate over time to obtain that On the other hand, we have that since ∂ t ϕ n is V n -valued. Thus, we readily conclude that Fifth uniform estimate. Thanks to (4.36), the same argument, applied to equation (4.17) for w n , yields that (4.37) Passage to the limit. At this point, we can pass to the limit as n → ∞. Indeed, by recalling (4.29), (4.32) and (4.35)-(4.37), and applying well-known weak, weak star, and strong compactness results (for the latter see, e.g., [33,Sect. 8,Cor. 4]), we deduce that there exists a triple (ϕ ε , µ ε , w ε ) such that as n tends to infinity (at least for a subsequence which is not relabeled). Moreover, since β ε and π are Lipschitz continuous, we also have that β ε (ϕ n ) → β ε (ϕ ε ) and π(ϕ n ) → π(ϕ ε ) strongly in C 0 ([0, T ]; H). We claim that this triple is a (weak) solution to problem (4.7)-(4.10). Since ϕ n (0), w n (0) and ∂ t w n (0) are the H projections of ϕ 0 , w 0 and w 1 , they strongly converge in H to ϕ 0 , w 0 and w 1 , respectively. On the other hand, they converge to ϕ ε (0), w ε (0) and ∂ t w ε (0), respectively, strongly (at least) in V * , thanks to (4.38) and (4.40). Hence, the initial conditions (4.10) are satisfied. Now, we show that the variational equations (4.7)-(4.9) are satisfied as well. We recall Remark 4.2, fix any v ∈ L 2 (0, T ; V ), define v n ∈ L 2 (0, T ; V n ) by setting v n (t) := P n (v(t)) for a.a. t ∈ (0, T ), and observe that v n converges to v strongly in L 2 (0, T ; V ). Next, we test each of the equations (4.15)-(4.17) by v n and integrate in time over (0, T ). At this point, on account of the convergence properties proved or mentioned, it is straightforward to pass to the limit as n → ∞ in the equalities we obtain. The resulting equalities are the same equations with (ϕ ε , µ ε , w ε ) in place of (ϕ n , µ n , w n ), i.e., the time-integrated versions of (4.15)-(4.17) with arbitrary time-dependent test functions v ∈ L 2 (0, T ; V ), which are equivalent to (4.15)-(4.17) themselves.

Existence and regularity
This final part of the paper is devoted to prove the existence and regularity results stated in Theorems 2.2 and 2.5.

Proof of Theorem 2.5
Following the line of arguments of the proof of Theorem 2.2, we use the estimates already established for the discrete solution (ϕ n , µ n , w n ) and the approximating solution (ϕ ε , µ ε , w ε ) and perform further estimates. So, we want to show that ϕ n H 1 (0,T ;V ) + µ n L ∞ (0,T ;V ) + w n H 2 (0,T ;H)∩W 1,∞ (0,T ;V ) ≤ c , (5.6) at least for sufficiently large n ∈ N, as well as β ε (ϕ ε ) L ∞ (0,T ;L 6 (Ω)) ≤ c and ϕ ε L ∞ (0,T ;W 2,6 (Ω)) ≤ c , with a constant c that has the same dependence as the constant K 3 in the statement. To prove (5.6), we first observe that the component µ n of the discrete solution (ϕ n , µ n , w n ) is more regular than required: indeed, it is Lipschitz continuous, as follows from looking at µ n in equation (4.20). Hence, we are allowed to take t = 0 in (4.16). It results that µ n (0) = P n (−∆ϕ 0 + β ε (ϕ n (0)) + π(ϕ n (0)) + a − bw 1 ) , where P n : H → V n is the orthogonal projection operator. Recall that µ n (0) also depends on ε, of course, despite of the used notation. It is convenient to first establish an estimate for µ n (0). holds true for every ε ∈ (0, 1) and every n ≥ n 0 .
In particular, the sequence {β ε (ϕ n (0))} is uniformly bounded in H. On the other hand, since the restriction of β to the interior of D(β) is a C 1 function by (2.25), the following inequality holds: |β ′ (s)| =: C for every r ∈ [r * , r * ].
For a detailed proof (with a different notation) see, e.g., [13, formula (5.2)]. Then, we have that ∇β ε (ϕ n (0)) = β ′ ε (ϕ n (0))∇ϕ n (0) ≤ C ∇ϕ 0 , so that (5.10) follows. At this point, we easily derive (5.9). By also accounting for assumption (2.26) and Remark 4.2 once more, we have indeed Let us now continue with the proof. It is understood that n ≥ n 0 (given by the lemma) from now on. In order to make the argument more transparent, it is convenient to prepare an auxiliary estimate depending on a positive parameter M whose value will be chosen later on.
Auxiliary estimate. We repeat part of the argument used to arrive at (4.34), but this time we avoid time integration. We account for (4.24) in order to apply (4.22) once more. We test (4.16) a.e. in (0, T ) by M(ϕ n − ϕ n ). Then, we invoke the Poincaré inequality (2.30) and the Young inequality (2.29) with δ = (8MC Ω ) −1 . By also taking advantage of (4.32), we find (a.e. in (0, T )) that with the positive constant C 0 arising from (4.22). Since (4.16) and the already known estimates for ϕ n and ∂ t w n imply that |µ n | ≤ |β ε (ϕ n )| + c a.e. in (0, T ), we deduce that δ 0 M|Ω||µ n (t)| ≤ 1 8 Ω |∇µ n (t)| 2 + c M for a.a. t ∈ (0, T ). (5.11) Sixth uniform estimate. By virtue of the already proved regularity of µ n , we can now take ∂ t µ n as a test function in (4.15) and, at the same time, we can differentiate (4.16) with respect to time and then test the resulting equality by −∂ t ϕ n . We do this and also test (4.17) by (b/λ)∂ 2 t w n . Then, we sum up and notice that four terms cancel each other. Finally, we integrate with respect to time and add (5.11) to the resulting equality. Collecting the terms, we obtain that 2λ Ω |∇∂ t w n (0)| 2 − bκ 2 λ Qt ∇w n · ∇∂ 2 t w n + b λ Qt g ∂ 2 t w n + 1 8 Ω |∇µ n (t)| 2 + c M , (5.12) where all of the terms on the left-hand side are nonnegative. The first term on the righthand side is estimated by the above lemma. The other term involving an initial value is bounded by the V -norm of w 1 . As for the first volume term on the right-hand side, we integrate by parts in time and have that The volume integral on the right is estimated by ∂ t f − γ∂ t ϕ n L 2 (0,T ;V * ) µ n L 2 (0,T ;V ) , which is bounded on account of (2.26), (4.35) and (4.36). The last term is easily treated once again with the help of the lemma. The remaining term is dealt with by using the Young and Poincaré inequalities: where we have used the special symbol C * to mark the constant in front of |µ n (t)| for future reference. Notice that C * is a multiple of an upper bound for the norm of f −γϕ n in L ∞ (0, T ; V * ), which is known by (2.26) and (4.29). Next, it turns out that − Qt π ′ (ϕ n )|∂ t ϕ n | 2 ≤ 1 2 Qt |∇∂ t ϕ n | 2 + c ∂ t ϕ n 2 L 2 (0,T ;V * ) ≤ 1 2 Qt |∇∂ t ϕ n | 2 + c , thanks to the Lipschitz continuity of π, the compactness inequality (2.32), and (4.36).
It remains to estimate the volume integrals involving w n that appear on the righthand side of (5.12). The last one is trivially treated via Young's inequality. The other can be dealt with as follows, using integration by parts, (4.32), the Young inequality, and Remark 4.2. Indeed, we have that − bκ 2 λ Qt ∇w n · ∇∂ 2 t w n = bκ 2 λ Qt |∇∂ t w n | 2 − bκ 2 λ Ω ∇w n (t) · ∇∂ t w n (t) + bκ 2 λ Ω ∇w n (0) · ∇∂ t w n (0) At this point, we can easily conclude. Indeed, if we choose M in order that δ 0 M|Ω| = C * + 1, and collect (5.12) and all the above estimates, then we obtain ϕ n H 1 (0,T ;V ) + µ n L ∞ (0,T ;V ) + w n H 2 (0,T ;H)∩W 1,∞ (0,T ;V ) ≤ c (5.14) with a constant c that has the same dependence on the structure and the data as required, since even M has this property.
Conclusion of the proof. We are then left with checking (5.7). To this end, it suffices to come back to the nonlinear elliptic problem (4.42) and the corresponding estimates (4.43)-(4.44) and to argue as we did to prove (4.45)-(4.46), in this case avoiding time integration.
At this point, we can easily conclude. As the discrete solution (ϕ n , µ n , w n ) converges as n → ∞ to the solution (ϕ ε , µ ε , w ε ) to the approximating problem, it is clear that the analogue of (5.6) for (ϕ ε , µ ε , w ε ) holds true with the same constant, by the semicontinuity of norms. We conclude that the convergence properties (5.2)-(5.5) can be improved on account of the estimate just mentioned and (5.7). On the other hand, the previous proof ensures that the limiting quadruplet (ϕ, µ, ξ, w) is a solution to problem (2.17)-(2.20), and the estimates proved for the approximating solution are conserved in the limit. Therefore, the proof of Theorem 2.5 is complete.
other support from the MIUR-PRIN Grant 2020F3NCPX "Mathematics for industry 4.0 (Math4I4)" and their affiliation to the GNAMPA (Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni) of INdAM (Istituto Nazionale di Alta Matematica).