Deposited solid thin films are unstable and could dewet to form isolated islands on the substrate in order to minimize the total surface energy (Leroy et al.2016; Thompson2012). This phenomenon is known as solid-state dewetting (SSD), since the thin films remain in a solid state during the process. SSD has attracted a lot of attention recently and is emerging as a promising route to produce patterns of arrays of particles used in sensor technology, optical and magnetic devices, and catalyst formations, see, e.g., Armelao et al. (2006), Benkouider et al. (2015), Schmidt et al. (2009), Bollani et al. (2019), Salvalaglio et al. (2020), Backofen et al. (2017).

The dominant mass transport mechanism in SSD is surface diffusion (Srolovitz and Safran1986). This evolution law was first introduced by Mullins (1957) to describe the mass diffusion within interfaces in polycrystalline materials. For surface diffusion, the normal velocity of the interface is proportional to the surface Laplacian of the mean curvature. In the case of SSD, the evolution of the interface that separates the thin film from the surrounding vapor also involves the motion of the contact line, i.e., the region where the film/vapor interface meets the substrate. The equilibrium contact angle is given by Young’s law which prescribes a force balance along the substrate. Many efforts have been devoted to SSD problems in recent years. For example, a large body of experiments have revealed that the pattern formations could depend highly on the crystallographic alignments, the film sizes and shapes, as well as the substrate topology, see, e.g., Ye and Thompson (2011), Amram et al. (2012), Thompson (2012), Naffouti et al. (2017), Bollani et al. (2019). In addition, mathematical studies based on different models have been considered in Srolovitz and Safran (1986), Burger (2005), Dornel et al. (2006), Fonseca et al. (2012), Jiang et al. (2012), Naffouti et al. (2017), Boccardo et al. (2022), Jiang and Zhao (2019), Wang et al. (Jan 2015), Jiang et al. (2020), Dziwnik et al. (2017).

In this work, we aim to study the SSD problem with anisotropic surface energies in the diffuse-interface framework. In the isotropic case, diffuse-interface models are based on the Ginzburg–Landau energy

where\(\Omega \subset {\mathbb {R}}^d\) is a given domain with\(d\in \{2,3\}\),\(\varphi :\Omega \rightarrow {\mathbb {R}}\) is the order parameter,\(\varepsilon >0\) is a small parameter proportional to the thickness of the interfacial layer, and\(F(\varphi )\) is the free energy density. The following three choices forF are mainly used in the literature:

1. (i)

   The smooth double-well potential (Taylor and Cahn1994)

   $$\begin{aligned} F(\varphi )=\frac{1}{2}(1-\varphi ^2)^2, \end{aligned}$$

   (1.2a)

   which has two global minimum points at\(\varphi =\pm 1\) and a local maximum point at\(\varphi =0\);

2. (ii)

   The logarithmic potential (Cahn and Hilliard1958)

   $$\begin{aligned} F(\varphi )=\frac{1}{2}\theta \, [(1+\varphi )\,\log (1+\varphi ) +(1-\varphi )\,\log (1-\varphi )]+\frac{1}{2}(1-\varphi ^2),\nonumber \\ \end{aligned}$$

   (1.2b)

   where\(\theta >0\) is the absolute temperature. This potential has two minima\(\varphi = \pm (1-{\tilde{k}}(\theta ))\), where\({\tilde{k}}(\theta )\) is a small positive real number satisfying\({\tilde{k}}(\theta )\rightarrow 0\) as\(\theta \rightarrow 0\), and its usage enforces\(\varphi \) to attain values within\((-1,1)\);

3. (iii)

   The double-obstacle potential (Blowey and Elliott1991)

   $$\begin{aligned} F(\varphi )=\left\{ \begin{array}{ll} \frac{1}{2}(1-\varphi ^2) &{}\quad \text{ if }\quad |\varphi |\le 1,\\ \infty &{}\quad \text{ otherwise }. \end{array}\right. \end{aligned}$$

   (1.2c)

   It can be characterized via thedeep quench limit of the logarithmic potential, i.e., the limit of (1.2b) as\(\theta \rightarrow 0\).

The (isotropic) Cahn–Hilliard equation can be interpreted as a weighted\(H^{-1}\)-gradient flow of the free energy (1.1). It reads as

where\(m(\varphi )\) is a mobility function, together with Neumann boundary conditions for\(\mu \) and\(\varphi \). The Cahn–Hilliard equation was first introduced to study the spinodal decomposition in binary alloys (Cahn and Hilliard1958; Cahn1961) and has since then been used to model many other phenomenon, e.g., Abels et al. (2012), Garcke and Novick-Cohen (2000), Khain and Sander (2008), Bertozzi et al. (2006). We note that the double-obstacle potential is not differentiable at\(\varphi =\pm 1\), and the definition of the generalized chemical potential in this case becomes

where\(\partial F(\varphi )\) is the Fréchet sub-differential ofF at\(\varphi \) and\(\Delta \varphi \) has to be understood in a weak sense, see Blowey and Elliott (1991), Barrett et al. (2013). In the case of a constant mobility\(m(\varphi )\equiv 1\), (1.3) converges to the Mullins–Sekerka problem (Mullins and Sekerka1963) as\(\varepsilon \rightarrow 0\) (Pego1989; Alikakos et al.1994). In order to obtain the surface diffusion equation in the sharp-interface limit, a degenerate mobility needs to be chosen. For example, it was shown in Cahn et al. (1996) by a formal asymptotic analysis that the surface diffusion flow is recovered by considering a slow time scale\(\tau =O(\varepsilon ^{-1}t)\) of (1.3) with\(m(\varphi )= (1-\varphi ^2)_+\) and with the potential\(F(\varphi )\) either chosen as in (1.2c), or as in (1.2b) with\(\theta = O(\varepsilon ^\xi )\),\(\xi >0\). When using the smooth double-well potential (1.2a), the situation is less clear. While the limiting motion of surface diffusion is obtained with the choice\(m(\varphi )= (1-\varphi ^2)^2\) (Voigt2016; Rätz et al.2006; Jiang et al.2012; Dai and Du2014), using the less degenerate mobility\(m(\varphi )= (1-\varphi ^2)_+\) may not lead to pure surface diffusion in the limit\(\varepsilon \rightarrow 0\), since an additional bulk diffusion term is conjectured to be present due to the nonzero-flux contributions (Dai and Du2014; Lee et al.2015,2016). However, in all these cases, no rigorous proof for the sharp-interface limit or the presence of nonzero-flux contributions is available so far.

A natural generalization of the free energy (1.1) to the case of anisotropic surface energies is given by

see, e.g., Kobayashi (1993), Elliott (1997). Here,\(\gamma : {\mathbb {R}}^d\rightarrow [0, \infty )\) is the anisotropic density function, which is positively homogeneous of degree one, and\(A: = \frac{1}{2}\gamma ^2\). This then gives rise to the anisotropic Cahn–Hilliard equation

where\(A^\prime \) represents the gradient of the map\(A:{\mathbb {R}}^d\rightarrow [0, \infty )\). In contrast to the isotropic case, diffuse-interface models based on (1.5) result in a nonuniform asymptotic interface thickness, which in fact depends on the anisotropic density function\(\gamma (\nabla \varphi )\), see Wheeler and McFadden (1996), Wheeler (2006), Bellettini and Paolini (1996), Elliott and Schätzle (1996), Alfaro et al. (2010). To remedy this issue, an alternative energy of the form

can be considered, see Torabi et al. (2009), Salvalaglio et al. (2015), so that a constant thickness of the asymptotic interface is achieved. However, the resulting diffuse-interface models based on (1.7) become more nonlinear and are singular at\(\nabla \varphi =0\), which poses great challenges in the mathematical analysis and the stable numerical approximation. Therefore, in this work, we will restrict ourselves to the classical energy in (1.5). We also note that to guarantee that (1.6) converges to the anisotropic surface diffusion flow as\(\varepsilon \rightarrow 0\), a rescaled anisotropic coefficient needs to be introduced to the degenerate mobility (Rätz et al.2006; Li et al.2009). We refer to Sect. 2 below for the precise details.

Sketch of the structure for SSD near the contact line (green point), where the vapor, film and substrate phases meet (Color figure online)

Full size image

When it comes to SSD, as shown in Fig. 1, the total surface energy of the system consists of the film/vapor interface energy\({\mathcal {E}}_{\mathrm{inf}}\) and the substrate energy\({\mathcal {E}}_{\mathrm{sub}}\),

where\(\Gamma (t)\) is the dynamic film/vapor interface with\(\varvec{\nu }\) being the interface normal pointing into the vapor phase,\(\Gamma _{{\scriptscriptstyle {FS}}}\) and\(\Gamma _{{\scriptscriptstyle {VS}}}\) are the interfaces between film/substrate and vapor/substrate, respectively, and\(\gamma _{_{\scriptscriptstyle {FS}}}\) and\(\gamma _{_{\scriptscriptstyle {VS}}}\) are the corresponding surface energy densities. In order to model SSD by the diffuse-interface approach, we associate the vapor phase with\(\varphi \approx 1\) and the film phase with\(\varphi \approx -1\). Then, the Ginzburg–Landau type energy (1.5), up to a multiplicative constant, will approximate the sharp-interface energy\({\mathcal {E}}_\mathrm{inf}\). Moreover, the contribution to the wall energy\({\mathcal {E}}_{\mathrm{sub}}\) can be approximated by

where\(G(\varphi )\) is a smooth function satisfying\(G(\pm 1)=\pm \frac{1}{2}\), see Jiang et al. (2012), Dziwnik et al. (2017), Huang et al. (2019), Backofen et al. (2017) for SSD and (Jacqmin2000; Qian et al.2003) for moving contact lines in fluid mechanics.

There are several results on the existence of weak solutions for the degenerate Cahn–Hilliard equation (1.3) with homogeneous boundary conditions or its variants with inhomogeneous boundary conditions, see Elliott and Garcke (1996), Dai and Du (2016), Yin (1992). However, little is known about the anisotropic case except the work in Dziwnik (2019) which focuses on a particularn-fold anisotropy in two space dimensions.

The main aim of this work is to develop a diffuse-interface approach to SSD in the case of anisotropic surface energies based on the energy contributions (1.5) and (1.9). The obtained diffuse-interface model consists of a degenerate anisotropic Cahn–Hilliard equation with appropriate boundary conditions. We study the sharp-interface limit and show the existence of weak solutions to the diffuse-interface model.

The rest of the paper is organized as follows. In Sect. 2, we review a sharp-interface model for SSD and then introduce a diffuse-interface model based on a gradient flow approach. We then derive the sharp-interface limit from a regularized model with the help of asymptotic expansions in Sect. 3. In Sect. 4, we prove the existence of weak solutions to the diffuse-interface model. Numerical tests are presented in Sect. 5, where a comparison between sharp-interface approximations and diffuse-interface approximations is made.

In this section, we first review a sharp-interface model for SSD with anisotropic surface energies in two or three space dimensions. Then, we propose a suitable diffuse-interface model to approximate this sharp-interface model. Here, we note that there exist several works on the modeling of SSD using a diffuse-interface approach in the literature. However, these works consider either the isotropic case, e.g., Jiang et al. (2012), Backofen et al. (2017), or the anisotropic case in 2d, e.g., Dziwnik et al. (2017).

2.1The Sharp-Interface Model

We consider the dewetting of a solid thin film on a flat substrate in\({\mathbb {R}}^d\) with\(d\in \{2,3\}\), as shown in Fig. 1. We parameterize the interface of\(\Gamma (t)\) over the initial interface as follows

where\(T>0\) is a prescribed final time. The induced velocity is then given by

where\(\Gamma (0)\) is a smooth hypersurface with boundary. The sharp-interface model for SSD (cf. Cahn and Taylor1994; Taylor and Cahn1994; Barrett et al.2010a; Jiang et al.2020) reads as

which has to hold for all\(t\in [0, T]\) and all points on\(\Gamma (t)\). Here,\({\mathcal {V}}=\mathcal {\mathbf {V}}\cdot \varvec{\nu }\) is the normal velocity,\(\varvec{\nu }\) is the unit normal to\(\Gamma (t)\) pointing into the vapor, and\(\nabla _s\) is the surface gradient operator on\(\Gamma (t)\). Besides,\(D(\varvec{\nu })\) is an orientation-dependent mobility (cf.  Taylor and Cahn1994). The functionD needs to be defined for unit vectors, but here, we extend its domain to\({\mathbb {R}}^d\) such that it is positively homogeneous of degree one. The term\(\varkappa _\gamma \) represents the anisotropic mean curvature, and\(\gamma ^\prime (\varvec{\nu })\) is the Cahn–Hoffman vector, where\(\gamma ^\prime \) denotes the gradient of\(\gamma \) (cf.  Hoffman and Cahn1972). The above equations are subject to the following boundary conditions at the contact line, where the film/vapor interface\(\Gamma (t)\) meets the substrate:

• Attachment condition

  $$\begin{aligned} \mathcal {\mathbf {V}}\cdot \mathbf {n}_w = 0, \end{aligned}$$

  (2.2a)

• Contact angle condition

  $$\begin{aligned} \gamma ^\prime (\varvec{\nu })\cdot \mathbf {n}_w + \sigma =0, \end{aligned}$$

  (2.2b)

• Zero-flux condition

  $$\begin{aligned} D(\varvec{\nu })\,\nabla _s\varkappa _\gamma \cdot \mathbf {n}_c=0, \end{aligned}$$

  (2.2c)

where

denotes the difference of the substrate energy densities across the contact line. Here,\(\mathbf {n}_w\) is the unit normal to the substrate and points in the direction of the substrate, and\(\mathbf {n}_c\) is the conormal vector of\(\Gamma (t)\), i.e., it is the outward unit normal to\(\partial \Gamma (t)\) and it lies within the tangent plane of\(\Gamma (t)\). We observe that (2.2b) enforces an angle condition between the Cahn–Hoffman vector\(\gamma ^\prime (\varvec{\nu })\) and the substrate unit normal\(\mathbf {n}_w\). For example, in the isotropic case,\(\gamma (\mathbf {p}) = |\mathbf {p}|\), the Cahn–Hoffman vector reduces to the normal\(\varvec{\nu }\), and so if\(\sigma =0\), the condition (2.2b) encodes a\(90^\circ \) contact angle between the film/vapor interface and the substrate.

We assume that the anisotropy function\(\gamma \) that belongs to\(C^2\big ({\mathbb {R}}^d \setminus \{\mathbf {0}\}\big ) \cap C({\mathbb {R}}^d,{\mathbb {R}}_{\ge 0})\) is convex and satisfies\(\gamma >0\) on\({\mathbb {R}}^d\setminus \{{\mathbf {0}}\}\). We further assume that\(\gamma \) is positively homogeneous of degree one, meaning that

This immediately implies\(\gamma ({\mathbf {0}}) = 0\) and the gradient of\(\gamma (\mathbf {p})\) satisfies

Similarly, the orientation-dependent mobility function\(D\in C^2\big ({\mathbb {R}}^d \setminus \{\mathbf {0}\}\big ) \cap C({\mathbb {R}}^d,{\mathbb {R}}_{\ge 0})\) is assumed to satisfy\(D>0\) on\({\mathbb {R}}^d\setminus \{\mathbf {0}\}\) and

Consequently, for the map

introduced in (1.5), we have\(A\in C^2\big ({\mathbb {R}}^d \setminus \{\mathbf {0}\}\big ) \cap C({\mathbb {R}}^d,{\mathbb {R}}_{\ge 0})\). It also follows directly from (2.4) that the relations

hold for all\(\mathbf {p}\in {\mathbb {R}}^d\setminus \{\mathbf {0}\}\) and all\(\lambda >0\). Here,\(A^{\prime }\) and\(A^{\prime \prime }\) denote the gradient and the Hessian ofA, respectively.

2.2The Diffuse-Interface Model

Geometric setup for SSD in a bounded domain\(\Omega \) with\(\Omega =\overline{\Omega _-}(t)\cup \overline{\Omega _+}(t)\), where\(\Omega _{-}(t):=\{\mathbf {x}\in \Omega : \varphi (\mathbf {x},t)<0\}\) and\(\Omega _{+}(t) := \{\mathbf {x}\in \Omega : \varphi (\mathbf {x},t)>0\}\)

Full size image

Let\(\varphi :\Omega \times [0, T]\rightarrow {\mathbb {R}}\) be an order parameter such that the zero level set\(\{\mathbf {x}\in \Omega : \varphi (\mathbf {x}, t) =0\}\) approximates the film/vapor interface\(\Gamma (t)\),\(\{\mathbf {x}\in \Omega : \varphi (\mathbf {x}, t) < 0\}\) corresponds to the region occupied by the thin film at timet, whereas\(\{\mathbf {x}\in \Omega : \varphi (\mathbf {x}, t) > 0\}\) represents the region occupied by the vapor at timet (see Fig. 2). In addition,\(\Gamma _w\subset \partial \Omega \) models the boundary of the substrate. As a combination of (1.5) and (1.9), the total free energy of the system is given by

where\({c_{_{F}}}=\int _{-1}^1\sqrt{2F(s)}\,\mathrm{d}s\) and\(|\Gamma _w|= \int _{\Gamma _w}\, \mathrm {d}S\). This choice of\({c_{_{F}}}\) ensures that

for sufficiently small\({\varepsilon }>0\). Besides, the constant term\(-\gamma _{_{\scriptscriptstyle {VS}}}|\Gamma _w|\) was added to the total energy such that\({\mathcal {E}}(\varphi )\) now only depends on the single parameter\(\sigma \) (see (2.3)) instead of on\(\gamma _{_{\scriptscriptstyle {VS}}}\) and\(\gamma _{_{\scriptscriptstyle {FS}}}\). We next derive the diffuse-interface model. To this end, we use the smooth double-well potential

This implies

We further choose

which yields\(G(\pm 1) = \pm \frac{1}{2}\) and\(G^\prime (\pm 1) = 0\). Let\(\psi :\Omega \rightarrow {\mathbb {R}}\) be a sufficiently smooth function. Then, the first variation in the total free energy (2.7) in the direction of\(\psi \) can be computed as

where\(\mathbf {n}\) is the outward unit normal to\(\partial \Omega \setminus \Gamma _w\) and\(\mathbf {n}_w\) is the outward unit normal to\(\Gamma _w\), as defined previously. The following diffuse-interface model for SSD can be interpreted as a weighted\(H^{-1}\)-gradient flow of the energy functional (2.7):

Here,\(\alpha >0\) is a time scaling coefficient,\(m(\varphi )\) is the degenerate mobility given by

and\(\beta (\nabla \varphi )\) is defined as

and so is positively homogeneous of degree zero.

We now write\(\Sigma = \partial \Omega \times (0,T]\) and\(\Sigma _w = \Gamma _w\times (0,T]\). On\(\Sigma _w\), we impose the boundary conditions

Here, the first equation is the zero-flux condition on the boundary, whereas the second equation guarantees the integral over\(\Gamma _w\) in (2.10) vanishes. Moreover, on\(\Sigma \setminus \Sigma _w\), we impose the natural boundary conditions

Remark 2.1

It is also possible to consider the double-obstacle potential (1.2c) along with the mobility\(m(\varphi ) = (1-\varphi ^2)_+\). Then, the corresponding diffuse-interface model consists of (2.11a) and the variational inclusion

instead of (2.11b).

We consider the smooth double-well potential introduced in (2.8) and regularize the coefficients\(m(\varphi )\) and\(\beta (\nabla \varphi )\) of the diffuse-interface model (2.11) with the help of the interfacial parameter\({\varepsilon }\) by defining

where\(r\ge 2\). The regularized diffuse-interface model is then given by

We note that the introduction of the three regularization terms\({\varepsilon }^r\) in (3.1) and (3.2) allows for a mathematical analysis of (3.3) in Sect. 4 below. In fact, on defining

we have

Moreover, by choosing\(r\ge 2\), we ensure that the sharp-interface limit of (3.3) is unchanged compared to the limit of (2.11).

We now formally derive the sharp-interface limit of the regularized model (2.11) via the method of matched asymptotic expansions. We suppose that for\({\varepsilon }>0\),\((\varphi ^{\varepsilon }, \mu ^{\varepsilon })\) is the solution of the regularized diffuse-interface model (3.3). Then, we write

to denote the interface and the contact line, respectively. We further assume that their limits as\(\varepsilon \rightarrow 0\) are given by\(\Gamma (t)\) and\(\Lambda (t)\), respectively. We introduce a local parameterization for\(\Gamma (t)\) on an open subset\({\mathcal {O}}\subset {\mathbb {R}}^{d-1}\) by

Our asymptotic analysis for the interface dynamics will follow similar techniques in the literature for degenerate Cahn-Hilliard equations, see, e.g., Cahn et al. (1996), Dai and Du (2014), for the isotropic case and Rätz et al. (2006), Dziwnik et al. (2017) for the anisotropic case in 2d.

3.1Outer Expansions

Away from the interface and the contact line, we assume that the following ansatz holds

Moreover, in view of (3.1) and (3.2), we know that

since\(r\ge 2\) and\(\beta ^{\varepsilon }(\mathbf {p}) = \beta (\mathbf {p}) + O({\varepsilon }^r)\), where\(\beta ^\prime \) denotes the gradient of\(\beta \). Plugging the expansions (3.6) and (3.7) into (3.3a) and (3.3b) gives

As the energy (1.5) is expected to be bounded at leading order, it needs to hold\(F(\varphi _0) = 0\). This means that\(\varphi _0\) attains only the values\(-1\) and 1. Hence,\(\mu _{-1}=0\). We now define

as the outer regions, meaning that\(\varphi _0 = \pm 1\) in\(\Omega _\pm (t)\).

3.2Inner Expansions

In the inner region near the interface\(\Gamma (t)\), we introduce the annular neighborhood

where\(\mathrm{d}(\mathbf {x}, t)\) represents the signed distance of\(\mathbf {x} \) to\(\Gamma (t)\), defined to be positive in\(\Omega _+(t)\). Assuming\(\Gamma (t)\) to be sufficiently smooth, we find a\(\delta >0\) such that for every\(\mathbf {x}\in {\mathcal {N}}(t)\), there exist unique vectors\(\mathbf {r}(\mathbf {x}, t)\) and\(\mathbf {s}(\mathbf {x}, t)\) such that

Here,\({\varvec{\nu }}(\mathbf {s} ,t)\) is the unit normal vector on\(\Gamma (t)\) at\(\mathbf {r}(\mathbf {s}, t)\) pointing into\(\Omega _+(t)\).

Due to rapid changes of\(\varphi ^{\varepsilon }\) in normal direction, we introduce the stretched variable\(\rho (\mathbf {x}, t) = {\varepsilon }^{-1} \mathrm{d}(\mathbf {x}, t)\). Any scalar function\(b(\mathbf {x}, t)\) can be expressed in the new coordinate system as\(b(\mathbf {x}, t) = {\overline{b}}(\mathbf {s}(\mathbf {x}, t),\rho (\mathbf {x}, t),t)\). For any vector field\(\mathbf {b}(\mathbf {x}, t)\), we use an analogous notation. Without loss of generality, we assume that\(\{\mathbf {t}_1,~\mathbf {t}_2,\cdots , \mathbf {t}_{d-1}\}\) forms an orthonormal basis of the tangent space of\(\Gamma (t)\) at the point\(\mathbf {r}(\mathbf {s}, t)\) such that

where\(\kappa _j\) is the principal curvature of\(\Gamma (t)\) at the point\(\mathbf {r}(\mathbf {s}, t)\) in the direction of\(\mathbf {t}_j\). As in Dai and Du (2014), we obtain the identities

Therefore, using the new coordinates, we calculate

where\(\nabla _s=\sum _{j=1}^{d-1}\mathbf {t}_j\,\partial _{s_j}\) denotes the surface gradient operator on\(\Gamma (t)\),

and\({\mathcal {V}}\) is the velocity of\(\Gamma (t)\) in the direction of\({\varvec{\nu }}\), i.e.,\({\mathcal {V}} = -\partial _t\mathrm{d} = - {\varepsilon }\,\partial _t\rho \).

In the inner region, we assume the following expansions

In particular, on assuming\(\partial _\rho \Phi _0>0\), we have, similarly to (3.7), that

where we have used the fact that\(\beta \) is positively homogeneous of order zero.

Plugging (3.10) and (3.11) into (3.3a), we obtain the leading order term

which implies that\(m(\Phi _0)\partial _\rho M_{-1}\) is independent of\(\rho \), i.e., it can be expressed as

In addition, using the matching condition

we infer\(J(\mathbf {s},t) = 0\) due to the degenerate mobility\(m(\Phi _0)\). Since\(m(s)>0\) if\(s\in (-1,1)\), we thus conclude that\(M_{-1}\) is independent of\(\rho \). By the matching condition\(\lim _{\rho \rightarrow \pm \infty }M_{-1}(\mathbf {s},t) = \mu _{-1}\), we obtain

For the terms of order\(O(\frac{1}{{\varepsilon }^3})\), we obtain

Repeating the above line of argument, we deduce

Using the fact that\(M_{-1}=0\) and\(\partial _\rho M_0 = 0\), we then have the following expansions

Considering the order\(O(\frac{1}{{\varepsilon }^2})\) of (3.3a), we obtain that

Similarly, by using the matching conditions, we arrive at

At\(O(\frac{1}{{\varepsilon }})\), using\(\partial _\rho M_1=0\) and (3.15), we have

We next consider the expansion of (3.3b). Using the identities in (2.6) and assuming\(\partial _\rho \Phi _0>0\), we expand the anisotropic term\(A^\prime (\nabla \Phi ^{\varepsilon })\) as follows:

This then yields

Now, plugging (3.10) into (3.3b), we obtain for the leading order term that

Using the translation identity\(\Phi _0(0) = 0\), we then obtain

Similarly, theO(1) term resulting from (3.3b) implies

Multiplying (3.21) by\(\partial _\rho \Phi _0\) and then integrating from\(-\infty \) to\(\infty \) with respect to\(\rho \) yields

Differentiating (3.19) with respect to\(\rho \) gives

Therefore, since\(\lim _{\rho \rightarrow \pm \infty }\partial _\rho \Phi _0=0\) and\(\lim _{\rho \rightarrow \pm \infty }\Phi _1=0\), we compute

via integration by parts. Then, using (3.20) and the matching condition in (3.13), we can reformulate (3.22) as

It further follows from (3.20) that\(\partial _\rho (\nabla _s\Phi _0)= \nabla _s(\partial _\rho \Phi _0)\). We thus have

which yields

where\(\varkappa _\gamma = - \nabla _s\cdot \gamma ^\prime ({\varvec{\nu }})\) is the weighted mean curvature defined in (2.1b).

We now return to (3.18) and integrate it with respect to\(\rho \) from\(-\infty \) to\(+\infty \). Using the fact that\(\lim _{\rho \rightarrow \pm \infty }m(\Phi _0)\partial _\rho M_2= 0\), we get

where we recall (3.20) and also use the identities

We thus obtain

3.3Expansions Near the Intersection with the Substrate

We next study the expansions near the intersection with the substrate using the technique discussed in Dziwnik et al. (2017), Owen et al. (1990).

3.3.1The Boundary Layer Near the Wall

In the boundary layer near\(\Gamma _w\), we first introduce the variable\(\eta = {\varepsilon }^{-1}\,\mathrm{d}_w(\mathbf {x})\), where\(\mathrm{d}_w(\mathbf {x})\) represents the distance from\(\mathbf {x}\) to the wall\(\Gamma _w\). Then, for a scalar function\(b(\mathbf {x}, t)\), we can write it as\(b(\mathbf {x}, t) = {\widehat{b}}(\eta , \mathbf {y}, t)\), where\(\mathbf {y}\) is the\((d-1)\)-dimensional coordinate system that is orthogonal to\(\eta \). This implies

We consider the expansions

and plug them into (3.3a) and (3.3b). The leading order terms yield

where\({\widehat{\beta }}_0 =\beta (-\partial _\eta {\widehat{\varphi }}_0\,\mathbf {n}_w)\). At the boundary\(\eta =0\), it holds

Thus, from (3.28a) and (3.29b), we obtain

Multiplying (3.28b) by\(\partial _\eta {\widehat{\varphi }}_0\) and using the identities in (2.6), we arrive at

Integrating (3.30) over\(\eta \) leads to

where\(c(\mathbf {y}, t)=0 \) due to the matching condition when\(\eta \rightarrow \infty \). This implies

3.3.2The Inner Layer Near the Contact Line

We assume that a local parameterization of the contact line\(\Lambda (t)\) is given by

where in the case\(d=2\), we simply set\({\mathcal {O}}_w = \{0\}\). For a contact point\(\mathbf {x}_c\in \Lambda (t)\), we then introduce an interior layer near it. Precisely, for any\(\mathbf {x}\) in the plane that contains\(\mathbf {x}_c\) and is spanned by\(\mathbf {n}_s\) and\(\mathbf {n}_w\), we write

where\(\mathbf {n}_s\) is the unit normal to\(\Lambda (t)\) on the wall\(\Gamma _w\) and pointing into\(\Omega _+(t)\). For a scalar function\(b(\mathbf {x}, t)\), we can rewrite it as\(b(\mathbf {x}, t) = {\widetilde{b}}(s_w, \xi , \eta , t)\). In a similar manner to (3.9), we compute

where\(\nabla _t^\Lambda {\widetilde{b}} = \partial _t{\widetilde{b}} +\partial _ts_w\cdot \nabla _{s_w}{\widetilde{b}}\). We then consider the expansions

and plug them into (3.3a) and (3.3b). By defining\(\nabla _c = \mathbf {n}_s\,\partial _\xi -\mathbf {n}_w\,\partial _\eta \), the leading order term yields

where\({\widetilde{\beta }}_0=\beta (\nabla _c{\widetilde{\varphi }}_0)\). Similarly, the leading order terms of the boundary conditions (3.3c) and (3.3d) give

Besides, we have the matching condition

Now, multiplying (3.39b) by\(\partial _\xi {\widetilde{\varphi }}_0\) and integrating the resulting equation in a box\(R:=[-\xi _1,\xi _1]\times [0,\eta _1]\), we get

which can be rewritten as

by using the identity

For the first integral in (3.42), applying Gauss’s theorem and using the matching condition in (3.41) as well as the fact\(\lim _{\xi \rightarrow +\infty }\partial _{\xi }{\widetilde{\varphi }}_0=0\), we have

where we have used (3.31) and (3.32).

We then apply Gauss’s theorem to the second integral in (3.42). Recalling the boundary condition (3.40b), we obtain

Sending\(\xi _1\rightarrow +\infty \) and recalling (2.9), we obtain

Sketch of the local coordinates\((\xi ,\eta )\) and\((\rho ,\zeta )\) at a contact point\(\mathbf {x}_c\), where\(\theta _d\in (0,\pi )\) is the contact angle

Full size image

Next, we rewrite the termI in terms of the new coordinate system\((\rho ,\zeta )\), which can be regarded as a transformation from\((\eta ,\xi )\) with a counterclockwise rotation of\(\theta _d\) in the plane (see Fig. 3). Precisely, it holds that

and thus,

Moreover, we have

where\(\mathbf {n}_c\) is the conormal vector of\(\Gamma (t)\) at\(\mathbf {x}_c\). By (3.46), we can recast the termI as

By the matching condition\(\lim _{\zeta \rightarrow +\infty }{\widetilde{\varphi }}_0 = \Phi _0(\rho )\), we have\(\lim _{\zeta \rightarrow +\infty }\partial _{\zeta }{\widetilde{\varphi }}_0=0\). Then, it follows directly that

Collecting the results in (3.43), (3.45) and (3.48) yields that

which is exactly the anisotropic Young’s law in (2.2b).

We next derive the zero-flux condition. Similarly to the above, we integrate (3.39a) over the boxR. Applying Gauss’s theorem and using the boundary condition (3.40a) gives rise to

Taking\(\xi _1\rightarrow \infty \) and using fact\(\lim _{\xi \rightarrow \pm \infty }{\widetilde{\varphi }}_0 = \pm 1\), as well as\(m({\widetilde{\varphi }}_0)=0\), we get\(III=0\). On recalling (3.46), as well as the matching conditions,

we get in the case of\(\xi _1, \eta _1\rightarrow \infty \) that

This yields the zero-flux condition

In addition, the attachment condition in (2.2a) follows naturally.

In summary, we thus obtain the following system of equations as the sharp-interface limit of the regularized diffuse-interface model (3.3):

Remark 3.1

In the case of the double-obstacle potential (1.2c) and the degenerate mobility\(m(\varphi )=(1-\varphi ^2)_+\), we could obtain (3.52a) in a similar manner. But the leading order inner solution (3.20) should be replaced by

This yields\({c_{_{F}}}=\frac{\pi }{2}\) in (3.52a). The boundary conditions in (3.52b) can be derived similarly. It is also possible to consider the logarithmic potential (see (1.2b)) along with the mobility\(m(\varphi )=(1-\varphi ^2)_+\). If\(\theta = O({\varepsilon }^\xi )\) for some\(\xi >0\), it can be shown by means of the techniques from Cahn et al. (1996) that the same desired sharp-interface limit is obtained.

In this section, we analyze a general class of diffuse-interface models of the type

where\(\alpha ,{\varepsilon },{c_{_{F}}}\in {\mathbb {R}}_{>0}\) and\(\sigma \in {\mathbb {R}}\) are given constants. In contrast to the previous sections, the potential\(F:{\mathbb {R}}\rightarrow {\mathbb {R}}\) and\(G:{\mathbb {R}}\rightarrow {\mathbb {R}}\),\(A:{\mathbb {R}}^d\rightarrow {\mathbb {R}}\) and\(M:{\mathbb {R}}^d \times {\mathbb {R}}\rightarrow {\mathbb {R}}\) are general functions satisfying certain conditions that will be specified in Sect. 4.1. IfA,F,G,\(m^{\varepsilon }\) and\(\beta ^{\varepsilon }\) are chosen as in (2.5), (2.8), (2.9), (3.1) and (3.2), respectively, and ifM is defined by\(M(\mathbf {p},s):= \beta ^{\varepsilon }(\mathbf {p}) m^{\varepsilon }(s)\) for all\(\mathbf {p}\in {\mathbb {R}}^d\) and\(s\in {\mathbb {R}}\), then the system (4.1) is exactly the model (3.3) that was introduced in Sect. 3. The total free energy functional\({\mathcal {E}}:H^1(\Omega ) \rightarrow {\mathbb {R}}\) associated with the system (4.1), up to an additive constant, reads as

It is also possible to consider the system (4.1) forF being the double-obstacle potential, which can be expressed as

where the function

represents its regular part, and

denotes the indicator functional of the interval\([-1,1]\). In this case, (4.1b) needs to be represented by a variational inequality, see (4.18b).

4.1Notation and Preliminaries

Notation. In this section, we use the following notation: For any\(1 \le p \le \infty \) and\(k \ge 0\), the standard Lebesgue and Sobolev spaces on\(\Omega \) are denoted by\(L^p(\Omega )\) and\(W^{k,p}(\Omega )\). Their standard norms are written as\(\left\| \,\cdot \, \right\| _{L^p(\Omega )}\) and\(\left\| \,\cdot \, \right\| _{W^{k,p}(\Omega )}\). In the case\(p = 2\), these spaces are Hilbert spaces, and we write\(H^k(\Omega ) = W^{k,2}(\Omega )\). Here, we identify\(H^0(\Omega )\) with\(L^2(\Omega )\). For the Lebesgue and Sobolev spaces on\(\partial \Omega \), we use an analogous notation. For any Banach spaceX, its dual space is denoted by\(X'\), and the associated duality pairing by\( \left\langle \cdot \, , \, \cdot \right\rangle _X\). IfX is a Hilbert space, we write\((\cdot , \cdot )_X\) to denote its inner product. We further define

as the generalized spatial mean off, where\(\left| \Omega \right| \) denotes thed-dimensional Lebesgue measure of\(\Omega \). With the usual identification\(L^1(\Omega ) \subset H^1(\Omega )'\), it holds that\(\left\langle f \right\rangle _\Omega = \frac{1}{\left| \Omega \right| } \int _\Omega f \, \mathrm {d}x\) if\(f \in L^1(\Omega )\). In addition, we introduce

We point out that for every\(m\in {\mathbb {R}}\),\(H^1_{(m)}(\Omega )\) is an affine subspace of the Hilbert space\(H^1(\Omega )\). In the case\(m=0\), it is even a closed linear subspace, meaning that\(H^1_{(0)}(\Omega )\) is also a Hilbert space.

General assumptions We make the following general assumptions that are supposed to hold throughout this section.

1. A1

   The set\(\Omega \subset {\mathbb {R}}^d\) with\(d \in \{2,3\}\) is a bounded Lipschitz domain. Moreover,\(T>0\) denotes an arbitrary final time.

2. A2

   The function\(G:{\mathbb {R}}\rightarrow {\mathbb {R}}\) is non-negative and twice continuously differentiable. Moreover, there exists an exponent\(q\in [2,4)\), as well as positive constants\(C_{G}\) and\(C_{G'}\), such that

   $$\begin{aligned} G(s) \le C_{G}(1+\left| s \right| ^q) \quad \text {and}\quad \left| G'(s) \right| \le C_{G'}(1+\left| s \right| ^{q-1}) \end{aligned}$$

   for all\(s\in {\mathbb {R}}\).

3. A3

   “The function\(A:{\mathbb {R}}^d\rightarrow {\mathbb {R}}\) is continuously differentiable, positively homogeneous of degree 2 and positive on\({\mathbb {R}}^{d} \ {0}\). Hence there exist constants\(A_0,A_1\in {\mathbb {R}}\) with\(0<A_0\le A_1\) such that”

   $$\begin{aligned} A_0 \left| \mathbf {p} \right| ^2 \le A(\mathbf {p}) \le A_1 \left| \mathbf {p} \right| ^2 \quad \text {for all}\,\, \mathbf {p}\in {\mathbb {R}}^d. \end{aligned}$$

   The gradient\(A':{\mathbb {R}}^d\rightarrow {\mathbb {R}}^d\) is strongly monotone, i.e., there exists a constant\(a_0>0\) such that

   $$\begin{aligned} \big (A'(\mathbf {p})-A'(\mathbf {q})\big )\cdot (\mathbf {p}-\mathbf {q}) \ge a_0 \left| \mathbf {p}-\mathbf {q} \right| ^2 \quad \text {for all}\,\, \mathbf {p},\mathbf {q}\in {\mathbb {R}}^d, \end{aligned}$$

   which implies thatA is strongly convex and thus strictly convex. Moreover, there exists a constant\(a_1>0\) such that

   $$\begin{aligned} \left| A'(\mathbf {p}) \right| \le a_1 \left| \mathbf {p} \right| \quad \text {for all}\,\, \mathbf {p}\in {\mathbb {R}}^d. \end{aligned}$$
4. A4

   The function\(M:{\mathbb {R}}^d\times {\mathbb {R}}\rightarrow {\mathbb {R}}\) is continuous and there exist constants\(M_0,M_1\in {\mathbb {R}}\) with\(0<M_0\le M_1\) such that

   $$\begin{aligned} M_0 \le M(\mathbf {p},s) \le M_1 \quad \text {for all}\,\, \mathbf {p}\in {\mathbb {R}}^d \,\,\text {and}\,\, s\in {\mathbb {R}}. \end{aligned}$$

Remark 4.1

1. (a)

   We point out that the choices

   $$\begin{aligned} G(s)&:= \frac{1}{4}(3s - s^3)&\quad \quad \text {for all}\,\, s\in {\mathbb {R}}&\quad \quad \big (\text {cf.}\ 2.9\big ),\\ M(\mathbf {p},s)&:= \beta ^{\varepsilon }(\mathbf {p})\, m^{\varepsilon }(s) \end{aligned}$$

   with

   $$\begin{aligned} m^{\varepsilon }(s)&:= {\varepsilon }^r + (1-s^2)_+^2&\quad \quad \text {for all}\,\, s\in {\mathbb {R}}&\quad \quad \big (\text {cf.}\ 3.1\big ),\\ \beta ^{\varepsilon }(\mathbf {p})&:= \sqrt{ \frac{d_1^2\,{\varepsilon }^{2r} + D^2({\mathbf {p}})}{\gamma _0^2{\varepsilon }^{2r} + \gamma ^2(\mathbf {p})} }&\quad \quad \text {for all}\,\, \mathbf {p}\in {\mathbb {R}}^d&\quad \quad \big (\text {cf.}\ 3.2\big ), \end{aligned}$$

   are admissible as they satisfy the conditions imposed in A2 (with\(q=3\)) and A4.

2. (b)

   Suppose that the function\(\gamma \) that was introduced in Sect. 2.1 additionally satisfies the following convexity condition: There exists a constant\(\alpha _0>0\) such that

   $$\begin{aligned} \gamma ^{\prime \prime }(\mathbf {p})\mathbf {q}\cdot \mathbf {q}\ge \alpha _0|\mathbf {q}|^2 \qquad \text{ for } \text{ all } \quad \mathbf {p},\mathbf {q}\in {\mathbb {R}}^d\;\text{ with }\quad |\mathbf {p}|=1\quad \text{ and }\quad \mathbf {p}\cdot \mathbf {q} = 0, \end{aligned}$$

   (4.6)

   where\(\gamma ^{\prime \prime }\) represents the Hessian of\(\gamma \). Thus, the function

   $$\begin{aligned} A:{\mathbb {R}}^d\rightarrow {\mathbb {R}},\quad A(\mathbf {p}) = \tfrac{1}{2} \gamma ^2(\mathbf {p}) \end{aligned}$$

   is admissible as it satisfies all conditions imposed in assumption A3. In particular, as shown in Gräser et al. (2013), the convexity condition (4.6) ensures that\(A'\) is strongly monotone.

A special inner product on\(H^{-1}_{(0)}(\Omega )\). We now introduce a certain inner product on the function space\(H^{-1}_{(0)}(\Omega )\) based on the solution operator of a suitable elliptic problem. Therefore, let\(a\in L^\infty (\Omega )\) be a uniformly positive function, i.e., there exist\(a_0, a_1 \in {\mathbb {R}}\) with\(0<a_0<a_1\) such that

Then, for every\(f\in H^{-1}_{(0)}(\Omega )\), there exists a unique weak solution\(u_f\in H^{1}_{(0)}(\Omega )\) of the elliptic boundary value problem

meaning that

We can thus define a solution operator

We next define the bilinear form

which defines an inner product on\(H^{-1}_{(0)}(\Omega )\) sincea is uniformly positive and\(\nabla S_a(f) = 0\) a.e. in\(\Omega \) already implies\(f=0\) a.e. in\(\Omega \). Its induced norm is given by

We point out that on the space\(H^{-1}_{(0)}(\Omega )\), the norm\(\left\| \,\cdot \, \right\| _{S_a}\) is equivalent to the standard operator norm\(\left\| \,\cdot \, \right\| _{(H^1(\Omega ))'}\,\). The bilinear form\( \left( \cdot \,{,}\, \cdot \right) _{S_a}\) also defines an inner product on the space\(H^{1}_{(0)}(\Omega )\). Moreover,\(\left\| \,\cdot \, \right\| _{S_a}\) is also a norm on\(H^{1}_{(0)}(\Omega )\) but the space is not complete with respect to this norm.

4.2Existence of Weak Solutions

For ease of presentation, in what follows, we simply fix\(\alpha ={\varepsilon }=\sigma ={c_{_{F}}}=1\), since the precise choice of these values has no impact on the mathematical analysis.

4.2.1Weak Solutions for Smooth Potentials

In this subsection, we make the following assumption on the potentialF:

1. F1

   The potential\(F:{\mathbb {R}}\rightarrow {\mathbb {R}}\) is continuously differentiable. Moreover, there exists an exponent\(p\in [2,6)\), as well as non-negative constants\(B_F\),\(C_{F}\) and\(C_{F'}\) such that

   $$\begin{aligned} - B_{F} \le F(s) \le C_{F}(1+\left| s \right| ^p) \quad \text {and}\,\,\quad \left| F'(s) \right| \le C_{F'}(1+\left| s \right| ^{p-1}). \end{aligned}$$

   for all\(s\in {\mathbb {R}}\).

Obviously, the smooth double-well potential introduced in (1.2a) fulfills F1 with\(p=4\). However, the logarithmic potential (see (1.2b)) and the double-obstacle potential (see (1.2c)) do not satisfy this assumption.

A weak solution of the general diffuse-interface model (4.1) is then defined as follows.

Definition 4.2

Suppose that the assumptions A1–A4 and F1 are fulfilled, and let\({\varphi _0\in H^1(\Omega )}\) be any initial datum. Then, the pair\((\varphi ,\mu )\) is called a weak solution to system (4.1) if the following properties hold:

1. (i)

   The functions\(\varphi \) and\(\mu \) have the following regularity:

   $$\begin{aligned} \varphi&\in C^{0,1/4}\big ([0,T];L^2(\Omega )\big ) \cap L^\infty \big (0,T;H^1(\Omega )\big ) \cap H^1\big (0,T;H^1(\Omega )'\big ),\\ \mu&\in L^2\big (0,T;H^1(\Omega )\big ). \end{aligned}$$
2. (ii)

   The pair\((\varphi ,\mu )\) satisfies the weak formulations

   $$\begin{aligned} \big < \partial _t\varphi \, , \, \zeta \big >_{H^1(\Omega )}&= - \int _{\Omega }M(\nabla \varphi ,\varphi )\, \nabla \mu \cdot \nabla \zeta \, \mathrm {d}x,\end{aligned}$$

   (4.12a)
   $$\begin{aligned} \int _{\Omega }\mu \, \eta \, \mathrm {d}x&= \int _{\Omega }A'(\nabla \varphi )\cdot \nabla \eta + F'(\varphi )\,\eta \, \mathrm {d}x+ \int _{\Gamma _w}G'(\varphi ) \, \eta \, \mathrm {d}S\end{aligned}$$

   (4.12b)

   a.e. on [0, T] for all test functions\(\zeta ,\eta \in H^1(\Omega )\). Moreover,\(\varphi \) satisfies the initial condition

   $$\begin{aligned} \varphi (0) = \varphi _0 \qquad \text {a.e. in}\,\, \Omega . \end{aligned}$$

   (4.13)
3. (iii)

   The pair\((\varphi ,\mu )\) satisfies the weak energy dissipation law

   $$\begin{aligned} {\mathcal {E}}\big (\varphi (t)\big ) + \frac{1}{2} \int _0^t \int _{\Omega }M(\nabla \varphi ,\varphi )\, \left| \nabla \mu \right| ^2 \, \mathrm {d}x\, \mathrm {d}t\le {\mathcal {E}}(\varphi _0) \qquad \text {for almost all}\,\, t\in [0,T]. \end{aligned}$$

   (4.14)

The existence of such a weak solution is ensured by the following theorem.

Theorem 4.3

Suppose that the assumptions A1–A4 and F1 are fulfilled, and let\({\varphi _0\in H^1(\Omega )}\) be any initial datum. Then, there exists a weak solution\((\varphi ,\mu )\) to the system (4.1) in the sense of Definition 4.2.

The proof of this theorem is presented in Sect. 4.3.

In the next subsection, we intend to prove the existence of a weak solution to the diffuse-interface model (4.1) for the double-obstacle potential (1.2c). Our strategy is to approximate the double-obstacle potential by a sequence of regular potentials. To this end, in Corollary 4.4, we will present an additional uniform estimate for\(F'(\varphi )\), where\((\varphi ,\mu )\) is a weak solution to (4.1) with a regular potentialF satisfying the following assumption:

1. F2

   The potential\(F:{\mathbb {R}}\rightarrow {\mathbb {R}}\) is twice continuously differentiable and there exist constants\(c_0,c_1 \ge 0\) such that

   $$\begin{aligned} - c_0 \le F''(s) \le c_1 \quad \text {for all}\,\, s\in {\mathbb {R}}. \end{aligned}$$

   (4.15)

We point out that if F2 is fulfilled, then F1 holds with\(p=2\).

Corollary 4.4

Suppose that the assumptions A1–A4 and F2 are fulfilled. Let\(\varphi _0\in H^1(\Omega )\) be any initial datum satisfying\(\left| \left\langle \varphi _0 \right\rangle _\Omega \right| \le 1 - \kappa \) for some\(\kappa \in (0,1]\), and let\((\varphi ,\mu )\) be a corresponding weak solution. Then, there exists a constant\(c>0\) depending only on\(\varphi _0\),\({\mathcal {E}}(\varphi _0)\),\(c_0\) and the constants in A1–A4, but not on\(c_1\), such that

where\( R:= \left| \left\langle \varphi _0 \right\rangle _\Omega \right| + \frac{\kappa }{2} < 1\).

Remark 4.5

In Corollary 4.4, the assumption\(\left| \left\langle \varphi _0 \right\rangle _\Omega \right| \le 1 - \kappa \) is made in order to ensure\(R\le 1\), which is crucial for later use. However, without this assumption, a similar estimate can be derived if\(R>1\) is allowed. For instance, choosing\(R:=\left| \left\langle \varphi _0 \right\rangle _\Omega \right| +1\), we obtain the estimate

instead of (4.16) even without the mean value assumption.

4.2.2Weak Solutions for the Double-Obstacle Potential

In this subsection, we assume that\(F = F_0 + I_{[-1,1]}\) is the double-obstacle potential as introduced in (4.3). Then, a weak solution of the general diffuse-interface model (4.1) is defined as follows.

Definition 4.6

Suppose that the assumptions A1–A4 are fulfilled, and let\(\varphi _0\in H^1(\Omega )\) be any initial datum satisfying\(\left| \varphi _0 \right| \le 1\) a.e. in\(\Omega \). Then, the pair\((\varphi ,\mu )\) is called a weak solution to system (4.1) if the following properties hold:

1. (i)

   The functions\(\varphi \) and\(\mu \) have the following regularity:

   $$\begin{aligned} \varphi&\in C^{0,1/4}\big ([0,T];L^2(\Omega )\big ) \cap L^\infty \big (0,T;H^1(\Omega )\big ) \cap H^1\big (0,T;H^1(\Omega )'\big ),\\ \mu&\in L^2\big (0,T;H^1(\Omega )\big ). \end{aligned}$$
2. (ii)

   It holds that\(\left| \varphi \right| \le 1\) a.e. inQ and the pair\((\varphi ,\mu )\) satisfies the weak formulation

   $$\begin{aligned} \big < \partial _t\varphi \, , \, \zeta \big >_{H^1(\Omega )}&= - \int _{\Omega }M(\nabla \varphi ,\varphi )\, \nabla \mu \cdot \nabla \zeta \, \mathrm {d}x, \end{aligned}$$

   (4.18a)

   for all\(\zeta \in H^1(\Omega )\) as well as the variational inequality

   $$\begin{aligned} \iint _Q \mu \, (\varphi -\eta ) \, \mathrm {d}x\, \mathrm {d}t&\ge \iint _Q A'(\nabla \varphi )\cdot (\nabla \varphi -\nabla \eta ) +F_0'(\varphi )(\varphi -\eta ) \, \mathrm {d}x\, \mathrm {d}t\nonumber \\&\qquad + \iint _{\Sigma _w}G'(\varphi ) \, (\varphi -\eta ) \, \mathrm {d}S\, \mathrm {d}t\end{aligned}$$

   (4.18b)

   for all\(\eta \in L^2(0,T;H^1(\Omega ))\) with\(\left| \eta \right| \le 1\) a.e. inQ. Moreover,\(\varphi \) satisfies the initial condition

   $$\begin{aligned} \varphi (0) = \varphi _0 \qquad \text {a.e. in}\,\, \Omega . \end{aligned}$$

   (4.19)
3. (iii)

   The pair\((\varphi ,\mu )\) satisfies the weak energy dissipation law

   $$\begin{aligned} {\mathcal {E}}\big (\varphi (t)\big ) + \frac{1}{2} \int _0^t \int _{\Omega }M(\nabla \varphi ,\varphi )\, \left| \nabla \mu \right| ^2 \, \mathrm {d}x\, \mathrm {d}t\le {\mathcal {E}}(\varphi _0) \qquad \text {for almost all}\,\, t\in [0,T]. \end{aligned}$$

   (4.20)

The existence of such a weak solution is ensured by the following theorem.

Theorem 4.7

Suppose that the assumptions A1–A4 are fulfilled, and let\(\varphi _0\in H^1(\Omega )\) be any initial datum satisfying\(\left| \varphi _0 \right| \le 1\) a.e. in\(\Omega \) and\(\left| \left\langle \varphi _0 \right\rangle _\Omega \right| \le 1-\kappa \) for some\(\kappa \in (0,1]\). Then, there exists a weak solution\((\varphi ,\mu )\) to the system (4.1) in the sense of Definition 4.6.

The idea behind the proof of Theorem 4.7 is to approximate the double-obstacle potential by a sequence\((F_n)_{n\in {\mathbb {N}}}\) of regular potentials where for each\(n\in {\mathbb {N}}\),\(F_n\) is a regular potential fulfilling the condition F2. Therefore, Corollary 4.4 can be applied to derive a suitable uniform bound on the terms involving\(F_n'\). We point out that the same strategy could be used to construct a weak solution to the diffuse-interface model (4.1) in the case thatF is the logarithmic potential (1.2b).

4.3Proofs

4.3.1Proof of Theorem 4.3

The proof is divided into five steps.

Step 1: Implicit time discretization. Let\(N\in {\mathbb {N}}\) be arbitrary. We define\(\tau := T/N\) as our time step size. Let now\(n\in \{0,...,N-1\}\) be arbitrary. We now define functions\(\varphi ^n\) with\(n=0,...,N\) by the following recursion:

• The zeroth iterate is defined as the initial datum, i.e.,\(\varphi ^0 := \varphi _0\).

• If for some\(n\in \{0,...,N-1\}\), then-th iterate\(\varphi ^n\) is already constructed, we choose\(\varphi ^{n+1}\in H^1_{(m)}\) as a minimizer of the functional

  $$\begin{aligned} J_n: H^1_{(m)}(\Omega ) \rightarrow {\mathbb {R}}, \quad J_n(\varphi ) := \frac{1}{2\tau }\left\| \varphi -\varphi ^n \right\| _{S_a}^2 + {\mathcal {E}}(\varphi ). \end{aligned}$$

  (4.21)

  Here,\({\mathcal {E}}\) is the energy functional defined in (4.2), with\({\varepsilon }=\sigma ={c_{_{F}}}=1\), and\(\left\| \cdot \right\| _{S_a}\) is the norm defined in (4.11) witha being chosen as

  $$\begin{aligned} a := M(\nabla \varphi ^n,\varphi ^n). \end{aligned}$$

  (4.22)

  This choice is actually possible since the functionM is assumed to be bounded and uniformly positive (see A4). The existence of a minimizer of the functional\(J_n\) will be established in Step 2.

The idea behind this construction is that the first variation in the functional\(J_n\) at the point\(\varphi ^{n+1}\) is zero since\(\varphi ^{n+1}\) is a minimizer of\(J_n\). This means that

for all test functions\(\eta \in H^1_{(0)}(\Omega )\). We now define

with

anda being chosen as in (4.22). Recalling the definition of the inner product\( \left( \cdot \,{,}\, \cdot \right) _{S_a}\) (see (4.10)), we infer from (4.23) that

for all\(\eta \in H^1_{(0)}(\Omega )\). Due to the choice of the constant\(c^{n+1}\), a straightforward computation reveals that (4.26) remains true even for all test functions\(\eta \in H^1(\Omega )\). This means that for every\(n\in \{0,...,N-1\}\), the pair\((\varphi ^{n+1},\mu ^{n+1})\) satisfies the equations

for all test functions\(\zeta ,\eta \in H^1(\Omega )\). Here, (4.27a) follows directly from the construction of\(\mu ^{n+1}\) in (4.24) and the definition of the solution operator\(S_a\) (see (4.9)). The system (4.27) can be interpreted as a time-discrete approximation of the weak formulation (4.12).

The time-discrete approximate solution now needs to be extended onto the whole time interval [0, T]. Thepiecewise constant extension\((\varphi _N,\mu _N)\) is defined as

whereas thepiecewise linear extension\((\overline{\varphi }_N,\overline{\mu }_N)\) is defined as

for\(t = \lambda n\tau + (1-\lambda )(n-1)\tau \) with\(n\in \{1,...,N\}\) and\(\lambda \in [0,1]\).

Henceforth, the letterC will denote generic positive constants that may depend only on\(\varphi _0\) and the constants introduced in A2–A4 and F1 but not onn,N or\(\tau \). These constants may also change their value from line to line.

Step 2: Existence of a minimizer to the functional\(J_n\). We now prove that the functional\(J_n\) introduced in (4.21) actually possesses a minimizer. Therefore, we employ the direct method of the calculus of variations.

For any\(\varphi \in H^1_{(m)}(\Omega )\), we obtain

by means of Poincaré’s inequality. This directly implies

for some positive constant\(c_*\) depending only onm and\(\Omega \). Recalling the assumptions onA (see A3), that\(F\ge -B_F\) (see F1) and that\(G\ge 0\) (see A2), we use Poincaré’s inequality to derive the estimate

for all\(\varphi \in H^1_{(m)}(\Omega )\). This means that\(J_n\) is coercive and bounded from below. Hence, the infimum

exists, and consequently, there also exists a corresponding minimizing sequence\((\varphi _k)_{k\in {\mathbb {N}}}\) with

Now, (4.31) directly implies that\((\varphi _k)_{k\in {\mathbb {N}}}\) is bounded in\(H^1_{(m)}(\Omega )\). Using the Banach–Alaoglu theorem, the compact embeddings\(H^1_{(m)}(\Omega ) \hookrightarrow L^p(\Omega )\) and\(H^1_{(m)}(\Omega ) \hookrightarrow L^q(\partial \Omega )\), we infer that there exists a function\(\overline{\varphi }\in H^1_{(m)}(\Omega )\) such that

along a non-relabeled subsequence. SinceA is continuous and convex (see A3), we infer

due to weak lower semicontinuity. Recalling the growth conditions onF andG (see F1 and A2) and the convergences in (4.32), we apply Lebesgue’s general convergence theorem (see Alt2016, [Section 3.25]) to conclude

as\(k\rightarrow \infty \). Combining (4.33) and (4.34), we obtain

This proves that\(\overline{\varphi }\) is a minimizer of the functional\(J_n\).

Step 3: A priori estimates for the piecewise constant extension. We now claim that the piecewise constant extension\((\varphi _N,\mu _N)\) fulfills the uniform priori estimate

To prove (4.35), we exploit the recursive construction of the time-discrete approximate solution. Since for any\(n\in \{0,...,N-1\}\),\(\varphi ^{n+1}\) was chosen to be a minimizer of the functional\(J_n\), we have

for all\(n\in \{0,...,N-1\}\). By a simple induction, we thus infer

Recalling the assumptions onA (see A3) and that the potentialsF andG are bounded from below (see A2 and F1), we use estimate (4.30) and (4.37) to obtain

for all\(n\in \{0,...,N-1\}\). By the definition of\(\varphi _N\), this directly implies

For any\(n\in \{1,...,N\}\), we now set\(t_n := n\tau \). By the definition of the piecewise constant extension, we have

for all\(t\in (t_{n-1},t_n]\). Recalling the priori estimate (4.36) and the definition of\(\mu ^{n}\) (see (4.24)), we obtain

for all\(n\in \{0,...,N-1\}\). Hence, by induction, we get

for all\(n\in \{0,...,N-1\}\). Now, for any\(t\in (0,T]\), we find an index\(n\in \{0,...,N-1\}\) such that\(t\in (t_{n-1},t_n]\). Recalling (4.40), we eventually conclude that

for all\(t\in [0,T]\). In particular, choosing\(t=T\), we obtain the uniform bound

We now test (4.27b) with the constant function\(\eta \equiv 1/\left| \Omega \right| \). Using the growth assumptions from F1, the continuous embeddings\(H^1(\Omega ) \hookrightarrow L^5(\Omega )\) and\(H^1(\Omega ) \hookrightarrow L^3(\partial \Omega )\), as well as the uniform bound (4.39), we derive the estimate

Applying Poincaré’s inequality, we thus obtain

Combining (4.42) and (4.43), this yields

Due to (4.39) and (4.44), the a priori estimate (4.35) is now established.

Step 4: A priori estimate for the piecewise linear extension. We next claim that for all\(s,t\in [0,T]\),

In particular, the first estimate means that the piecewise linear extension\(\overline{\varphi }_N\) is Hölder continuous in time.

To prove these inequalities, we first infer from (4.27a) and the definition of the piecewise linear extension (see (4.29)) that

for almost all\(\tau \in [0,T]\) and all\(\zeta \in H^1(\Omega )\). Let now\(\xi \in L^2(0,T;H^1(\Omega ))\) be arbitrary. We test (4.46) with\(\xi (\tau )\) and integrate the resulting equation with respect to\(\tau \) from 0 toT. Then, using Hölder’s inequality as well as the a priori estimate (4.35), we obtain

Taking the supremum over all\(\xi \in L^2(0,T;H^1(\Omega ))\) with\(\left\| \xi \right\| _{L^2(0,T;H^1(\Omega ))} \le 1\), this proves estimate (4.45c).

Next, let\(s,t\in [0,T]\) be arbitrary. Without loss of generality, we assume\(s<t\). Integrating (4.46) with respect to\(\tau \) froms tot, choosing\(\zeta = \overline{\varphi }_N(t) - \overline{\varphi }_N(s)\), and using Hölder’s inequality, we derive the estimate

In view of the a priori estimate (4.35), this proves (4.45a).

Let now\(t\in [0,T]\) be arbitrary. Then, we find\(\lambda \in [0,1]\) and\(n\in \{1,...,N\}\) such that\(t=\lambda n \tau + (1-\lambda ) (n-1) \tau \). We thus obtain

Applying (4.45a) with\(t=n\tau \) and\(s=(n-1)\tau \), we conclude (4.45b). This means that all estimates in (4.45) are established.

Step 5: Convergence to a weak solution. In view of the uniform a priori estimate (4.35), the Banach–Alaoglu theorem implies the existence of functions\(\varphi \in L^\infty (0,T;H^1(\Omega ))\) and\(\mu \in L^2(0,T;H^1(\Omega ))\) such that

as\(N\rightarrow \infty \), along a non-relabeled subsequence. We further know that

In combination with the uniform estimate (4.45c), we use the Banach–Alaoglu theorem to infer\(\varphi \in H^1(0,T;H^1(\Omega )')\) with

as\(N\rightarrow \infty \), up to subsequence extraction. Moreover, due to the compact embeddings\(H^1(\Omega ) \hookrightarrow L^p(\Omega )\) and\(H^1(\Omega ) \hookrightarrow L^q(\partial \Omega )\), we apply the Aubin–Lions lemma to obtain

By passing to the limit in estimate (4.45a), we conclude\(\varphi \in C^{0,1/4}([0,T],L^2(\Omega ))\). This means that the functions\(\varphi \) and\(\mu \) satisfy the regularity conditions of Definition 4.2(i). Using the estimate (4.45b), we directly deduce from (4.52) that

as\(N\rightarrow \infty \), after another subsequence extraction.

From the time-discrete weak formulation (4.27), we infer that the piecewise constant extension\((\varphi _N,\mu _N)\) and the piecewise linear extension\((\overline{\varphi }_N,\overline{\mu }_N)\) satisfy the approximate weak formulation

for all test functions\(\xi ,\vartheta \in L^2(0,T;H^1(\Omega ))\). Recalling the growth conditions on\(F'\) and\(G'\) from F1 and A2, as well as the priori estimate (4.35), we infer that the sequence\((F'(\varphi _N))_{N\in {\mathbb {N}}}\) is bounded in\(L^\infty (0,T;L^{6/5}(\Omega ))\) and the sequence\((G'(\varphi _N))_{N\in {\mathbb {N}}}\) is bounded in\(L^\infty (0,T;L^{4/3}(\partial \Omega ))\). Hence, according to the Banach–Alaoglu theorem, there exist functions\(f^* \in L^\infty (0,T;L^{6/5}(\Omega ))\) and\(g^* \in L^\infty (0,T;L^{4/3}(\partial \Omega ))\) such that

as\(N\rightarrow \infty \), along a non-relabeled subsequence. Moreover, the convergences in (4.53) directly imply\(F'(\varphi _N) \rightarrow F'(\varphi )\) a.e. in\(\Omega \) and\(G'(\varphi _N) \rightarrow G'(\varphi )\) a.e. on\(\partial \Omega \). By a convergence principle based on Egorov’s theorem (see DiBenedetto2002, [Proposition 9.2c]), we now infer\(f^* = F'(\varphi )\) a.e. in\(\Omega \) and\(g^* = G'(\varphi )\) a.e. on\(\partial \Omega \). This means that

as\(N\rightarrow \infty \). Testing the approximate weak formulation (4.54b) with\(\vartheta = \varphi _N - \varphi \) and employing the strong monotonicity condition on\(A'\) from A3, we obtain

Using the convergences (4.50), (4.53), (4.55) and (4.56) along with the weak-strong convergence principle, we infer that the right-hand side of the above estimate tends to zero. We thus conclude that

as\(N\rightarrow \infty \), up to subsequence extraction. In view of the growth condition on\(A'\) from A3, Lebesgue’s general convergence theorem further reveals that

Due to the convergences (4.50), (4.55), (4.56) and (4.59), we can now pass to the limit in (4.54b) to conclude that

holds for all\(\vartheta \in L^2(0,T;H^1(\Omega ))\).

We now fix an arbitrary time\(t_0\in (0,T]\). Since\(\tau = T/N \rightarrow 0\) as\(N\rightarrow \infty \), we may assume (without loss of generality) thatN is chosen large enough to ensure\(t-\tau \in [0,T]\) for all\(t\in [t_0,T]\). We have

for almost all\(t\in [t_0,T]\). Here, from the second to the third line, we used the change of variables\(s = t-\tau \) and the fact that\([t_0-\tau ,T-\tau ]\subset [0,T]\) to estimate the first summand. Now, as\(N\rightarrow \infty \), the first summand in the third line of (4.61) tends to zero because of (4.58), whereas the second summand tends to zero since due to mean-continuity in\(L^p(Q)\) (see, e.g., Alt2016, [Section 4.15]). This proves

as\(N\rightarrow \infty \). Since\(t_0\in (0,T]\) was arbitrary, we deduce

as\(N\rightarrow \infty \), after extracting a subsequence. Proceeding similarly, and using the strong convergence\(\varphi _N\rightarrow \varphi \) in\(L^2(Q)\) (which directly follows from (4.53)), we further obtain

as\(N\rightarrow \infty \). Using (4.63) and (4.64) along with Lebesgue’s dominated convergence theorem, we infer

strongly in\(L^2(Q)\), as\(N\rightarrow \infty \), up to subsequence extraction. Employing the weak-strong convergence principle, we can thus pass to the limit\(N\rightarrow \infty \) in the approximate weak formulation (4.54a) to obtain

for all\(\zeta \in L^2(0,T;H^1(\Omega ))\). Combining (4.60) and (4.66), we eventually conclude that the pair\((\varphi ,\mu )\) satisfies the weak formulation (4.12). Moreover, as a direct consequence of the convergence (4.52),\(\varphi \) satisfies the initial condition (4.13). This means that all conditions of Definition 4.2(ii) are fulfilled.

Recalling the growth conditions onF andG from F1 and A2, as well as the convergences in (4.53), we apply Lebesgue’s general convergence theorem (see Alt2016, [Section 3.25]) to conclude

Then, from the convergences (4.58),(4.67) and (4.68), we infer that

as\(N\rightarrow \infty \). Recalling (4.50) and (4.65), we use the weak-strong convergence principle to infer

as\(N\rightarrow \infty \). We now use the convergences (4.69) and (4.70), the weak lower semicontinuity of the\(L^2(Q)\)-norm, as well as the discrete energy inequality (4.41), to derive the estimate

for almost all\(t\in [0,T]\). This proves the weak energy dissipation law (4.14), and thus, the condition in Definition 4.2(iii) is fulfilled.

We eventually conclude that the pair\((\varphi ,\mu )\) is a weak solution to system (4.1) in the sense of Definition 4.2. Hence, the proof is complete.\(\Box \)

4.3.2Proof of Corollary 4.4

Let\((\varphi ,\mu )\) be a weak solution to the system (4.1), whose existence is guaranteed by Theorem 4.3. By a straightforward computation, we notice that

where

Hence, in the following, we intend to prove (4.16) by deriving suitable bounds on the terms\(I_1\) and\(I_2\). The letterC will denote generic positive constants depending only on\(\varphi _0\),\({\mathcal {E}}(\varphi _0)\),\(c_0\) and the constants in A1–A4, but not on\(c_1\).

Let\(\eta \in H^1(\Omega )\) be arbitrary. SinceA is convex (see A3), we know that

Testing the weak formulation (4.12b) with\(\eta -\varphi \) instead of\(\eta \) and using the above estimate, we thus infer that the variational inequality

holds a.e. in [0, T] for all\(\eta \in H^1(\Omega )\). Moreover, since\((\varphi ,\mu )\) is a weak solution of (4.1), it satisfies the weak energy inequality (4.14). Using Poincaré’s inequality, we infer

Step 1: We first derive an estimate for the term\(I_1\). Therefore, we choose

for sufficiently small\(\delta >0\) which ensures\(1 - \delta F''(\varphi ) > 0\). Since\(F'(\varphi ) \in L^\infty (0,T;H^1(\Omega ))\) due to (4.15), we know that\(\eta \in L^\infty (0,T;H^1(\Omega ))\). Recalling thatA is positively homogeneous of degree 2, we obtain

a.e. inQ. We now test the variational inequality (4.73) with\(\eta \). After dividing the resulting inequality by\(\delta \), we use (4.76) to deduce

a.e. in [0, T]. Recalling that F2 implies that F1 holds with\(p=2\), we derive the estimate

a.e. in [0, T]. Hence, using the growth condition on\(G'\) from A2 and the continuous embedding\(H^1(\Omega )\hookrightarrow L^4(\partial \Omega )\), we deduce

a.e. in [0, T]. Sending\(\delta \rightarrow 0\) and using the growth condition from A3, the condition\(-F'' \le c_0\) (cf. (4.15)) as well as Poincaré’s inequality and Young’s inequality, we infer

a.e. in [0, T]. Integrating this inequality with respect to time from 0 toT, and using estimate (4.74), we eventually conclude the bound

Step 2: We now derive a suitable estimate for the term\(I_2\). Let\(\lambda \in L^\infty ([0,T])\) be any function that will be fixed later. We set

for some\(\delta >0\). Testing the variational formulation with this\(\eta \), dividing the resulting equation by\(\delta \), and recalling thatA is positively homogeneous of degree 2, we derive the estimate

Since\(F'' +c_0 \ge 0\) due to (4.15), we know that the function\(s\mapsto F(s) + \tfrac{1}{2} c_0s^2 \) is convex. We thus have

a.e. inQ. Using this estimate, as well as Young’s inequality, we now get

almost everywhere in [0, T]. Sending\(\delta \rightarrow 0\) in (4.80) and using the above estimate, the growth conditions from A2 and A3, the continuous embedding\(H^1(\Omega )\hookrightarrow L^4(\partial \Omega )\), as well as Poincaré’s inequality, we infer

We now fix\(\lambda \) as

for all\(t\in [0,T]\). Testing (4.12a) with\(\zeta \equiv 1\) and integrating the resulting equation with respect to time, we infer\(\left\langle \varphi (t) \right\rangle _\Omega = \left\langle \varphi _0 \right\rangle _\Omega \) for all\(t\in [0,T]\). In view of (4.82), we thus get

a.e. in [0, T]. We now multiply this estimate by\(\tfrac{2}{\kappa }\) and take the square on both sides. Integrating the resulting inequality with respect to time and using the uniform estimate (4.74), we eventually conclude

We finally plug the estimates (4.78) for\(I_1\) and (4.85) for\(I_2\) into (4.72). This proves (4.16), and thus, the proof of Corollary 4.4 is complete.\(\Box \)

4.3.3Proof of Theorem4.7

The proof is split into three steps.

Step 1: Approximation of the double-obstacle potential by smooth potentials. To prove the assertion, we approximate the double-obstacle potentialF by a sequence of regular potentials\((F_n)_{n\in {\mathbb {N}}}\). Therefore, we define the function

and for any\(n\in {\mathbb {N}}\), we set

By this construction, we have\(J\in C^2({\mathbb {R}};[0,\infty ))\),J is convex, and\(F_n = F_0\) on\([-1,1]\) for all\(n\in {\mathbb {N}}\). It is straightforward to check that for all\(n\in {\mathbb {N}}\), the approximate potential\(F_n\) satisfies the assumption F2 with\(c_0 = 1\) and\(c_1 = 12n\). It thus follows that F1 is satisfied with\(p=2\) and\(B_F = \tfrac{3}{2}\). In the remainder of this proof, it will be crucial that the constants\(B_F\) and\(c_0\) are independent ofn. For any\(n\in {\mathbb {N}}\), we further define the approximate energy functional by

Step 2: A priori estimates for the sequence of approximate solutions. We now conclude from Theorem 4.3 that for every\(n\in {\mathbb {N}}\), there exists a weak solution\((\varphi _n,\mu _n)\) of the system (4.1) to the potential\(F_n\) in the sense of Definition 4.2. In the following, the letterC will denote generic positive constants that may depend on\(\varphi _0\),\(\kappa \) and the constants in A1–A4 but not on the approximation indexn.

As the weak solutions\((\varphi _n,\mu _n)\) satisfy the weak energy dissipation law (4.14) written for\({\mathcal {E}}_n\), we deduce the estimate

for almost all\(t\in [0,T]\) and all\(n\in {\mathbb {N}}\). As\(B_F\) is independent ofn, we use Poincaré’s inequality to conclude the uniform bound

Integrating the weak formulation (4.12a) written for\((\varphi _n,\mu _n)\) with respect to time from 0 toT, we now use (4.88) to derive the uniform estimate

Furthermore, Corollary 4.4 provides the estimate

where\( R= \left| \left\langle \varphi _0 \right\rangle _\Omega \right| + \frac{\kappa }{2} < 1\). Here, the constantc depends only on\(\varphi _0\),\(\mathcal E_n(\varphi _0)\),\(c_0=1\) and the constants in A1–A4. Since\(F_n = F_0\) on\([-1,1]\), we know that\(F_n(\varphi _0) = F_0(\varphi _0)\) for all\(n\in {\mathbb {N}}\). Consequently,\({\mathcal {E}}_n(\varphi _0)\) does not depend onn, and thus,c is independent ofn. We infer the uniform bound

Using (4.88), we further get

Combining (4.91) and (4.92), we now conclude

Step 3: Convergence to a weak solution. In view of the uniform estimates (4.88) and (4.89), we now use the continuous embedding\(H^1(\Omega )\hookrightarrow L^4(\partial \Omega )\), the Banach–Alaoglu theorem, and the Aubin–Lions lemma along with the compact embeddings\(H^1(\Omega )\hookrightarrow L^2(\Omega )\) and\(H^1(\Omega )\hookrightarrow L^r(\partial \Omega )\) for\(r\in [1,4)\) to conclude the existence of functions\(\varphi \) and\(\mu \) such that

for all\(r\in [1,4)\) as\(n\rightarrow \infty \), after extraction of a subsequence. Using the uniform bound (4.93) along with the Banach–Alaoglu theorem, as well as the pointwise–a.e. convergence stated in (4.95), we deduce

as\(n\rightarrow \infty \). As the strong limit in\(L^2(Q)\) and the pointwise limit coincide, we have\(J'(\varphi ) = 0\) a.e. inQ. Since\(J'(s)=0\) if\(\left| s \right| \le 1\) and\(\left| J'(s) \right| >0\) if\(\left| s \right| >1\), we conclude

As\(F_0'(\varphi _n) = - \varphi _n\), the convergence

follows directly from (4.95). Moreover, using the growth condition on\(G'\) (see A2), (4.95) and Lebesgue’s general convergence theorem (see Alt2016, [Section 3.25]), we obtain

as\(n\rightarrow \infty \), after another subsequence extraction. Arguing as in the proof of Theorem 4.3, we exploit the strong monotonicity condition on\(A'\) from A3 to further derive the convergences

as\(n\rightarrow \infty \), up to subsequence extraction. Combining (4.95) and (4.100), we eventually get

Let now\(n\in {\mathbb {N}}\) be arbitrary and let\(\zeta \in H^1(\Omega )\) and\(\eta \in L^2(0,T;H^1(\Omega ))\) with\(\left| \eta \right| \le 1\) a.e. inQ be an arbitrary test functions. This already implies that\(J'(\eta ) = 0\) a.e. inQ. Moreover, sinceJ is convex its derivative\(J'\) is monotonically increasing. We thus have

We now recall that the weak solution\((\varphi _n,\mu _n)\) satisfies the weak formulation (4.12) written for\((\varphi _n,\mu _n)\) instead of\((\varphi ,\mu )\). The weak formulation (4.12a) written for\((\varphi _n,\mu _n)\) and tested with\(\zeta \) reads as

Testing the weak formulation (4.12b) written for\((\varphi _n,\mu _n)\) with\(\varphi _n-\eta \), integrating with respect to time from 0 toT, and employing estimate (4.103), we obtain

Using the convergences (4.94)–(4.102), Lebesgue’s dominated convergence theorem, as well as the weak-strong convergence principle, we pass to the limit\(n\rightarrow \infty \) in (4.104) and (4.105). This proves that the pair\((\varphi ,\mu )\) satisfies the weak formulation (4.18a) for all\(\zeta \in H^1(\Omega )\), as well as the variational inequality (4.18b) for all\(\eta \in L^2(0,T;H^1(\Omega ))\) with\(\left| \eta \right| \le 1\) a.e. inQ. Moreover, (4.95) directly implies that\(\varphi \) satisfies the initial condition (4.19). This means that all conditions of Definition 4.6(ii) are verified.

Proceeding similarly as in Step 4 of the proof of Theorem 4.3, and using the weak formulation (4.18a), we can show a posteriori that\(\varphi \) is Hölder continuous in time in the sense that\(\varphi \in C^{0,1/4}([0,T];L^2(\Omega ))\). In combination with (4.94)–(4.96), this proves that all conditions of Definition 4.6(1) are fulfilled.

Recalling that\(\left| \varphi \right| \le 1\) a.e. inQ, we use (4.95) along with Lebesgue’s general convergence theorem (see Alt2016, [Section 3.25]) to conclude

a.e. in [0, T]. Using the convergences (4.95), (4.96), (4.100) and (4.102), we now proceed similarly as in Step 5 of the proof of Theorem 4.3 (cf. (4.71)) to verify that the pair\((\varphi ,\mu )\) satisfies the weak energy dissipation law (4.20). This means that Definition 4.6(iii) is also fulfilled.

In summary, we conclude that the pair\((\varphi ,\mu )\) is a weak solution to system (4.1) (withF being the double-obstacle potential) in the sense of Definition 4.6. Hence, the proof of Theorem 4.7 is complete.\(\Box \)

In this section, we present numerical comparisons between the diffuse-interface model (3.3) and its sharp-interface limit (3.52).

For the sharp-interface computations, (SI), we employ the parametric finite element approximation from Bao et al. (2023), which uses piecewise linear finite elements and relies crucially on the stable approximation of the anisotropy introduced in Barrett et al. (2008a,2008b), see also Bao and Zhao (2022). Here, we recall that this stable approximation is designed for anisotropy functions of the form

where\(\Lambda _\ell \),\(\ell =1,\ldots ,L\), are symmetric and positive definite matrices. We refer to Barrett et al. (2008a), Barrett et al. (2008b), Barrett et al. (2010b), Bao and Zhao (2022), Bao et al. (2023) for details. Clearly, for (5.1), the assumption A3 is satisfied, recall Remark 4.1.

For the diffuse-interface approximations, (DI), we adapt the finite element discretizations from Barrett et al. (2014) to the system (3.3). To this end, we assume that\(\Omega \) is a polyhedral domain and let\({\mathcal {T}}_{h}\) be a regular triangulation of\(\Omega \) into disjoint open simplices. Associated with\({\mathcal {T}}_h\) is the piecewise linear finite element space

where we denote by\(P_{1}(o)\) the set of all affine linear functions ono, cf. Ciarlet (1978). We also let\((\cdot ,\cdot )\) denote the\(L^{2}\)-inner product on\(\Omega \), and let\((\cdot ,\cdot )^{h}\) be the usual mass lumped\(L^{2}\)-inner product on\(\Omega \) associated with\({\mathcal {T}}_{h}\). In a similar fashion, we let\(\langle \cdot ,\cdot \rangle _{\Gamma _w}^h\) denote the mass lumped\(L^{2}\)-inner product on\(\Gamma _w\). Finally,\(\Delta t\) denotes a chosen uniform time step size.

Our fully discrete finite element approximation of (3.3) is then given as follows. For\(n\ge 0\), let\(\varphi _h^{n} \in S^h\) be given. Then, find\((\varphi _h^{n+1}, \mu ^{n+1}_h) \in S^h \times S^h\) such that

for all\((\chi , \eta ) \in S^h \times S^h\). The above scheme utilizes the linearization\(B(\mathbf {p}) \mathbf {p} = A'(\mathbf {p})\) for anisotropies of the form (5.1), which was first introduced in Barrett et al. (2013). In particular, the symmetric positive definite matricesB are defined by

We stress that the induced semi-implicit discretization of\(A'(\nabla \varphi )\) in (5.2b) ensures that our numerical method is stable. In fact, using the techniques in Barrett et al. (2013), Barrett et al. (2014), and on employing semi-implicit approximations of\(F'(\varphi )\) and\(G'(\varphi )\) based on convex/concave splittings ofF andG, an unconditional stability result can be shown. However, for the purposes of this paper, we prefer the simpler approximation (5.2). We also note that extending the scheme (5.2) to the case of the obstacle potential (1.2c), when (5.2b) needs to be replaced with a variational inequality, is straightforward. We refer to Barrett et al. (2013), Barrett et al. (2014) for the precise details.

We implemented the scheme (5.2), and its obstacle potential variant, with the help of the finite element toolbox ALBERTA, see (Schmidt and Siebert2005). To increase computational efficiency, we employ adaptive meshes, which have a finer mesh size\(h_{f} = \frac{\sqrt{2}}{N_f}\) within the diffuse interfacial regions and a coarser mesh size\(h_{c} = \frac{\sqrt{2}}{N_c}\) away from them, with\(N_f, N_c \in {\mathbb {N}}\), see Baňas and Nürnberg (2008), Barrett et al. (2004) for a more detailed description. The nonlinear systems of equations arising from (5.2) at each time step are solved with a Newton method, where we employ the sparse factorization package UMFPACK, see Davis (2004), for the solution of the linear systems at each iteration. In the case of the double-obstacle potential, we employ the solution method from Baňas and Nürnberg (2008), Barrett et al. (2014).

In all our computations, we fix the mobility\(D({\varvec{\nu }}) =1\) and, up to possible rotations, use the anisotropy

which can be regarded as a smoothed\(\ell ^1\)-norm, with a small regularization parameter\(\delta >0\). Note that, (5.4) is a special case of (5.1). For the (DI) computations, we choose for the potentialF either (1.2a), so that\({c_{_{F}}}= \frac{4}{3}\), or (1.2c), so that\({c_{_{F}}}= \frac{\pi }{2}\). We letG be defined by (2.9), while the regularized mobility functions are defined via (3.1) and (3.2), with\(r=2\) and\(\gamma _0 = d_1 = 1\). We also choose\(\alpha =\frac{{c^2_{_{F}}}}{4}\) so that (3.25) is consistent with (2.1a). Finally, unless otherwise stated, we use the smooth potential (1.2a) for our (DI) computations.

Evolution of small island films toward the equilibrium (red line) for the SI approximations.a Plots at\(t=0,0.002,0.01,0.02,0.030878,0.0319, 0.0339, 0.1\), where the island occurs pinch-off at\(t=0.030878\);b–d are the plots at\(t=0,0.01,0.02,\cdots ,0.1\) (Color figure online)

Full size image

5.12d Results

In numerical simulations of solid-state dewetting problems, it is often of interest whether a thin film of material breaks up into islands. For example, in two space dimensions and in the isotropic case with a\(90^\circ \) contact angle condition, it has been observed that elongated films undergo pinch-off once the aspect ratio of length versus height goes beyond a critical value\(R_0 \approx 127.9\), (Dornel et al.2006; Wang et al.Jan 2015). For nonzero values of\(\sigma \), the critical value behaves like\(R_0 \approx 96.6/\sin (\frac{1}{2}\arccos \sigma )-8.66\), (Dornel et al.2006).

it turns out that the anisotropy\(\gamma \) can have a dramatic influence on the critical value\(R_0\). To investigate this numerically, we simulate the evolution of small thin films, starting from an initial interface in the form of the upper half of a tube with aspect ratio\(R = L / H\), and fix\(H = 0.3\). We consider the following four example situations:

1. (a)

   an island of\(R=15\) with anisotropy\(\gamma ({{\mathcal {R}}}(\frac{\pi }{4})\mathbf {p})\) and\(\sigma =\cos \frac{5\pi }{6}\);

2. (b)

   an island of\(R=15\) with anisotropy\(\gamma (\mathbf {p})\) and\(\sigma =\cos \frac{5\pi }{6}\);

3. (c)

   an island of\(R=15\) with anisotropy\(\gamma ({{\mathcal {R}}}(\frac{\pi }{4})\mathbf {p})\) and\(\sigma =\cos \frac{\pi }{2}\);

4. (d)

   an island of\(R=13\) with anisotropy\(\gamma ({{\mathcal {R}}}(\frac{\pi }{4})\mathbf {p})\) and\(\sigma =\cos \frac{5\pi }{6}\),

where\({{\mathcal {R}}}(\theta )\) is the rotation matrix with an angle\(\theta \), and\(\gamma (\mathbf {p})\) is given by (5.4) with\(d=2\),\(\delta = 0.1\). We note that anisotropies with a fourfold symmetry like our choices above are often used in two-dimensional models for materials with a cubic crystalline surface energy (Liu and Metiu1993; McFadden et al.2000; Zhang and Gladwell2003).

Plots of the interface profiles for the SI approximations are presented in Fig. 4a–d for the four examples, respectively, where the approximated polygonal curve consists of 2048 line segments, and the time step size is fixed as\(10^{-6}\). From these figures, we can observe the influence of the anisotropy\(\gamma \), the contact energy density difference\(\sigma \), and the aspect ratioR of the thin film on the evolution. In particular, comparing the evolutions in Fig. 4a and d, we see that the critical value\(R_0\) for break-up to occur appears to satisfy\(13 < R_0 \le 15\), which is much smaller than in the isotropic case. Moreover, we see that either rotating the anisotropy, Fig. 4b, or changing the contact angle, Fig. 4c, ensures that no break-up occurs, meaning that\(R_0 > 15\) in both cases.

Let us remark that the pinch-off observed in Fig. 4a represents a singularity for the parametric description on which the SI approximations are based. Hence, we perform a heuristical topological change, from a single curve to two separate curves, once an inner vertex of the polygonal curve touches the substrate. In what follows, we will use the computations in Fig. 4 as reference solutions for our DI approximations, in order to empirically confirm our theoretical results from Sect. 3.

The time history of the energy for the DI and SI approximations in the four different examples using the double smooth potential

Full size image

Left panel: The time history of the energy for the DI and SI approximations in Example (d) using the obstacle potential (1.2c). Right panel: The errors\({\mathcal {E}}_\Delta \) of the energy at the final time\(T=0.1\) between the SI and DI approximations plotted against\({\varepsilon }\). Here, “(d)-obstacle” refers to Example (d) with the obstacle potential, while the remaining graphs are for Examples (a)–(d) with the smooth potential

Full size image

For our DI approximations, we consider the computational domain\(\Omega = [0,3]\times [0,1]\), on which for symmetry reasons, we only compute the right half of the evolving thin film. As interfacial parameters, we consider\({\varepsilon }= 1 / (2^{4+i}\pi )\), for\(i=0,\ldots ,2\), and choose the discretization parameters as\(N_f = 2^{8+i}\),\(N_c = 2^{5+i}\),\(\Delta t = 10^{-3} / 2^{4+2i}\). These spatial adaptive discretization parameters allow for a sufficient resolution of the diffuse-interface, while the temporal discretization parameters yield an excellent agreement with the SI approximations. In fact, in Fig. 5, we show the energy plots of the DI approximations and compare them with the corresponding SI approximations for the four different examples from Fig. 4. We observe that for sufficiently small values of\({\varepsilon }\), there is excellent agreement between the SI and DI evolutions, in line with our asymptotic analysis in Sect. 3. What is interesting to note is that for Example (a), the pinch-off time predicted by the DI computations is too early when\({\varepsilon }\) is not small, and this can be explained by the fact that the wider interfacial region “sees” contact with the substrate earlier, leading to the break-up into two islands. For the same reason, in Examples (c) and (d), the DI computations for\({\varepsilon }=1/(16\pi )\) erroneously predict a pinch-off, leading to a larger final energy. But once\({\varepsilon }\) is sufficiently small, no pinch-off occurs, in agreement with the SI evolutions.

We note that using the double-obstacle potential (1.2c) leads to very similar results. As an example, we show the evolution of the discrete energies for Example (d) in Fig. 6. In addition, in order to also have a quantitative comparison between our SI and DI computations, in the same figure, we also present plots of the energy difference\({\mathcal {E}}_\Delta \) between the final SI and DI solutions against\({\varepsilon }\). The presented results suggest that the DI energies of the final states approach the corresponding SI energy with\(O({\varepsilon })\). Note that, the three instances where\({\mathcal {E}}_\Delta \ge 10^{-1}\) correspond to cases where the DI computations wrongly predict a pinch-off. Moreover, in practice, we observe that the contact angles between DI and SI at the final time agree very well, with the error being of order\(10^{-3}\) throughout.

The qualitative behavior of the DI and SI approximations is compared in Figs. 7,8,9,10. In all four examples, we note an excellent agreement between the two different approaches. This is particularly noteworthy in Example (a) with the occurrence of a topological change, which is not covered by our asymptotic analysis.

[Example (a)] Interface profiles at times\(t=0,0.01,0.02,0.03,0.04,0.1\) for the DI approximations with\({\varepsilon }= 1/(64\pi )\), and the red dash line represents the SI approximations (Color figure online)

Full size image

[Example (b)] Interface profiles at times\(t=0,0.01,0.02,0.03,0.04,0.1\) for the DI approximations with\({\varepsilon }= 1/(64\pi )\), and the red dash line represents the SI approximations (Color figure online)

Full size image

[Example (c)] Interface profiles at times\(t=0,0.01,0.02,0.03,0.04,0.1\) for the DI approximations with\({\varepsilon }= 1/(64\pi )\), and the red dash line represents the SI approximations (Color figure online)

Full size image

[ Example (d)] Interface profiles at times\(t=0,0.01,0.02,0.03,0.04,0.1\) for the DI approximations with\({\varepsilon }= 1/(64\pi )\), and the red dash line represents the SI approximations (Color figure online)

Full size image

5.23d Results

In 3d, we compare our SI and DI approximations for the evolution of an initially spherical island for the anisotropy\(\gamma ({{\mathcal {R}}}_x(\frac{\pi }{4}){{\mathcal {R}}}_y(\frac{\pi }{4})\mathbf {p})\), where\(\gamma (\mathbf {p})\) is given by (5.4) with\(d=3\),\(\delta =0.1\), and where\({{\mathcal {R}}}_x(\theta ), {{\mathcal {R}}}_y(\theta )\) are rotation matrices which rotate a vector through an angle\(\theta \) within the (y, z)- and (x, z)-planes, respectively. The initial interface is chosen to be a semisphere of radius 0.4, attached to the (x, y)-plane, and we let\(\sigma = \cos (\frac{5\pi }{6})\).

For the SI computation, we consider a polyhedral surface with 8256 triangles and 4225 vertices, and a time step size\(10^{-4}\). For our DI approximations, on the other hand, we consider the computational domain\(\Omega = [-\frac{1}{2},\frac{1}{2}]^3\) and as interfacial parameters consider\({\varepsilon }= 1 / (2^{2+i}\pi )\), for\(i=0,\ldots ,2\), with the corresponding discretization parameters\(N_f = 2^{5+i}\),\(N_c = 2^{2+i}\),\(\Delta t = 10^{-3} / 2^{2i}\). In Fig. 11, we show the energy plots of the DI approximations and compare them with the corresponding SI simulation, noting once again an excellent agreement when\({\varepsilon }\) is sufficiently small. We also present a plot of the error in the energy between the DI and SI approximations against\({\varepsilon }\). Note that, the large error for\({\varepsilon }= 1/(4\pi )\) is due to that DI simulation wrongly predicting a pinch-off.

Moreover, a qualitative comparison between the evolutions of the interface for both approaches is shown in Fig. 12. In particular, at the bottom of Fig. 12, we see that the sharp-interface approximation agrees very well with the zero level set from the DI computation, underlining once more our asymptotic analysis in Sect. 3.

Left panel: The time history of the energy for the DI and SI approximations for the semisphere experiment in 3d. Right panel: The error\({\mathcal {E}}_\Delta \) of the energy at the final time between the DI and SI approximations plotted against\({\varepsilon }\)

Full size image

A visualization of the zero level sets of the DI approximations for\({\varepsilon }=(16\pi )^{-1}\) at times\(t=0,0.01,0.1\), together with a slice through the adaptive mesh. Below a comparison between the DI and the SI computation at time\(t=0.01\)

Full size image