A vicinal surface model for epitaxial growth with logarithmic free energy

We study a continuum model for solid films that arises from the modeling of one-dimensional step flows on a vicinal surface in the attachment-detachment-limited regime. The resulting nonlinear partial differential equation, $u_t = -u^2(u^3+\alpha u)_{hhhh}$, gives the evolution for the surface slope $u$ as a function of the local height $h$ in a monotone step train. Subject to periodic boundary conditions and positive initial conditions, we prove the existence, uniqueness and positivity of global strong solutions to this PDE using two Lyapunov energy functions. The long time behavior of $u$ converging to a constant that only depends on the initial data is also investigated both analytically and numerically.


1.
Introduction. Epitaxial growth of crystal surfaces below the roughing temperature has attracted extensive interest. Unlike traditional modeling for fluid and solid mechanics where one starts with continuum theories for macroscopic variables, the modeling of crystal films was first established from the atomic perspective. At the nanoscale, crystal surfaces consist of basic structures such as interacting line defects (steps) and flat surface regions (facets). With adatoms detaching from one step and reattaching to another step after traveling along the facets, a step flow describes the mass transport along the crystal surface. For broader physical surveys of crystal growth we refer readers to [3,10,21].
From both discrete and continuum viewpoints, different models have been constructed to characterize step flows. From the discrete perspective, the dynamics of the steps are described by the step velocities using a Burton-Cabrera-Frank (BCF) type framework [3], which was investigated by Ozdemir and Zangwill [20]. The motion of these individual steps are usually modeled by systems of differential equations for step locations. At the macroscopic scale the continuum limit of step evolution is generally reformulated as nonlinear PDEs [5, 9, 12, 17-20, 22, 26]. Based on the conservation of mass, the dynamic equation for the surface height of a solid film, h(t, x), can be written as where the mobility function M (∇h) is a functional of the gradient of h and G(h) represents a surface energy. The mobility function takes distinctive forms in different limiting regimes. In the diffusion-limited (DL) regime, where the dynamics is dominated by the diffusion across the terraces, M is a constant, M ≡ 1. While in the attachment-detachmentlimited (ADL) case, the dominant processes are the attachment and detachment of atoms at step edges, and the mobility function takes the form [1] M (∇h) = |∇h| −1 . (1.2) One form of the surface energy associated with the dynamical model (1.1) was identified [4] as where the first term represents the step-formation energy and leads to a conical singularity at ∇h = 0, and the second term is the step-interaction energy and is consistent with the discrete model for the crystal surface which will be discussed later. In the DL regime, Giga and Kohn [14] rigorously showed that with periodic boundary conditions on h, finite-time flattening to a spatially-uniform solution, h ≡ C, occurs for α = 0. A heuristic argument provided by Kohn [11] indicates that the flattening dynamics is linear in time. However, in the ADL regime with the nonlinear mobility given by (1.2), the dynamics of the surface height equation (1.1) is still an open question [11].
In the ADL case, using (1.2) and (1.3), the evolution equation for the surface dynamics becomes To simplify the problem, consider a one-dimensional monotone vicinal surface corresponding to a monotone step train in the discrete model. Without loss of generality, assume that h x > 0 for the whole domain. Under this assumption the α term drops out and (1.4) can be rewritten as (1. 5) This equation can also be derived from the dynamic equation (1.1) using the surface energy (1.3) with α = 0. Following the method introduced by Ozdemir and Zangwill [20] and Shehadeh et al. [1], the new variable u(t, h) = h x (t, x) > 0 can be adopted to rewrite the PDE as an evolution equation for the surface slope, This is a fourth-order degenerate parabolic PDE and the existence of global weak solutions for this problem has been shown previously in [8]. Specifically, it was proved that the solution u is positive almost everywhere and that u converges to a constant solution in the limit of long times. The central difficulty is how to deal with singularities if u touches down to zero. 1 A method for degenerate parabolic equations used in [2] was adapted by neglecting a zero measure set P T = {h ∈ [0, 1]; u(h) = 0}, which is a closed set so that we can define a distribution on (0, T ) × (0, 1)\P T and formulate the definition of weak solutions to (1.6). Using this definition, the existence and almost-everywhere positivity of global weak solutions of (1.6) was proved. Note that the slope equation (1.6) can also be derived from discrete models in the BCF framework [3] by showing the convergence of the discrete model to its continuum limit [7]. Let {x i (t), i ∈ Z} be the step locations at time t and let the height of each step be a constant a = H/N , where x i+N (t) − x i (t) = L is specified for a train of steps with a length scale L, and H represents the total height of N steps. In the ADL regime, the BCF model without deposition flux is expressed by the step-flow differential equations, where µ i is a chemical potential giving the gradient of the free energy with respect to changes in the steps, µ i = ∂E/∂x i . When the free energy involves only local contributions due to the interaction among steps, f (r), it can be written as Then, defining the step slopes u i (t) = a/(x i+1 (t) − x i (t)), we obtain the differential equations for slopes u i for i = 1, 2, · · · , N − 1, (1.7) Note that the right hand side of (1.7) is equivalent to a centered finite difference discretization in the limit of the step height a → 0, or equivalently, as the number of steps N → ∞. Therefore the solution of the slope ODE (1.7) should converge to the solution u(t, h) of the continuum model where the step slope u is considered as a function of h. For vicinal surface models with entropic and elastic-dipole interactions, the continuum model (1.8) is equivalent to the PDE (1.6) in [8] with the local contribution f being f 1 (r) = 1 2 r −2 . In 1 In [8] the solution u(t, h) of (1.6) was regularized by a small O( ) perturbation, and it has been shown that the solution has a positive lower bound u (h) > . However, the main obstacle remains as the perturbation approaches zero. Since after taking the limit → 0 it is only guaranteed that u(h) > 0 almost everywhere, we can not prevent u from touching down to zero where u −1 tends to infinity. More explicitly, 1/u is in L 1 space, which is non-reflexive, and there is no weak compactness in L 1 space. Thus as we take the limit → 0, 1 u is in Radon space M([0, 1]), rather than L 1 space. Hence we can not obtain 1 u t = (u 3 ) hhhh by directly taking the limit in the regularized problem. [26] investigated a continuum model in the DL regime with nonlocal contributions with f 2 (r) = − ln |r|.

2002, Xiang
In order to improve the results for (1.6) and explore other interesting dynamics in solid films, we modify the surface energy (1.3) to incorporate a logarithmic factor, This modification is motivated by kinetic theory and related energy techniques, and the contribution from the new energy enables us to gain weak compactness in the proof of global strong solutions. This modified surface energy is comparable to (1.3) since the logarithmic correction is negligible for small surface gradients. With this new energy, the one-dimensional evolution equation, restricted to the case when h x > 0 on 0 ≤ x ≤ L, is This evolution equation is comparable to (1.5) except for the logarithmic term due to the difference between the energy functionals (1.3) and (1.9). The rate of dissipation of the surface energy (1.9) for this model is Applying the change of variables u(t, h) = h x (t, x) to (1.10) yields the slope equation This one-dimensional fourth-order nonlinear PDE (1.11) is the main focus of this paper. Note that the previously studied model (1.6) corresponds to the case α = 0 in (1.11). Moreover, the continuum model (1.11) can also be derived from (1.8) if the same r −2 and ln |r| elastic contributions in [8,26] with only local step interactions are considered, with the local contributions in (1.8) set to be This form of f (r) is consistent with the choice of the surface energy in [26] where elastic interactions among steps are incorporated in the derivation of chemical potentials from a discrete BCF model. For the case α > 0, we introduce the scalings for the governing slope equation (1.11), which after dropping the tildes leads to This rescaling simplifies the discussion of the original problem (1.11).
To model a monotone step train with periodic slope where the maximum and minimum heights in each period do not change in time, we impose periodic boundary conditions on u, u(h) = u(h + H), (1.14) corresponding to Dirichlet boundary conditions on h with a fixed height difference H in each period, Specifically, we investigate the regularity of solutions of (1.13) associated with the periodic boundary condition (1.14) and positive initial data In this work, we obtain a strictly positive lower bound for u, which prevents 1/u from blowing up, and prove the existence and uniqueness of global strong solutions to (1.13) with periodic boundary conditions. We also study the long time behavior of u by investigating two related Lyapunov functionals for (1.13). In the long time, u(t, h) converges to a spatially-uniform solution which only depends on the initial data u 0 . Fig. 1 shows typical evolution of the h-equation (1.10) solved numerically using finite differences with α = 1, L = 1 and H = 2 starting from the monotone initial data h 0 (x) = 0.5 tanh(10(x − 0.5)) + 200(x − 0.5) + 1 on 0 ≤ x ≤ 1. It is observed that the monotone profile of h converges to the line h = 2x as t → ∞, with the corresponding slope profile u approaching to a spatially-uniform solution u = 2 as t → ∞.
Since the model (1.13) for 0 ≤ h ≤ H is equivalent to the model with 0 ≤ h ≤ 1 up to a rescaling, for the rest of this paper we set H = 1.
The structure of the rest of the paper is as follows. In section 2 we introduce two useful Lyapunov functions. In section 3 we show the regularity of solutions of PDE (1.13) and the long time behavior of the PDE solution is further investigated in section 4. Numerical verification of these analytical results are then provided in section 5 with a brief discussion of the effect of the logarithmic term on the transient dynamics of the slope equation.

Two Lyapunov functions.
It is important to note that if u is strictly positive, then (1.13) can be written as (u −1 ) t = (u 3 + u) hhhh , which yields the conservation law, In particular, using the relationship u(t, h) = h x (t, x) and boundary conditions (1.14) and (1.15), we conclude that which gives an underlying connection between the surface height equation (1.10) and the slope equation (1.13). Two Lyapunov energy functions of (1.13) will now be introduced. Specifically, there is a special bi-variational structure embedded in the slope equation (1.13) that is useful in the proof of the regularity of solutions. With the energy density function f in (1.12), the first Lyapunov function is given by Then the variational structure of (1.13) can be written as which leads to the dissipation inequality for F (u), In addition, we define the Lyapunov function and direct calculation gives that δE/δu = (2 + 6u 2 )(u 3 + u) hhhh . Hence the PDE (1.13) actually has a "bi-variational structure" which is given by This is the key point in proving the existence and positivity of u. From (2.6), we also obtain the dissipation inequality for E(u), Moreover, from (2.4) and (2.5), the relation between F (u) and E(u) is This, together with (2.7), gives that and thus where C 0 is a constant depending only on initial data u 0 , and u min is the minimal This formal observation is important for studying the long time behavior of strong solutions to (1.13). From the energy dissipation (2.7), to obtain an exponential decay for energy E, it remains to prove that the dissipation rate of E should be bounded above by −E, i.e., dE dt ≤ −cE, which is usually shown by a logarithmic Sobolev inequality. One example for this is [25]. After proving the global positive lower bound for u, we will use Poincare's inequality to obtain this relation and an exponential decay rate for the energy E. See the analysis for the long time behavior of strong solutions in section 4.
3. Global strong solution. In the following, with standard notations for Sobolev spaces, we denote First we give the definition of a strong solution to PDE (1.13).
Definition 1. For any T > 0, we define a strong solution u(t) to PDE (1.13) to be a positive function that satisfies We now state the main result, the global existence of strong solutions to (1.13) as follows. with initial data u 0 and the following two energy-dissipation equalities hold Further, the lower bound for u is obtained We also point out the following two important facts: If it is assumed that the initial data u 0 is smooth in addition to the conditions in Theorem 1, then with the lower bound in (3.7) we can use standard arguments to obtain higher-order estimates for u, and this therefore ensures positive smooth solutions to PDE (1.13).
in Theorem 1, we will have an asymptotic lower bound for u, First we introduce two lemmas which will be used later.
Lemma 1. For any 1-periodic function v(h), we have the following relation Proof. Notice that Integrating from 0 to 1, we obtain (3.9).

Lemma 2. For any function
Hence we have Proof of Theorem 1. Recall that L = due to Lemma 1. In Step 1, we first test some estimates under the a priori assumption In Step 2, we will show that this a priori assumption leads to a better lower bound. Then we can take the limit using those a priori estimates in Step 1. In Step 3, we prove the energy-dissipation equalities (3.5), (3.6) and thus obtain the existence result to (1.13). In Step 4, we prove that the solution obtained above is unique.
Step 1: A priori estimates. First, we obtain the higher order estimate. Multiplying (1.13) by (u 3 + u) hhhh and integrating by parts leads to (3.13) Then multiplying (1.13) by 3u 2 (u 3 + u) hhhh and integrating by parts yields (3.14) Combining (3.13) and (3.14), we have Then for any T > 0, for any t ∈ [0, T ]. Then by Poincare's inequality we obtain and it remains to estimate u . Using Poincare's inequality and Sobolev embedding again, we have On the other hand, since u has the lower bound assumption (3.12), we have (2.1) and (2.2). Hence This, together with (3.20), yields Thus for any T > 0 we have the lower order estimate (3.24) Third, notice that the function u 3 + u is increasing. From Lemma 2, we have Step 2: Verify the a priori assumption. From (3.25), we have for any t ∈ [0, T ] and for any 0 < β ≤ 1 2 , using the relationship in (2.2) we obtain Hence we have and for 0 < β ≤ 1 2 the right hand side attains its maximum value at Then direct calculation leads to a lower bound u m that depends on E(t) (3.32) Therefore we have and taking the minimal value with respect to t ∈ [0, T ], we have which verifies the a priori assumption (3.12) and shows that u min has a positive lower bound, Thus we obtain (3.7).
Step 3: Proof of energy-dissipation equalities. After obtaining the a priori estimates in Step 1, using standard compactness argument we obtain that for any T > 0 and any  Step 4: Uniqueness of the solution to (1.13). Assume that u, v are two solutions of (1.13). Then we have Combining (3.39) and (3.38), we have Recall that for any p ≥ 0, u p is increasing with respect to u, so from (3.7) there exist constants m, M > 0 whose values depend on u 0 H 2 ([0,1]) , p and L and satisfy Again from (3.7), we have Proof.
Step 1: Exponential decay of the free energy E as t → ∞. By (3.15) we know that E is decreasing with respect to t. Combining (3.15) and the lower bound estimate (3.7), we have = −(2u 2 min + 6u 4 min )E = −βE with the constant β = 2u 2 min + 6u 4 min > 0. Thus E(u(t)) converges to 0 as t → ∞ with an exponential decay rate which shows that (4.1) holds.
Step 2: Showing convergence to the stationary solution u = 1 L . Since we have a positive lower bound for u (3.7), the relations (2.1) and (2.2) hold. Recall the notations in the proof of Theorem 1, u = 1 0 u dh andū = u− u . From Hölder's inequality, This, together with (3.22), gives which, combining with the decay of E in step 1, completes the proof.   Fig. 2 showing convergence to a spatially-uniform solution u = u as t → ∞. (Bottom) A plot showing that u m (t) = min h u(t, h) is bounded below by J (E(t)) given by (3.32) which is in line with the conclusion of Theorem 1, and that the asymptotic lower bound J (E(t)) → 1/(2L) for t → ∞ as in (3.8).
Typical numerical simulations of (1.6) and periodic boundary conditions are plotted in Figure 2, where the spatial variation grows starting from the initial data and leads to a near-rupture transient behavior in the early stage before the solution converges to a constant solution u = u for t → ∞. While for equation (1.13) (or equation (1.11) with α = 1), starting from the same initial data, monotone behavior of the PDE solution u converging to u = u is presented in Fig. 3. With the initial condition (5.1), we have L = 1 0 1 u0 dh = 3.7, and the final state u = 0.27 can be obtained from the formula u = 1/L. The energy-dependent lower bound estimate (3.31) is also plotted in a dashed curve in Fig. 3 (right) in comparison to the numerical results. These numerical studies are conducted using centered finite differences in a Keller box scheme, where the fourth-order PDE (1.13) is decomposed into a system of first-order differential equations  We now revisit the long-time behavior of the solutions of (1.13). It has been shown in the previous section that as t → ∞, the PDE solution u(t, h) approaches to a spatially-uniform solution u = u . To study the long time behavior of the solution, using u(t, h) = u + v(t, h) we obtain the linearized equation for (1.13) as where the function r is defined as r(u) = u 2 (3u 2 + 1). Note that the positivity of solutions to the fourth-order linear PDE (5.3) is not guaranteed. Furthermore, we perturbū by individual Fourier mode disturbances u(t, h) = u + δe ikh e σt + O(δ 2 ), (5.4) where k is the wavenumber and σ represents the rate of the PDE solution converging to u . Substituting (5.4) into (1.13) and linearizing about u = u leads to the dispersion relation 5) which indicates that the steady state solution u ≡ u is stable with respect to any Fourier mode perturbations. As the spatially-uniform solution is approached for long times, the energy E(t) decays exponentially in the form of where C depends on the initial conditions and other system parameters. We note that the energy decay rate 2r(u )k 4 obtained from the linear analysis is inline with the estimated bound of energy decay rate β = 2r(u min ) obtained in Theorem 2. Similar linear analysis for PDE (1.6) leads to the energy decay rate (5.6) where the energy is given by and the dispersion relation σ 0 = −3u 4 k 4 . It is shown in Fig. 4 that the exponentially decay rates of the energy E(t) for the PDE simulations in Fig. 2 and Fig. 3 both agree with the above linearization result as t → ∞. In addition to the evolution of monotone trains of steps with boundary conditions (1.14) and (1.15) described above, we are also inspired by Li and Liu's work on thin film epitaxy [15,16] to further investigate the surface height evolution h in the periodic setting, which is described by the height equation with associated periodic boundary conditions h(t, x) = h(t, x + L). Note that in this regime the solution h is allowed to have both positive and negative slopes  Figure 6. Evolution of the surface height h(t, x) and slope h x (t, x) for equation (6.2) with α = 1 starting from identical initial data used in Fig. 5, showing convergence to a piece-wise constant profile in h and jump in h x .
which is not covered by the argument provided in section 3. While the regularity of solutions of (6.1) in this regime is beyond the scope of this article, we shall present some preliminary numerical results and show the potential connection between (6.1) and thin film epitaxy models with slope selections [15]. Since the logarithmic term in the chemical potential and the mobility function both become singular as |h x | → 0, a direct finite-difference discretization of the equation leads to numerical difficulties. Therefore, we introduce regularization for both the logarithmic term and the mobility function in the model with small coefficients and δ, and numerically study the following regularized equation, (6.2) With = δ = 10 −5 , numerical simulations for (6.2) with both α = 0 and α = 1 starting from identical initial data h 0 (x) = sin(2πx) are presented in Fig. 5 and Fig. 6. It is observed that for the case α = 0, the height profile flattens out with a shock formed in the slope profile. This is consistent with the numerical results for the two-dimensional height profile without line tension in [24]. However, for α = 1, due to the contribution of the logarithmic term in (6.1), a shock is formed in h at the extrema of the height profile, and the solution h approaches to a periodic piecewise linear function with positive and negative slopes. Although the proof provided in Section 3 and Section 4 cannot be applied to this periodic setting directly, we plot the asymptotic bounds ±H/2L for the slope function from (3.8) to compare it with the slope profile obtained from the simulation in Fig. 6, where H = 0.2936 and L = 1, and the step height H is obtained from the difference between the two extrema in the final stage in Fig. 6 (left). It is suggestive that the asymptotic lower bound H/2L < |h x | away from the jump in h x and is still a good estimate of the slope profile.
A number of questions remain to be answered. In the periodic surface height scenario, the interesting dynamics from numerical observation of the regularized problem (6.2) with the non-zero logarithmic contribution need further investigation. While we use classical finite difference methods to numerically simulate the dynamics of the regularized problem (6.2), more sophisticated numerical schemes can be developed to capture the dynamics of the unregularized problem (6.1). For instance, some numerical studies [4,23,24] of the well-known equation (1.4) use a nonlinear Galerkin scheme and finite-element schemes have been applied to the corresponding DL regime [13]. In addition, we are also interested to study whether the logarithmic term in the surface energy contribute to the singularity at the edge of the facets and steps in the modeling of crystal surface evolution. Moreover, without the monotone assumption for initial data, the analysis of the general h-equations (1.4) and (6.1) is very challenging and still open. Due to the form of the mobility function |∇h| −1 , new techniques are needed to rigorously describe the singularity that occurs at |∇h| = 0.
Inspired by Xiang's work in [26] which incorporates nonlocal interactions in the DL regime, we are also interested in the nonlocal contributions in the ADL regime both analytically and numerically. Under the monotone assumption, the h-equation can be written as where H represents the Hilbert transform. While in this work only the onedimensional models are investigated, the equations in higher dimensions are also interesting future subjects. Furthermore, since the general PDE (1.1) can be obtained from linearizing the exponential in the Gibbs-Thomson relation by taking the approximation e δG δh ≈ 1 + δG δh , the corresponding model ∂h ∂t = ∇ · M (∇h)∇e δG δh also needs further investigation.