Computational synthesis of 2D materials grown by chemical vapor deposition

The exotic properties of 2D materials made them ideal candidates for applications in quantum computing, flexible electronics, and energy technologies. A major barrier to their adaptation for industrial applications is their controllable and reproducible growth at a large scale. A significant effort has been devoted to the chemical vapor deposition (CVD) growth of wafer-scale highly crystalline monolayer materials through exhaustive trial-and-error experimentations. However, major challenges remain as the final morphology and growth quality of the 2D materials may significantly change upon subtle variation in growth conditions. Here, we introduced a multiscale/multiphysics model based on coupling continuum fluid mechanics and phase-field models for CVD growth of 2D materials. It connects the macroscale experimentally controllable parameters, such as inlet velocity and temperature, and mesoscale growth parameters such as surface diffusion and deposition rates, to morphology of the as-grown 2D materials. We considered WSe2 as our model material and established a relationship between the macroscale growth parameters and the growth coverage. Our model can guide the CVD growth of monolayer materials and paves the way to their synthesis-by-design.


Introduction
Two-dimensional (2D), also called single (or mono-) layer, materials are crystalline materials of only a single layer of atom thick. Their unique quantum properties are critical for technological applications in topographical computing, quantum informatics, and optoelectronic devices that have recently attracted Invited Feature Paper both size-and shape dependence as well [1]. Thus, controlling the growth conditions and kinetics can further customize their properties. For example, we may refer to altering the conductivity and ferromagnetic response of MoS 2 by shifting from a zigzag to an armchair structure in nanoribbons [2].
CVD method-and its variants-is a versatile technique that can control the chemistry and morphology of 2D materials. It has several advantages with respect to other synthesis methods, such as mechanical and chemical exfoliation techniques, including high-quality synthesis of monolayer materials with minimum contamination. However, it involves coupled physical processes with different length and time scales that hinder its understanding. Furthermore, due to the complexity of the material growth problem, a subtle change in any of the control parameters can obstruct the synthesis process. For example, a slight variation of the Mo:S ratio along the CVD reactor leads to different growth morphologies [3]. Also, a change in the growth direction from in-plane to normal to the substrate's plane upon changing the substrate's location has been reported triggered by the gradient of the precursor concentration [4,5]. Different 2D materials have already been synthesized using the CVD technique by trial-and-error experimentation. However, the main obstacle in their industrial application is the controllable and reproducible synthesis.
Experimental setups are commonly expensive, inflexible, and lack the versatility for changing and investigating the effect of different growth parameters and physics. Theoretical and computational methods, however, allow separating different growth parameters and physics. An industrially applicable model needs to: (i) capture the growth conditions/morphology/ property relationship, (ii) have a high-fidelity concerning experimental observations, and (iii) be computationally efficient. Furthermore, high-throughput computational simulations of such growth models combined with emerging machine learning technologies provide the opportunity for synthesis-by-design of new layered materials. Recent developments in modeling 2D materials and their growth processes are reviewed in Ref. [6], and a perspective on their controlled synthesis is given in Ref. [7].
Isolated modeling techniques have been used to investigate the growth of 2D materials [8][9][10][11][12]. Molecular dynamics (MD) simulations based on empirical [13,14] and reactive [15][16][17] potentials are utilized to study the growth of these monolayered materials. Ab initio MD simulations are utilized to study nucleation mechanisms [18]. High-throughput DFT simulations have also been used to determine the formation energy of various 2D materials [19], and Phase-Field (PF) simulations are utilized to investigate different competing mechanisms governing their morphology [20][21][22]. Most of these studies focus on a particular aspect of growth within a limited length and time scale.
Here, we developed a multiscale/multiphysics model for the synthesis of CVD-grown monolayer materials, connecting macroscale growth parameters with the mesoscale morphology of the as-grown material. The proposed model establishes a direct connection between experimentally controllable growth parameters and the final morphology of the synthesized 2D materials. Thus, it provides an alternative to costly and timeconsuming experimental investigations.

Macro-to mesoscale model of CVD growth
Several attempts have been made to capture the various physics involved in the growth of monolayer materials across different length-and time scales. A coupled MD and finite element model was used to study thermal transport in large-area MoS 2 monolayers [23]. The multiscale Kinetic Monte Carlo (KMC) model was used to study the epitaxial growth kinetics of graphene [24]. A combination of DFT, KMC, and continuum models has also been used to understand the growth kinetics of 2D materials [25]. A coupled DFT and continuum model has been developed to study the growth of multilayer 2D materials [26]. However, a model connecting experimentally controllable parameters to final morphology and characteristics of the as-grown monolayers is still missing.
We developed a hierarchical multiscale and multiphysics model connecting the macroscale heat and mass transport in the CVD chamber with the mesoscale growth of the 2D materials. Our model at the macroscale captures the heat transfer (conduction and convection), the flow of precursor gases, flowassisted diffusion of precursors, substrate rotation, and buoyancy forces. At mesoscale, the proposed model considers the effect of orientation-dependent edge energies, deposition rate, edge diffusion, and temperature on the growth morphology of the synthesized monolayer materials.

Macroscale model of heat and mass transport
The Navier-Stokes equation models the mass transfer in the growth chamber driven by the pressure difference and buoyancy forces, Here, ρ is the mass density of the gas mixture, v the velocity field, p the pressure, I the unit tensor, µ the dynamic viscosity, g the gravitational acceleration, and superscript T the transpose matrix. The dependable variables in this equation are pressure, p , and velocity field, v.

Invited Feature Paper
In this study, we assumed ρ to be a function of temperature only, as the concentration of the precursors is low in the gas phase. Flow-assisted diffusion is governed by, where, Here, c i is the concentration of the ith species in mol/m 3 , D i is the diffusion flux for the ith component in m 2 /s, R i is the source term for the ith precursor in mol/m 2 s, v is the mass averaged velocity vector in m/s, and J i is the diffusive flux vector in mol/m 2 s.
The equations for the conservation of energy governing the convection and conduction heat transfer are where the subscript i of ϕ i denotes the islands with different orientations, ǫ is an energy scale parameter, λ is a coupling coeff icient b et we en ϕ i and u ; f is an energy functional with α being a positive constant to ensure the correct equilibrium values of ϕ i and avoid the coexistence of two different island orientations at the same location; g(φ i ) = is the local surface orientation angle, W 0 = 12 μm is a reference interface thickness, and a(θ i ) describes the anisotropy of interface thickness.
Assuming a linear relationship between the thermodynamic driving forces and the rate of change of variables, the governing kinetic equations are is the overall diffusivity, including the contribution from both the substrate diffusivity D s and edge diffusivity D e , with h({φ i }) = −1 + i (1 + φ i ) ; and F is the deposition rate. The nucleation of 2D materials is realized by introducing random circular nucleus seeds. The anisotropic kinetic coefficient L(θ i ) , similar to the anisotropic interface thickness W(θ i ) , has the form of L(θ i ) = L 0 a L (θ i ) where L 0 = 4.0E4m 2 /J s. is a reference interface kinetic coefficient while a L (θ i ) describes the anisotropy of interface mobility. Considering the different energetics and kinetic properties of the armchair, zigzag and antenna edges of WSe 2 , the interfacial thickness anisotropy a(θ i ) takes the form of The fitting coefficients a 0 , a 1 , a 2 and f 0 in Eq. (10) are fitted to align with the calculated anisotropic edge energies of WSe 2 under Se-rich conditions [27], with a 0 = 0.2653, a 1 = 2.0, a 2 = 1.073 and f 0 = 0.0001. Similarly, the interfacial mobility anisotropy a L (θ i ) has the form of Here C p is the specific heat of the gas mixture, T the temperature, and k the thermal conductivity.

Mesoscale model of the growth
We pursued a PF approach to model the morphology evolution of 2D materials. We made a set of assumptions to develop a computationally tractable model, i.e., we neglected the gas phase and surface chemical reactions and multilayer growth. We consider that the growth is realized by the continuous deposition of precursors from the gas phase to form the crystalline monolayers. Therefore, two sets of PF variables are used to describe the microstructure evolution of the monolayer during its growth. The first set is the non-conserved order parameter ϕ to distinguish the substrate (ϕ = − 1) and a crystalline island with different orientations (ϕ = 1). The second set is the conserved order parameter representing the reduced saturation field u = (c − c eq )/c eq to describe concentration distribution. Here c is precursor concentration, and c eq is its equilibrium value at given temperature and pressure. The total free energy F , as a function of these PF variables, is The fitting coefficients a L0 , a L1 , a L2 and f L0 in Eq. (11) should be fitted to align with the anisotropic edge mobility of WSe 2 under Se-rich conditions. However, due to the lack of reliable edge mobility data, we adopt the following values in the simulation: a L0 = 0.0025, a L1 = 1000.0, a L2 = 0 and f L0 = 0.000001. The phase-field simulations are performed under zero flux boundary conditions.

Application to WSe 2 growth
We adjusted the developed model to study the synthesis of WSe 2 monolayers grown by a Metal Organic Chemical Vapor Deposition (MOCVD) furnace at the 2D Crystal Consortium-Materials Innovation Platform (2DCC-MIP) [28]. Dimensions of this furnace and its different components are shown in Supplementary Fig. S1. Here, precursor gases are W(CO) 6 and H 2 Se, and the carrier gas is H 2 . We studied the effects of several growth parameters on the coverage and uniformity of the as-grown WSe 2 monolayers, including substrate rotation, diffusion coefficients, and precursors' flow rate. The boundary conditions used for the macroscale model of the heat and mass transport in the growth chamber are listed in Table 1 for the reference case study. The corresponding parameter values are listed in Table S1 in the Supplemental Materials. Experimentally measured wall temperature is used as the boundary condition on the furnace wall ( Supplementary Fig.  S2). Temperature-dependent carrier gas properties are also  shown in Fig. S3, and the PF model parameters are also listed in Table S2.

Effect of substrate rotation
The velocity profile plays a significant role in determining the precursor distribution on the substrate due to the velocityassisted diffusion mechanism. Figure 1 shows the velocity profile and streamlines within the volume above the substrate. The flow velocity is high upstream, close to the inlets, decreasing over the substrate and increasing downstream closer to the outlet. The increase in the gas velocity downstream is due to the reduction in the cross-section area (see Fig. S1 in Supplementary Materials) and domination of buoyancy forces as the gas is heated, passing over the susceptor. Comparing the streamlines for the two cases of the stationary and rotating substrate, we revealed a collected set of streamlines in the latter case. Furthermore, the substrate's rotation deflected the streamlines toward the substrate, potentially increasing the effective adhesion of precursors. The streamlines, Fig. 1, also indicate a better mixing when the substrate is rotating.   The concentration distribution of precursors, H 2 Se and W(CO) 6 , and the WSe 2 growth coverage over the substrate at steady-state conditions is shown in Fig. 2, where arrows indicate the flow direction. Our results revealed a more uniform precursor distribution over the substrate when the substrate is rotating and a higher concentration of precursors at the center that potentially leads to multilayer growth. In contrast, most of the precursor concentration in the case of the stationary substrate is concentrated at its front edge, leading to large concentration gradients, which promotes the formation of out-of-plane growth and agglomeration by Mullins-Sekerka instability [4,5]. Substrate's rotation also enhances the intermixing of precursors, providing a more uniform W: Se ratio and a more consistent final synthesis morphology over the substrate, as shown in our PF simulations in Fig. 2c-d. Reduction in the concentration of precursors downstream can be explained by the domination of buoyancy forces as the carrier gas warms up when it passes over the heated substrate. The concentration gradient of W(CO) 6 is higher than the H 2 Se, which is due to the relative position of the inlets concerning the substrate. The H 2 Se inlet is at a higher height compared to W(CO) 6 . Thus, diffusion of H 2 Se plays a more dominant role than flow-assisted diffusion, resulting in its more uniform distribution. The temperature profile of the gas close to the substrate is shown in Fig. 3, indicating no significant difference in steady-state temperature profile between the two stationary and rotating substrates.

Effect of diffusion coefficient
One of the key growth parameters is the diffusion coefficient of precursors within the gas phase. Although it cannot be directly controlled, it is a strong function of temperature and pressure. Here, we performed a systematic study to understand the role  of the diffusion coefficient. We considered four cases with diffusion coefficients ranging from 10 -9 to 10 -6 m 2 /s. The precursor concentration profiles on the substrate are shown in Figs. 4 and 5 for the stationary and rotating substrates, respectively. The corresponding PF simulations showing the growth coverage are shown in Supplementary Fig. S4.
Our results indicate that the conical deposition region of the stationary substrate widens as the diffusion coefficient increases. We also did not see a significant change in the profile of precursor concentration on the rotating substrate upon variation in the diffusion coefficient. Figure 6 shows the average concentration of precursors normalized by their corresponding mean values for the reference case D = 10 -9 m 2 /s, c * = c/ c| D=10 −9 . Furthermore, the average concentration of precursors deposited on the stationary substrate increases by increasing diffusion coefficient, Fig. 6. Comparing the normalized mean values for the stationary and rotating substrates, we concluded that the mean value of precursor species on the stationary substrate is more sensitive to the diffusion coefficient. Furthermore, the mean value of H 2 Se concentration on the substrate is less sensitive to the diffusion coefficient than the one for W(CO) 6 . Thus, the position of inlets relative to the substrate plays a significant role in reducing the sensitivity of the final growth morphology to unavoidable variations in the growth conditions.

Effect of inlet velocity
The inlet velocity of precursors and carrier gas plays a major role in determining the concentration profile of precursors on the substrate via the flow-assisted diffusion process. Our simulation results are shown in Figs. 7 and 8, and corresponding growth morphologies are depicted in Supplemental Fig.  S5. We revealed a higher sensitivity of the stationary substrate  to variations of the inlet velocity. Furthermore, reduction of the inlet velocity has a more profound effect on the distribution of the W(CO) 6 , where deposition on the front edge of the substrate drastically reduces. It is due to the positioning of the substrate at the end of a ramp and the relative location of the W(CO) 6 inlet with respect to the substrate, i.e., a lower distance with the plane of the substrate. As a result, the inlet flow of W(CO) 6 and carrier gas will deflect upon hitting the ramp, forming a boundary layer. The thickness of this boundary layer is an inverse function of the inlet and carrier gas velocity. Thus, reducing the inlet gas velocity results in a thicker boundary layer, which increases the mean free path of the precursor atoms for reaching the substrate, reducing the concentration of W(CO) 6 .
In contrast to W(CO) 6 , the H 2 Se inlet is located at a higher height. Thus, no boundary layer forms, and reduction in the H 2 Se inlet velocity directly reduces the flow-assisted diffusion of precursors and a more uniform concentration distribution of the precursors on the substrate. In the case of a rotating substrate, the formed boundary layer results in lower concentration of W(CO) 6 on the outer section of the substrate, which may hinder the growth of 2D materials in this section.
The average precursor concentration as a function of the inlet velocity normalized by the one for the reference simulations (see Table S1 in the Supplemental Materials) is shown in Fig. 9. Our results indicate a direct correlation between the average precursor concentration and the inlet velocity for a stationary substrate and higher sensitivity of the W(CO) 6 to the inlet velocity variations. Also, the variation of the precursor concentration is larger for reduction of the inlet velocity than its increase. The average precursor concentration is proportional to the inlet velocity but less sensitive to variations when the substrate is rotating.

Conclusion
We developed a coupled fluid mechanics/PF multiscale model for the growth of monolayer materials, connecting macroscale synthesis parameters to the mesoscale morphology. The model considers various physics, including heat and mass transfer in the CVD chamber, flow-assisted diffusion of precursors, buoyancy forces, surface and edge diffusion, orientation-dependent edge energies, deposition rate, temperature, and pressure. It is unique in its versatility for modeling CVD growth of various 2D materials and establishing a correlation between experimentally controllable parameters, such as precursors' and carrier gas flow rates, with mesoscale morphology of the as-grown materials.
Our results indicate a more uniform growth upon rotating the substrate, with a higher chance of multilayer growth at the center. We revealed a higher concentration of precursors on the front edge of the stationary substrate, which may lead to the formation of vertical nanofins or agglomerates due to Mullins-Sekerka instability and the Volmer-Weber mechanism. The precursor distribution on the stationary substrate is more sensitive to the buoyancy forces, diffusion coefficient, and inlet velocity. It also reduces downstream, disrupting the growth of 2D materials. We revealed a lower growth coverage on the stationary substrate and a higher growth density on a rotating substrate upon increasing the diffusion coefficient. Our results also indicate higher sensitivity of the average precursor concentration to the diffusion coefficient for the stationary substrate than the rotating substrate. We further revealed a higher deposition rate and growth coverage upon increasing the inlet velocity by increasing the velocity-assisted diffusion of precursors. The average precursor concentration is directly proportional to the inlet velocity and monotonously varies by the inlet velocity. Furthermore, the deposited precursor concentration is less sensitive to the inlet velocity for a rotating substrate. We established correlations between different experimentally controllable parameters and the coverage and synthesis quality of 2D materials. We analyzed the sensitivity of the final growth quality to experimentally controllable growth parameters and proposed several ways for reducing this sensitivity. In the future, we will couple our model with the atomistic simulations calculating the orientation-dependent edge energies as a function of gas phase chemistry. This capability will allow capturing the variation in the growth morphology as a function of precursor concentration ratios.

Methods
The macro-and mesoscale simulations are coupled via the precursor concentration distribution over the substrate. The solution domain of the macroscale model is the entire growth chamber, while it is limited to the substrate for the PF simulations. The macroscale model of the flow is solved using the finite volume method, and the PF model is solved using the finite difference technique. We have neglected the depletion of precursor on the substrate due to the island growth. Because the growth is localized on the surface and the as-grown 2D materials are only a single atom thick and thus have a subtle effect on the macroscopic transport phenomena. We further assume the deposition process from supersaturated precursors vapor.
The macroscale model is meshed using free tetrahedral elements. We refined the mesh on the substrate, the inlets, the outlet, and the adjacent volumes to improve simulations' accuracy and get mesh-independent results. We utilized the Galerkin least-squares method to discretize the governing coupled system of equations. We choose a test function that stabilizes the hyperbolic terms and the pressure term in the transport equations. Furthermore, we used the Galerkin formulations to conserve momentum, mass, and energy. Finally, we utilized the damped Newton method and time-stepping techniques with automatic time-stepping and automatic polynomial orders to resolve the velocity and pressure fields. The mesoscale PF simulations are performed using our in-house finite difference code. A minimum of five elements per interface width is considered to have mesh-independent results.