{'fraction_non_alphanumeric': 0.04662100833402096, 'fraction_numerical': 0.017947025332123112, 'mean_word_length': 4.182488774855677, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling the vertical growth of van der Waals stacked 2D materials using the diffuse domain method

Vertically-stacked monolayers of graphene and other atomically-thin 2D materials have attracted considerable research interest because of their potential in fabricating materials with specifically-designed properties. Chemical vapor deposition has proved to be an efficient and scalable fabrication method. However, a lack of mechanistic understanding has hampered efforts to control the fabrication process beyond empirical trial-and-error approaches. In this paper, we develop a general multiscale Burton-Cabrera-Frank (BCF) type model of the vertical growth of 2D materials to predict the necessary growth conditions for vertical versus in-plane (monolayer) growth of arbitrarily-shaped layers. This extends previous work where we developed such a model assuming the layers were fully-faceted (Ye et al., ACS Nano, 11, 12780-12788, 2017). To solve the model numerically, we reformulate the system using the phase-field/diffuse domain method that enables the equations to be solved in a fixed regular domain. We use a second-order accurate, adaptive finite-difference/nonlinear multigrid algorithm to discretize and solve the discrete system. We investigate the effect of parameters, including the van der Waals interaction energies between the layers, the kinetic attachment rates, the edge-energies and the deposition flux, on layer growth and morphologies. While the conditions that favor vertical growth generally follow an analytic thermodynamic criterion we derived for circular layers, the layer boundaries may develop significant curvature during growth, consistent with experimental observations. Our approach provides a mechanistic framework for controlling and optimizing the growth multilayered 2D materials.


Introduction
Two-dimensional (2D) materials including graphene and transition metal dichalcogenides (TMDs) have garnered unprecedented interest in pursuit of unique electronic, optical, mechanical, and thermal properties [4,5,11,20,8]. Compared to homogeneous monolayers, multilayered heterostructures contain many more degrees of freedom and thus can be ideal platforms for electronic structure engineering of atomically thin 2D semiconducting materials for novel applications. A key challenge in the realization of vertically integrated 2D layers is their synthesis [21,4,7]. Chemical vapor deposition (CVD) has proved to be an efficient and scalable method to grow monolayer 2D materials on a variety of metal substrates [15,23,9,10]. CVD, however, is a complex process that contains many parameters that influence growth. For example, the growth temperature and the deposition flux have been found to be critical parameters for switching from inplane (monolayer) to vertically-stacked multilayer growth. In W S 2 /M oS 2 heterostructures on SiO 2 /Si substrates, high temperatures favor the growth of vertically-stacked multilayers while low temperatures favor monolayer growth [9]. In graphene, a lower deposition flux (e.g., higher concentrations of H 2 in the gas) also tends to favor multilayer growth [26,3]. Determining proper growth parameters is clearly a multivariable problem that until recently was tackled using empirical trial-and-error approaches.
In recent work, we developed a multiscale model of the growth of vertically-stacked 2D materials on a substrate using CVD [24]. The model, which is of Burton-Cabrera-Frank (BCF) type [2], accounts for attachment and diffusion of adatoms, van der Waals (vdW) interactions between the layers and the substrate, and edge energies of the layers. To simplify the system, the layers were assumed to be fully-faceted and so their shapes were constrained to be equilateral polygons (e.g., triangles and hexagons). This work predicted the thermodynamic requirements for growth of vertically-stacked faceted layers. The vdW-BCF model predictions on monolayer vs. multilayer morphologies were validated by comparison with a variety of CVD-synthesized M X 2 (M = Mo, W; X = S, Se, Te) single-species samples grown under conditions of varying temperature and precursor flux.
However, as seen in the experiments in [24,15] and in other references, the layers need not to be faceted and can develop significant, and even negative, curvatures. Because the layer morphologies influence growth and the material properties, it is important to accurately predict the layer shapes as well. In this paper, we extend the vdW-BCF model in [24] to account for arbitrary layer shapes. The resulting system is a highly nonlinear free boundary problem. We analyze the model and derive an analytic thermodynamic criterion for vertical growth assuming the layers are circular. To simulate the model when the layer geometries are unconstrained, we develop a secondorder accurate phase-field/diffusion-domain method (DDM) that enables us to simulate the system by solving a reformulated system (vdW-BCF-DDM equations) in a fixed regular domain.
The diffuse-domain, or smoothed boundary, method is an attractive approach for solving partial differential equations in complex geometries because of its simplicity and flexibility. In this method the complex geometry is embedded into a larger, regular domain. The original PDE is reformulated using a smoothed characteristic function of the complex domain and source terms are introduced to approximate the boundary conditions. An advantage of this approach is that the reformulated equations can be solved by standard numerical techniques without requiring body-fitted meshes, additional interfacial meshes or special stencils and the same solver can be used for any geometry. The diffuse-domain method (DDM) was introduced in [12] to solve diffusion equations with Neumann (no-flux) boundary conditions, to PDEs with Robin and Dirichlet boundary conditions in [14] and to cases in which bulk and surface equations are coupled [22]. Later, in [25] and [17] alternate derivations of diffuse-domain methods for such problems were presented. In [13] a matched asymptotic analysis for general DDMs with Neumann and Robin boundary conditions showed that for certain choices of the source terms, the DDMs were second-order accurate in and in the grid size h in both the L 2 and the L ∞ norms, taking ∝ h, see the recent paper [1] for a rigorous proof.
In [19], a DDM was proposed to solve a BCF model of epitaxial growth of thin, crystalline films that combined a DDM reformulation of the adatom diffusion equations together with a Cahn-Hilliard-type equation to model the dynamics of the films. This approach considered only isotropic edge energies and kinetic coefficients and did not consider van der Waals interactions. Further, the DDM used in [19] did not use a second-order accurate formulation and thus was only first order accurate in (and h assuming ∝ h).
Here, we combine and extend the approaches from [13] and [19] to develop a second-order accurate adaptive finite-difference/nonlinear multigrid method to discretize and solve the vdW-BCF-DDM equations numerically. We investigate the effect of parameters, including vdW interaction energies between the layers, kinetic attachment rates, edge-energies and deposition flux, on layer growth and morphologies. While the conditions that favor vertical growth generally follow the thermodynamic criterion we derived for circular layers, the layer boundaries may develop significant curvature during growth, consistent with experimental observations, that can also influence the growth kinetics.
The outline of the paper is as follows. In Sec. 2, we present and analyze the vdW-BCF model for arbitrary layer shapes. In Sec. 3 we present the phase-field/DDM reformulation of the vdW-BCF model and briefly describe the numerical methods used. In Sec. 4, we present numerical simulation studies and in Sec. 5 we present conclusions and discuss future work. Additional details are provided in the Appendices.

The vdW-BCF model for the growth of vertically-stacked multilayers
Let Ω 0 denote the substrate, Ω 1 denote a layer of atomic height 1 and Ω 2 be a layer of atomic height 2 with boundaries Γ 0 , Γ 1 and Γ 2 , respectively. See the diagrams in Figs. E1 and E3 (left column). The system free energy is taken to be: where E i is the binding energy of layer i that accounts for in-plane bonding and any corresponding vdW interactions. In addition, γ i = γ i (θ i ) is the edge energy of layer i, θ i is the normal angle of layer i (e.g., the angle between the normal vector n Γ i , which points into Ω i−1 , and the x-axis). The function ρ i is the adatom concentration on layer i and ρ ref = Ω s is the concentration of atomic sites (assumed to be the same on the layers). Further, k B is Boltzmann's constant, T is the temperature and the third term in Eq. (1) represents the regular solution model free energy.

Model equations
By requiring mass to be conserved and that the free energy is non-increasing in time, we can derive a thermodynamically-consistent Burton-Cabrera-Frank (BCF)-like system of equations that govern the dynamics of the adatom densities and the layer morphologies and sizes. Here, we only present the nondimensional equations that include several simplifications. A detailed derivation of the equations, a description and justification of the simplifications and the nondimensionalization are given in Appendix A. The nondimensional adatom concentrations satisfy the diffusion equations where D i > 0 is a dimensionless diffusion coefficient, F i is a dimensionless deposition flux and τ −1 d,i a dimensionless desorption rate. These are all assumed to be constant. At the layer boundaries Γ 2 and Γ 1 mass conservation is imposed, which yields the kinetic boundary conditions: Here, q ± i are the diffusion fluxes of adatoms to the layer boundaries, with the " + " and " − " subscripts denoting limits from the i th and i − 1 layers, ρ * is a nondimensional measure of the thermodynamic equilibrium density,γ i = γ i (θ i )+γ i (θ i ), where the primes denote derivatives with respect to θ i , denotes the layer boundary (edge) stiffness, and κ i is the curvature of the edge Γ i (i = 1, 2). The constants k ± i are the dimensionless rates for attachment of adatoms to the edges from the ith (k + i ) and i − 1 (k − i ) layers, respectively. The normal velocity of each layer boundary Γ i is given by where the dimensionless constant β is related to the mobility of an adatom along a curved edge. At the boundary of the substrate, we assume there is no flux of adatoms: ∇ρ 0 · n Γ 0 = 0.

Analysis of vdW-BCF model: Radial solutions and growth criteria
For simplicity, we consider a configuration in which the two layers and substrate are circular and centered at the same point O. We assume that the edge energy and the kinetic coefficients are isotropic. We solve the system (2)-(7) analytically to derive necessary and sufficient conditions for the growth of layer 2. The layers Ω 1 and Ω 2 have radii R 1 (t) and R 2 (t). The substrate has radius R ∞ , which is fixed. We assume that initially 0 < R 2 (0) < R 1 (0) < R 0 and that the dynamics are dominated by diffusion so that the time derivative on the left hand side of Eq. (2) is set to zero (quasi-steady case). We further assume the desorption of adatoms is small and so we set τ −1 d,i = 0. The reduced system can be solved analytically. Here, we present only the results, a complete derivation of the solutions is provided in Appendix B.
The analytical solutions for the densities ρ i are: where A i and B i are given in Appendix B. When the flux of adatoms is only non-zero on the substrate (e.g., F 0 ≥ 0, F 2 = F 1 = 0), which reflects the catalytic decomposition of CH 4 vapor on the substrate surface into mobile radicals (e.g., CH and C) that can attach to the graphene layers [16], the normal velocities of the layer boundaries are given by where γ 1 and γ 2 are isotropic edge energies. The velocities for the more general case with F 1 and F 2 not necessarily equal to zero can be found in Appendix B. Define E 2,1 = E 2 −E 1 to be the binding energy density between the two layers and E 1,0 = E 1 to be the binding energy density between layer 1 and the substrate. The difference between these two energies, is the gain in energy by adding atoms to layer 2 instead of layer 1. An analysis of v 2 in Eq. (10) reveals a sufficient condition for the growth of layer 2: since the denominator in Eq. (10) is always non-positive. This is analogous to the growth criterion derived in [24] for faceted layers. This condition states that the difference between the binding energies, ∆E, must be large enough to overcome the energy penalty of increasing the layer perimeter. It follows that if R 2 > R 2,c := γ 2 /∆E, then layer 2 always grows, regardless of the size of layer 1. This is analogous to a critical nucleation size. Further, if R 1 > R 1,c := γ 1 γ 2 /R 2 − ∆E then layer 2 always shrinks. When R 2 is close to R 2,c , layer 2 may grow due to kinetic effects. That is, R 2 may surpass R 2,c before R 1 surpasses R 1,c . Whether this occurs depends on the values of the parameters. For example, slowing down the growth of the first layer (e.g., by decreasing F 0 ) or increasing the rate of growth of the second layer (e.g., by increasing D 1 , k − 2 or k + 1 ) increases the region of kinetically-driven growth. We call R 2,k the kinetic critical radius-that is, if R 2,c > R 2 (0) > R 2,k , then the second layer grows due to the kinetics of the system.
By solving for the radii R 1 and R 2 numerically and varying the initial radii, we can estimate R 2,k numerically and construct a phase diagram for the growth of the 2nd layer. As an example, we fix the parameters ∆E = 0.05, F 0 = 0.1, γ 1 = γ 2 = 0.01, ρ * = 0.5, D 1 = 1, k − 2 = k + 1 = 0.5 and R ∞ = 3.8. We then vary the initial sizes of the layers R 1 (0) and R 2 (0), keeping R 2 (0) ≥ R 1 (0). The resulting phase diagram is shown in Fig. E2(a). Also observe that for R 2 in between R 1,c and R 2,k , the 2nd layer grows transiently before shrinking to zero size. Example trajectories of the layer dynamics are shown in Fig. E2(b).

Reformulation of the vdW-BCF model of multilayer growth using the diffuse domain method
To solve the vdW-BCF equations for unconstrained layer geometries, we reformulate the system using the diffuse domain method (DDM). Here, we combine and extend the approaches from [13] and [19] to develop a fully-second order accurate DDM for the vdW-BCF system. We embed the substrate and layer domains into a larger, rectangular domainΩ and we introduce a diffuse domain function ϕ to mark the locations of the layers and substrate (e.g., approximate atomic height). In particular, ϕ ≈ 0 in the substrate (Ω 0 ), ϕ ≈ 1 in layer 1 (Ω 1 ) and ϕ ≈ 2 in layer 2 (Ω 2 ).
In order to facilitate comparisons with theory from the previous section, we assume that the outer boundary of the substrate is circular and so we introduce another diffuse domain function ϕ ∞ to identify the deposition domain Ω = Ω 0 ∪ Ω 1 ∪ Ω 2 , where ϕ ∞ ≈ 1, within the larger domainΩ. See Fig. E3(a). The diffuse domain variables change rapidly but smoothly across the boundaries (e.g., steps) as shown in Fig. E3(b). The width of these narrow transition layers is ≈ , a small parameter. The boundaries of the substrate and layers 1 and 2 correspond to ϕ ≈ 0.5 and ϕ ≈ 1.5, respectively.
The kinetic boundary conditions are incorporated via source terms and the dynamics of the layers are captured by evolving the diffuse domain function ϕ. In addition, we follow [19] and solve only two adatom diffusion equations in the extended domainΩ. A brief description of the derivation and an asymptotic analysis of the vdW-BCF-DDM, which demonstrates that the vdW-BCF-DDM system approximates the sharp interface vdW-BCF model to O( 2 ), are given in Appendix C. Here, we present only the resulting equations: where the kinetic boundary conditions (3) and H 1 is the approximate characteristic function of layer 1: Further, the flux F 0 (ϕ) corresponds to the flux on the substrate and D 0 (ϕ) corresponds to the adatom diffusion coefficients on the substrate and layer 2 Analogously, the extended vdW energies and kinetic attachment rates are defined as and The evolution of the layers is implicitly captured by evolving ϕ: where the right hand side of Eq. (22) models the normal velocity from Eq. (7). Note that since the outer boundary of the substrate does not change we do not need to pose an evolution equation for ϕ 0 . In Eqs. (22) and (23), G(ϕ) = 2B(ϕ) is an extended double well potential: As shown in Appendix C, and confirmed by our numerical results in the next section, the vdW-BCF-DDM system is second order accurate with respect to the interface thickness . Moreover, our diffuse interface model can be extended to simulate the more nonlinear model derived in Appendix A and to simulate an arbitrary number of vertically-stacked layers (see Appendix A.5). Finally, at the boundary of the larger domain ∂Ω, we take the conditions ∇ρ 0 · n = ∇ρ 1 · n = ∇ϕ · n = ∇µ · n = 0.
The model is insensitive, however, to the choice of boundary conditions on ∂Ω.

Numerical Results
To solve the vdW-BCF-DDM system (13)-(25) numerically, we develop a massconservative, semi-implicit, second-order accurate, adaptive finite-difference method using Crank-Nicholson discretization in time and centered differences in space, by extending our previous work, e.g. [6]. To solve the nonlinear discrete system at the implicit time level, we use a full approximation storage (FAS) nonlinear multigrid method. Block-structured adaptive mesh refinement is utilized to efficiently discretize the system. The details of the method are provided in Appendix D.
We begin by considering the isotropic, quasi-steady case so we may compare our numerical results to the analytical solutions presented in Sec. 2.2 to validate the accuracy of our approximations. We then consider time-dependent diffusion and anisotropic edge energies and kinetic coefficients. We perform parametric studies to determine the effect of parameters on the growth and morphologies of the layers.

Quasi-stationary dynamics
We consider the same set up as in Sec. 2.2. Initially, two layers are centered at the origin with different radii R 1 and R 2 and the edge energy and kinetic coefficients are isotropic. The two islands are bounded by a larger circular substrate with radius R 0 . The initial condition for the diffuse domain variable is such that ϕ ≈ 1 approximates layer 1 and ϕ ≈ 2 approximates layer 2. We take which corresponds to the region containing the substrate and the two layers where deposition and growth take place. The parameter is the thickness of the layer and substrate boundaries. The initial radii of the layers are R 1 (0) = 1.2 and R 2 (0) = 0.6. The outer radius of the substrate is R ∞ = 1.8. The physical parameters are taken to be The computations are carried out on a square domain [−2, 2]×[−2, 2]. A 4-level adaptive mesh is employed, which consists of a root level with mesh size h 0 and three refinement levels above it so that the finest mesh size h 3 = h 0 /8. In order to test the convergence rate corresponding to different values of , we refine the root level grid size h 0 and together, and hence all the finer level grid sizes h 1 , h 2 and h 3 are refined as well. In particular, we set h 3 = 6.4 . The mesh is refined according to values of |∇ϕ| + |∇ϕ 0 | over the entire domain (see Appendix D). The time step is taken to be ∆t = 0.8 × 10 −4 to ensure that the time errors are small compared to spatial errors; the method is stable (and accurate) for larger time steps.
Five different values of are used for the convergence test, namely, 1 = 0.8, where k = 0 denotes the substrate and k = 1, 2 denote the layers. The convergence rate is obtained by  Fig. E4(a). We can observe that the numerical results approach the analytical solution as decreases. The corresponding errors and rates of convergence are presented in Table 1, which indicates that the numerical method converges to the analytic solution with an overall second order convergence rate in both the 2 and ∞ norms, as predicted by the asymptotic analysis in Appendix C.

Fully time-dependent case
Next, we include the time derivatives in the adatom diffusion equations. The physical parameters, the computational domain and the numerical parameters are the same as in Table 1. Convergence test for the adatom concentrations ρ 2 , ρ 1 and ρ 0 under quasisteady dynamics from §4.1.
the previous section. Since we do not have an analytic solution in this case, we compare the results obtained using different (and hence h 0 ) with each other. The horizontal slices of the adatom concentrations ρ 2 , ρ 1 and ρ 0 are shown at time t = 0.1 in Fig.  E4(b). Compared to the quasi-steady case, the adatom concentrations in each layer are smaller and there is less variation across the layers. Correspondingly, the layers do not move as rapidly in the time-dependent case with the first layer growing more slowly than the second, compared to the quasi-steady case ( Fig. E4(c)). Fig. E4(b) also shows that the results converge as is decreased. To estimate the accuracy and quantify the rate of convergence, we define the consecutive errors as where the ϕ j , with j = 0, 1 and 2 are the approximate characteristic functions on the substrate, layer 1 and layer 2 respectively. They are defined as and these functions are evaluated at = i . The errors and rates of convergence, which are calculated from the consecutive errors at time t = 0.1 in an analogous way as in the previous section, are presented in Table 2. As in the quasi-steady case, we observe that the results converge with second order accuracy in both the 2 and ∞ norms. t=0.10

Anisotropic dynamics
We now consider the case in which the edge energies and kinetic coefficients are anisotropic: where is the kinetic coefficient anisotropy function, θ is the normal angle (e.g., angle between the normal vector and the x-axis), and θ 0 is a reference angle which is taken to be θ 0 = π/n. The edge energies are defined analogously: where ξ s is the edge energy anisotropy function. The coefficients k,n and s,n measure the anisotropy strengths. In this paper, we only consider 3-fold (n = 3) and 6-fold (n = 6) anisotropies, which reflect the symmetries of M oS 2 and graphene multilayers, respectively. The trigonometric functions are calculated using ϕ: , and sin(θ) = ϕ y where we introduce a small parameter δ = 10 −6 to avoid singularities. Then, cos(3θ) and cos(6θ) can be calculated using the trigonometric identites: We next consider the quasi-steady dynamics of two anisotropic layers. Initially, the layers are taken to be circular with radii R 1 (0) = 1.0 and R 2 (0) = 0.2. The outer boundary of the substrate is R ∞ = 3.8. The physical parameters are taken to be: Note that unlike the previous examples, the only non-zero flux is on the substrate F 0 , which as discussed before reflects the assumption that the reactions to produce the attaching species occur only on the substrate surface [16]. Note that ∆E = 0.4 and γ 2 /R 2 (0) − γ 1 /R 1 (0) = 0.04 so that growth would occur under isotropic, quasi-steady dynamics (recall the growth condition in Eq. (12)). The parameters for the anisotropy are set as n = 6, k,n = 0.3, s,n = 0.01, The morphologies of the growing layers are shown in Fig. E5(a). In both the 6-fold and 3-fold anisotropic cases, layers 1 and 2 grow. In the 6-fold case, the layers are nearly faceted at early times while the corners are smoothed slightly from the surface diffusion. At later times, both layers develop negative curvature. In the 3-fold case, layer 1 evolves to a convex triangular shape at early times while layer 2 develops negative curvature early on. At later times, the corners of layer 1 somewhat elongate with their curvature being set by the surface diffusion coefficient (see Fig. E9(c)). The corresponding adatom concentrations are shown in Fig. E5(b) where we see the adatoms diffusing toward both layers driving their growth. In Fig. E5(c), the adaptive mesh is shown for the 6-fold anisotropic case. Observe that there is a fine mesh near the outer boundary of the substrate, which does not change. The mesh near the boundaries of layers 1 and 2 is dynamically refined and and the mesh in the bulk regions is coarsened. In the anisotropic case, we also observe second-order accurate convergence in 2 and ∞ , see Appendix E.

Parameter studies
We next investigate the effects of the physical parameters on the growth of the layers.
In particular, we consider the binding energy differences ∆E, the edge energy γ and the surface diffusion β, flux F 0 and the kinetic attachment rates k − 2 and k + 1 . We fix all the other parameters as in Eq.(42) and describe only those parameters that are changed.
Binding energy differences. We first investigate the effects of ∆E on the growth rate of layer 2. The morphologies and adatom concentrations for 6-fold anisotropic layers obtained from the quasi-steady dynamics are shown in Figs. E6(a) and (b), respectively. Consistent with theory (Sec. 2.2), the vertical growth of layer 2 is only preferable when ∆E > 0.04, based on Eqs. (12) and (42), and that growth rate increases with ∆E.
Further, the growth of layer 2 occurs at the expense of that of layer 1; the size of layer 1 is a decreasing function of ∆E. In all the cases, layer 1 is nearly faceted at early times, and develops negative curvatures at late times as layer 1 increases in size. Similar morphologies are observed for layer 2 with negative curvatures occurring when layer 2 is large enough (e.g., ∆E = 0.8).
For comparison, the morphologies for 6-fold anisotropy obtained from the fully time-dependent dynamics are shown in E7. Compared to the quasi-steady case, we observe that the growth of layer 1 is significantly slower but that layer 2 actually grows more rapidly. Further, layer 2 grows even when ∆E = 0. This reflects the fact that vertical growth is more favorable when the growth rate of layer 1 is decreased, which is suggested by the theory in Sec. 2.2.
In Figs. E8(a) and (b), the morphologies and adatom concentrations are shown, respectively, for 3-fold anisotropic layers using the fully time-dependent dynamics. Qualitatively, the results are similar to the 6-fold case in Fig. E7 although we observe that negative curvature occurs first in layer 2 before being manifest in layer 1.
Edge energy, surface diffusion and flux. In Fig. E9(a), we show the effects of edge energy γ on the growth of the layers in the fully time-dependent case. In both 6-fold and 3-fold anisotropies, we see that the growth rate of layer 2 decreases as we increase γ, and the layer 2 even shrinks when γ is large enough (γ = 0.16 or larger). The size of layer 1 is also decreased and the layer morphologies are smoother and the negative curvature on the layers disappears as γ is increased.
As seen in Fig. E9(b), surface diffusion also decreases the sizes of layer 2 and smoothens the layer corners although the negative curvature of the layers remains. In the 6-fold anisotropic case, layer 1 is also decreased in size as β increases while in the 3-fold anisotropic case, layer 1 is actually a little larger due to the decreased curvature at the vertices.
Next, we examine the effects of the adatom flux F 0 on the layer dynamics. As shown in Fig. E9(c), decreasing the supply of adatoms on the substrate (F 0 ) benefits the growth of second layer, which agrees with reported experimental observations for vertical growth of 2D materials (e.g., [24]). Moreover, in the case of 6-fold anisotropy, we see that both layers develop negative curvatures at small F 0 , but as F 0 is increased the shapes become more facetted. Similar features are observed in the 3-fold anisotropic case, except when F 0 = 100, where kinks with negative curvature develop at the boundary of layer 1. This feature persists under mesh refinement and seems to be associated with deposition only occurring on the substrate. If adatoms are deposited on all the layers, then layer 1 is convex at an equivalent size.
Kinetic coefficients. In Fig. E10(a), the kinetic parameter k − 2 is varied from 0.5 to 4.0 for layers with 6-fold anisotropies. As predicted by the theory in Sec. 2.2, increasing k − 2 favors the growth of layer 2 at the expense of layer 1. Both layers acquire negative curvature as they grow. In Fig. E10(b), we take k + 1 = k − 1 and vary this value from 0.5 to 4.0. In this case, the growth of layer 2 is insensitive to these changes, which is surprising because theory suggests that increasing k + 1 increases layer 2 growth (Eq. (10)). The reason for the discrepancy is that a morphological instability occurs on layer 1 that accelerates its growth relative to that of layer 2. Because layer 1 grows faster, this reduces the number of adatoms available for layer 2 growth.
The growth of 3-fold anisotropic layers subject to the same changes in the kinetic parameters shows somewhat different results. As seen in Figs. E11(a) and (b), increasing k − 2 and k + 1 both favor the growth of layer 2. Further, when k − 2 is increased, only layer 2 acquires negative curvature while layer 1 remains convex, in contrast to the results found for 6-fold anisotropy. In addition, when k + 1 is increased, the morphological instability of layer 1 found in the 6-fold case is not present in the 3-fold case. Because of this layer 1 in the 3-fold case does not grow as rapidly, relative to that of layer 2, which enables more adatoms to be available to drive the growth of layer 2.

Conclusions
Epitaxial growth of 2D materials is a complex process, influenced by thermodynamic, kinetic and growth parameters, often leading to diverse and complex growth morphologies determined both by atomic-scale phenomena and by the elastic interactions of surface features and defects and transport of diffusing molecules over length scales of hundreds of nanometers. No single model can describe all the processes involved. In this paper, we derived a general continuum vdW-BCF model to describe the growth of vertically-stacked, arbitrarily-shaped multilayered 2D materials. The model accounted for (i) energy changes upon incorporation of adatoms into the growing 2D layers, (ii) kinetic barriers to attachment, (iii) distinct vdW interactions between the 2D layers and the substrate, (iv) energy penalties associated with the layer edges, and (v) the entropy of the adatoms. This is an extension of our previous work where we developed and analyzed an analogous model for faceted layers where the layer dynamics was much simpler [24]. The vdW-BCF system presented here represents a highly nonlinear free boundary problem.
We analyzed a nondimensional version of the vdW-BCF model and derived an analytic thermodynamic criterion for vertical growth of stacked 2D materials assuming the layers are circular. To solve the system numerically, we used a second-order accurate phase-field/diffusion-domain method (DDM) that enabled us to solve the dynamic equations in a fixed regular domain. To discretize and solve the vdW-BCF-DDM reformulated system, we developed a second-order accurate finite-difference/nonlinear multigrid method using adaptive, block-structured Cartesian mesh refinement. We demonstrated convergence of the numerical methods and investigated the effect of parameters on the layer growth and morphological evolution. While the conditions that favor vertical growth generally follow the thermodynamic criterion we derived for circular layers, the layer boundaries may develop significant curvature during growth and even morphological instabilities. These deviations from faceted shapes can alter the growth dynamics of the layers and can hinder or enhance vertical growth.
Experiments show a wide variety of layer morphologies, including layers with negative curvature, which our model is capable of reproducing. A small sample of experimental layer morphologies are shown in Fig. E12 together with our numerical simulations. Fig. E12(a) shows a SEM image of bilayer graphene from [15] (left) that exhibits a star-shaped layer 1 and a nearly circular layer 2. The image on the right in Fig. E12(a) is a numerical simulation at time t = 4 with the parameters from Eq. (42) except that k − 1 = 10, F 0 = 1 and E 2 = 1.0. Fig. E12(b) shows a SEM image of bilayer graphene with a twisted layer 2 from [15] (left). This experiment was motivated by the observation that electronic structure of bilayer graphene can be altered by changing the relative twist angle, yielding a new class of low-dimensional carbon systems. To simulate twisted bilayer graphene, we modify the reference angle θ 0 of the kinetic coefficient ξ k (θ) in Eq. (36). In particular, we set Although we performed our study using a range of nondimensional parameters, atomistic and mesoscale models can be used to provide specific material parameters. For example, DFT simulations can provide estimates for vdW interaction energies as well as edge energies and kinetic barriers for attachment [3,18,24]. Incorporating such parameter estimates will be explored in future work.
Further, in this paper we have focused on single material homostructures due to perfect lattice matching and hence there are no interior strains. In the TMD family, one can go further and consider M oX 2 /W X 2 heterostructures (M = Mo, W; X = S, Se, Te) without introducing lattice mismatch. However, taking full advantage of the device properties accessible through marriage of disparate 2D materials requires understanding the role of strain in the competition between vertical and in-plane lateral growth. We expect that strain-driven defect formation and stacking-site symmetry breaking will significantly modify the potential energy surface, affecting the thermodynamics of monolayer vs. multilayer morphologies and the kinetics of adatom attachment. Such effects will also be considered in future work. vertically-stacked multilayer growth Appendix A. 1

. Mass Conservation
We define the total mass to be: where Ω s,i are the concentrations of atomic sites in the layers (i = 1, 2). Then, mass conservation requires where F i is the deposition flux on layer i and τ −1 d,i are desorption rates. Combining these two equations and using the Reynolds transport theorem gives: where v i are the velocities of the adatoms on the layers and substrate. For simplicity, we assume that Ω s,i = Ω s . We also assume that the boundary of the substrate Γ 0 does not move. Therefore, combining Eqs. (A.2) and (A.3) and using the divergence theorem we obtain where v i = v · n Γ i is the normal velocity of layer i, and ρ + i = ρ i | Γ i , ρ − i = ρ i−1 | Γ i are the boundary conditions for the densities at the ith layer from the step up and down respectively. Next, assuming that then the last term in Eq. (A.4) can be written as where J ± i (for i = 1, 2) denote the fluxes at the ith layer from a step up and down, respectively, with and we have assumed that there is no flux at the substrate boundary: J 0 ·n ∂Γ 0 = J 0 = 0. Further, the boundary conditions for Eq. (A.5) on Γ i are taken to be In order to satisfy mass conservation, we then have where ∂ s denotes the arclength derivative and J i represents surface fluxes (e.g., arising from the diffusion of adatoms along the layer edges). To obtain constitutive laws for the fluxes q i , J i and J i , we require that the system dissipates the free energy when the deposition flux F i = 0 and desorption coefficient τ −1 d,i = 0.

Appendix A.2. Free Energy Dissipation
Taking the time derivative of the free energy E from Eq. (1) and using the Reynolds transport theorem, we obtain whereγ i = γ i (θ) + γ i (θ) and the primes denote derivatives with respect to θ, the normal angle (e.g., angle that the normal vector makes with the x-axis). Defining the free energy density f and the chemical potential µ to be and applying the divergence theorem, we obtain where µ i = µ(ρ i ). Next, using Eq. (A.5) in Eq. (A.17) we obtain dE dt = Integrating by parts and using the divergence theorem, we obtain dE dt = where we have defined µ ± i = µ(ρ ± i ). See the previous subsection for the definitions of where we have integrated by parts on the edges Γ 1 and Γ 2 and defined Hence, to have energy dissipation (in the absence of flux and desorption), we may take the constitutive relations for the fluxes: where the β i are related to the mobility of an edge atom along a curved step, and the (linear) kinetic boundary conditions: where k ± i are kinetic attachment coefficients.

Appendix A.3. Model simplification
Since Further, µ ± i can be approximated as We further neglect the effects of anisotropy in the surface fluxes (e.g., we assume that the edge energy anisotropy is smallγ i (θ) ≈γ i ), although we keep the effects of anisotropy in the kinetic coefficients and in ρ BC i . Surface diffusion anisotropy will be considered in future work. It follows that the diffusional and surface fluxes can be approximated by Ωs the velocities can be approximated as and the kinetic boundary conditions can be approximated as s . Finally, Eq. (A.5) can be approximated by Let ρ ref = Ω s , L be the characteristic size of layer 1 and take the time scale to be T = L 2 Ω 2 s / (4Dk B T ), where D is a characteristic diffusion constant. Define the nondimensional density ρ i = ρ i /Ω s and the nondimensional flux F i = T F i /Ω s , where the nondimensional desorption coefficient is τ d,i = τ d,i /T . Then, the nondimensional adatom density equation (A.39) becomes: where D i = D i /D is the nondimensional diffusion coefficient. The kinetic boundary conditions become and E is a characteristic value of the binding energies. Finally, the nondimensional velocities are: where β i = β iγi 4Dk B T L 2 are nondimensional edge diffusion coefficients. Dropping the primes, this is the system given in Sec. 2.1.

Appendix A.5. The vdW-BCF model equations for an arbitrary number of vertically-stacked layers
One can extend the vdW-BCF model derived in the previous sections to describe the dynamics of an arbitrary number of layers. The resulting (nondimensional) system is where n is the number of layers. The boundary conditions at the boundary of the first layer with the substrate, Γ 1 , are given as and for all the layer boundaries (e.g., steps) Γ i (for i = 2, . . . , n) are: where γ i denotes the step stiffness and κ i is the curvature of the ith step Γ i , for i = 1, 2 . . . n. The normal velocity of each step Γ i is given by

Appendix B. Details of the derivation of radial solutions to the vdW-BCF model
We now derive the analytic solutions ρ i (r, t) in the quasi-steady state limit. That is, we drop the time derivatives in the adatom diffusion equations. We first rewrite Eq. (2) as where we have also neglected desorption and taken τ −1 d,i = 0. Integrating twice we obtain: where A i and B i are unknown constants. For r < R 2 (t), the solution in Eq. (B.2) satisfies the following boundary conditions: We then obtain For R 2 < r < R 1 , the solution in Eq. (B.2) satisfies At r = R 2 , we obtain At r = R 1 , we obtain such that such that Summarizing, we obtain the analytic solution The corresponding velocities of the layer boundaries are

Appendix C. The diffuse domain method: Details and asymptotic analysis
For simplicity, consider the problem with a single layer: where i = 0, 1 denote the substrate and layer, respectively. The kinetic boundary conditions are: with the normal velocity of Γ 1 (t) = ∂Ω 1 (t) given by In the above, κ is the curvature of Γ 1 . Next, following [14,13], we can reformulate Eqs. (C.1)-(C.3) as where ϕ = ϕ(x, t) is a phase-field function that approximates the characteristic function of Ω 1 (t), ϕ c = 1 − ϕ approximates the characteristic function of the substrate Ω 0 , and µ = B (ϕ) − 2 ∆ϕ is the chemical potential where B(ϕ) = 18ϕ 2 (1 − ϕ) 2 is a double well free energy. Eqs. (C.5) and (C.6) are solved in a large rectangular domainΩ that contains Ω 1 and Ω 2 . For simplicity, we do not include ϕ ∞ to specify that the deposition domain on the substrate is a circle and we assume that the kinetic parameters and edge energies are isotropic. The evolution of the layer is captured by the Cahn-Hilliard-like model: Below, we demonstrate using the method of matched asymptotic expansions that the DDM (C.5)-(C.11) yields a second-order accurate approximation of the sharp interface system (C.1)-(C.4). The analysis can easily be extended to the more complete model presented in the main text in Sec. 3 where two layers are considered and the substrate geometry is circular (implemented via ϕ ∞ ).
Matched asymptotic expansions. Away from the layer 1 boundary Γ 1 (t), we assume that all variables are smooth and have regular expansions in , e.g., while away from Γ 1 , ϕ = 1 inside Ω 1 and ϕ = 0 outside Ω 1 to all orders. Accordingly, we see that ρ (0) i satisfies Eq. (C.1), while the first order perturbations satisfy: To provide the boundary conditions for the diffusion equations, we need to analyze the behavior of the system near Γ i . To argue that ρ i is a second order approximation to the sharp interface solution ρ i , we need to demonstrate that ρ (1) i = 0. Near Γ i , we introduce a stretched, local coordinate system: x(s, t; ) = X(s, t) + zn(s, t) (C.14) where X(s, t) is a parameterization of Γ 1 (t), s is arclength, n(x, t) is the normal vector that points out of Ω 1 , z = r(x, t)/ is a stretched normal coordinate and r(x, t) is the signed distance from x to Γ 1 (t). In the local coordinate system, derivatives become: where the time derivative on the left hand side of Eq. (C.17) is the full time derivative and the time derivative on the right hand side is the time partial derivative in the inner variables, and v 1 is the effective diffuse interface normal velocity of Γ 1 . Note that n = −∇ϕ/|∇ϕ|. We assume that near Γ 1 (t), the inner variables can be expressed as, ρ i (z, s, t) = ρ i (X(s, t) + zn(s, t), t).
(C. 18) We assume that in the inner expansion, all variables have a regular expansion in the stretched coordinates, e.g., To match the inner and outer expansions, we assume that there is a region of overlap where both expansions are valid and must match. In particular, if we evaluate the outer solution in the inner variables, this must match the limits of the inner solutions away from the interface. That is, as z → ±∞ and → 0 with z → 0 ± . Using the inner and outer expansions and equating the powers of , we obtain i (s, t) + zn · ∇ρ i (s, t) + zn · ∇ρ , t), t).
Next, transforming the equations, plugging in the inner expansions and equating powers of we derive equations governing the inner solutions. At leading order O( −2 ), we obtain From these equations (and the matching conditions), we conclude that Integrating these equations from −∞ to ∞ in z, using that v (0) 1 is independent of z, ϕ (0) (+∞) = 0 andφ (0) (+∞) = 1, we obtain where we have additionally used Eq. (C.26) and assumed that ∂ zĝ (0) = 0, a fact that will be justified later. From the matching conditions Eqs. (C.21) and (C.22), we obtain where, as stated earlier, ρ (0) i are the limiting values of the leading order outer solution on Ω i and we have defined g (0) =ĝ (0) . Now, using that g (0) = ρ * (−E 1 −γ 1 κ), another fact we will demonstrate later, then we recover the kinetic boundary conditions Eqs. To justify the assumptions forĝ (0) and to determine the normal velocity v 1 , we need to analyze the Cahn-Hilliard-like system (C.8) and (C.9). Before doing this, however, we proceed to the next order in the inner expansion for the adatom diffusion equations in order to determine the boundary conditions for Eq. (C.13) for the outer solution at the next order, ρ where we have used which follow from Eqs. (C.27) and (C.28) and using the matching conditions. Next, we observe that on Γ i : where we have used that ρ This motivates us to rewrite Eq. (C.33) as where we have used that ∂ t ρ (0) 1 . Next, after a series of calculations, we rewrite Eq. (C.40) as where we have also assumed that ∂ tφ (0) = 0, which will be shown later. Integrating Eq. (C.41) in z from −∞ to +∞, using the matching conditions and that ∂ zρ 1 − g (1) .
Here, we obtain where we have used that G(0) = G(1) = G (0) = G (1) = 0 and To make further progress, we observe that Using these in Eq. (C.60), together with the matching conditions, we obtain: Combining this with Eq. (C.65) we obtain Next, from Eqs. (C.35) and (C.36), and the matching conditions, we havê We can therefore conclude that ρ which demonstrates that the DDM (C.5)-(C.9) provides a 2nd order accurate approximation in to the sharp interface model.

Appendix D. Details of the numerical method and implementation
Numerical method We use the Crank-Nicolson scheme to discretize the fully timedependent system Eqs. (13)- (25) in time on larger square domainΩ. In particular, we let δt > 0 denote the time step, and assume that ρ n 0 , ρ n 1 , ρ n 0 and µ n are the solutions at time t = nδt. We then find the solutions at time t = (n + 1)δt: ρ n+1 0 , ρ n+1 1 , ρ n+1 0 and µ n+1 by solving with the following boundary conditions ∇ρ n+1 0 · n = ∇ρ n+1 1 · n = ∇ϕ n+1 · n = ∇µ n+1 · n = 0 on ∂Ω. Implementation Standard, cell-centered central-difference finite difference methods are used, together with a block-structured adaptive mesh, to discretize the equations in space. The nonlinear equations at the implicit time level are solved using an efficient nonlinear FAS multigrid solver. See [6] for details. Here, we use a 4-level blockstructured adaptive mesh, which consists of one root level (grid size h 0 ) and three refinement levels (grid size h i ) with refinement ratio of 2. For each adaptive mesh level, we refine the grid cell (i, j) wherever h i |∇ϕ i,j | > q tol . Here, we set q tol = 0.01.

Appendix E. Convergence of anisotropic layer dynamics
Here we present the convergence analysis using the fully time-dependent dynamics.
The results for quasi-steady dynamics are similar (not shown). Using the parameters in Sec. 4.3, we analyze the convergence of our schemes at time t = 0.1. The consecutive errors (e.g., Eq. (30)) and convergence rates for the adatom concentrations are given in Tables E1 and E2 for 6-fold and 3-fold symmetric anisotropic edge energies and kinetic coefficients, respectively. The results suggest the scheme is second-order convergent in both the 2 and ∞ norms. Table E1. Convergence test for adatom concentrations ρ 2 , ρ 1 and ρ 0 using 6-fold symmetric anisotropic edge energies and kinetic coefficients in §4.3.   Table E2. Convergence test for adatom concentrations ρ 2 , ρ 1 and ρ 0 using 3-fold symmetric anisotropic edge energies and kinetic coefficients in §4.3.    Figure E2. (a). Morphology diagram, assuming that the layers are circular, showing the dependence of layer 2 growth on the size of layers 1 and 2. In particular, the sign of v 2 is shown for different sizes of the layers (R 1 , R 2 ). When R 2 > R 2,c layer 2 always grows. When R 1 > R 1,c layer 2 always shrinks. When R 1,c < R 2 < R 2,k layer 2 grows transiently before shrinking. When R 2 > R 2,k layer 2 grows because R 2 increases past R 2,c sooner than R 1 crosses R 1,c . See text for details on R 1,c , R 2,k and R 2,c . (b). Sample trajectories of the layer radii R 1 and R 2 in time, starting from different initial radii. The parameters are as in Eq. (42) except with E 2 = 1.05 so that ∆E = 0.05.        Figure E10. The effects of kinetic attachment coefficients on the dynamics of 6-fold anisotropic layers 1 and 2. See text for parameters. In (a), only k − 2 is varied. In (b), k + 1 = k − 1 are varied together. The kinetic parameter k − 2 promotes layer 2 growth, as predicted by theory. The growth of layer 2 is insensitive to simultaneous changes in k + 1 and k − 1 although layer 1, however, is significantly affected and undergoes a morphological instability. Figure E11. The effects of kinetic attachment coefficients on the dynamics of 3-fold anisotropic layers 1 and 2. See text for parameters. In (a), only k − 2 is varied. In (b), k + 1 = k − 1 are varied together. In contrast to the 6-fold case shown in Fig. E10, both kinetic parameters k − 2 and k + 1 = k − 1 promote the growth of layer 2 in the 3-fold case. Further, the morphological instability observed in the 6-fold case is not present in the 3-fold case.