Multiscale computational understanding and growth of 2D materials: a review

The successful discovery and isolation of graphene in 2004, and the subsequent synthesis of layered semiconductors and heterostructures beyond graphene have led to the exploding ﬁ eld of two-dimensional (2D) materials that explore their growth, new atomic-scale physics, and potential device applications. This review aims to provide an overview of theoretical, computational, and machine learning methods and tools at multiple length and time scales, and discuss how they can be utilized to assist/guide the design and synthesis of 2D materials beyond graphene. We focus on three methods at different length and time scales as follows: (i) nanoscale atomistic simulations including density functional theory (DFT) calculations and molecular dynamics simulations employing empirical and reactive interatomic potentials; (ii) mesoscale methods such as phase-ﬁ eld method; and (iii) macroscale continuum approaches by coupling thermal and chemical transport equations. We discuss how machine learning can be combined with computation and experiments to understand the correlations between structures and properties of 2D materials, and to guide the discovery of new 2D materials. We will also provide an outlook for the applications of computational approaches to 2D materials synthesis and growth in general.


INTRODUCTION
The perfection and physical properties of atomically thin twodimensional (2D) materials are extremely sensitive to their synthesis and growth process. Achieving desired characteristics such as structural uniformity, high carrier mobility 1 , strong light-matter interactions, tunable bandgap, and flexibility is the main challenge for the synthesis and growth of next generation, electronics-grade 2D materials. A reliable and optimized growth and manufacturing process is essential for the synthesis of 2D materials with uniform properties at the wafer scale, e.g., for application in flexible and transparent optoelectronics.
Two main approaches have been employed for the synthesis of 2D materials, i.e., (i) top-down approaches such as mechanical 2 and liquid-phase exfoliation that allows scalability 3 , and (ii) bottom-up approaches such as chemical vapor deposition (CVD) and atomic layer deposition techniques 4 . The former approaches are suitable for mass production of 2D materials but with typically lower quality, whereas the latter approaches can produce high-quality 2D materials but in small amounts. For both types of approaches, the morphology and characteristics of the synthesized 2D materials are very sensitive to the thermodynamic or kinetic conditions 5 of the growth processes, e.g., heat transfer and mass transfer of source chemical species, chemical reaction kinetics, adsorption of reaction product species on a substrate surface, and nucleation and growth of the resulting 2D materials. The goal of this review is to provide an overview of the main theoretical and computational methods for understanding the thermodynamics and kinetics of mass transport, reaction, and growth mechanisms during synthesis of 2D materials. We will discuss the possibility of synthesis-by-design of new 2D materials guided by computation to reduce the number of expensive and time-consuming trial-and-error experimentations.
The critical challenge for developing theoretical and computational design tools for the synthesis of 2D materials is the broad range of length and temporal scales involved in their growth process. For example, it may require quantum mechanical and atomistic reactive force-field calculations to determine the activation energies for atomic migration on a surface 6 and understand the atomistic surface reaction mechanisms 7 , and then a finite element method (FEM) to model the mesoscale mass transport phenomena 8 . Other challenges include incorporating the effects of substrates including the types of substrate defects 9 , the possible wrinkling of 2D films 10 , the effect of van der Waals (vdW) interactions at the mesoscale 11 , and the growth kinetics unique to atomically thin materials 12 . Also reproducing quadratic dispersion for the flexural acoustic modes of 2D materials using classical or reactive potentials may not be straightforward, it has been already formulated 13 . Furthermore, a practically useful multiscale model should be computationally efficient, numerically accurate, and, more importantly, able to capture the multi-physical governing relationships among the growth conditions, growth morphology, and materials properties. The eventual goal of developing multiscale computational models is to guide the design of new growth chambers to produce uniform large-area 2D materials.

First-principles calculations
This section discusses the computational approaches based on the density functional theory (DFT). The primary application of DFT in modeling 2D materials is to determine the relative thermodynamic stabilities of different crystal structures of a 2D material by computing their chemical potentials or to identify kinetic pathways by analyzing the energetics of potential transient structures from one stable equilibrium structure to another as thermodynamic conditions change. Examples of applications of DFT to understanding 2D materials and their growth, as well as the corresponding DFT-based first-principles methodologies and the corresponding experimental observables for validation are summarized in Table 1.
Thermodynamic stability. The thermodynamic stability of a 2D material requires its formation energy to be negative. The formation energies of various 2D materials have been efficiently calculated by DFT in combination with high-throughput screening platforms 14 . The success of DFT in predicting energies relies on the accuracy of well-tempered approximate exchangecorrelation functionals. Standard functionals, e.g., the Perdew-Burke-Ernzerhof parametrization of the generalized gradient approximation exchange-correlation functional (GGA-PBE), accompanied by appropriate corrections 15 , describe well the formation and atomization energies. However, they are limited by intrinsic delocalization errors when bonds are stretched, resulting in underestimated reaction barriers 16 . A more realistic criterion for a stable 2D material is a low "above-hull" energy 17 . Exfoliation energies within 0.2 eV/atom are suggested as a general rule of thumb for the stability in 2D form. A slightly stricter criterion to have low surface energy (<20 meV/Å) is suggested to rule out potential high-surface energy but multiatom-thick sheets.
Kinetically stabilized 2D Materials. As-grown products may be metastable structures rather than ground-state polytypes 18 . The well-known phase-stability competition between the 2H and 1T′ types in transition metal dichalcogenides (TMDs) has been studied extensively using DFT. Several possible mechanisms for the 1T′phase stabilization in the disulfide, including growth-related factors such as the presence of point defects and residual strain, have been elucidated through DFT calculations 19 . Phase-stability analysis for the entire group -IV metal dichalcogenide family 20 revealed that the 1T′ phases are generally more stable for ditellurides. In addition to pure compounds, mixing/segregation behavior and order-disorder transitions in alloys have also been discussed 21 . Within the W/Mo + S/Se/Te combinations, mixing either chalcogens or metals is thermodynamically favorable 22 , which can be experimentally quantified using the Warren-Cowley order parameters 23 . For TMDs, theory predicts that at finite temperatures, the formation of random alloys is favorable. However, kinetically stabilized atomically thin strips of alternating W and Mo in a sulfide alloy have been reported 24 .
Interlayer/substrate interactions. Accurate description of vdW interactions becomes critical when interlayer sheet-sheet or sheet-substrate interactions are considered, e.g., for identifying the growth orientation. As a non-local effect, vdW forces are not correctly described by DFT with semi-local exchange-correlation functionals and are often corrected by adding pairwise interatomic terms from empirical fittings 25 , based on charge densities 26 , or by introducing non-local exchange-correlation functionals, i.e., van der Waals Density Functional 27 . Thorough testing of these methods has been examined in ref. 28 . One common technical problem in this computational approach is that substrates are modeled in a slab geometry where the one surface not in contact with the 2D sheet may host surface states, unless it is passivated. To suppress artificial charge transfer that these may induce, capping using pseudohydrogen 29 with fractional charge is often performed. With the advent of new software platforms automating the generation of solid surfaces in combination with stacked 2D sheets and adsorption geometries (e.g., MPInterfaces 30 ), highthroughput screening of possible substrates may lead to a systematic approach to substrate engineering.
Precursor chemistry and kinetics. Transient intermediate states are commonly calculated using transition state theory and the nudged elastic band (NEB) method 31 by calculating bond dissociation energies or corresponding activation barriers 32 . One recently attempted approach to capture precursor reaction kinetics is constrained molecular dynamics (MD) 33 , where slowly varying coordination constraints are enforced on reacting species to map out free energy barriers, as implemented in VASP 34 . It is different from NEB calculations in that (1) it calculates free energies at finite temperatures; (2) the final state of the reaction need not to be known; and (3) it has a better numerical stability 35 . Constrained MD was used to study a sulfur precursor, S 2 , reacting with a MoO 3 surface 36 , concluding that MoO 3 surface vacancies favor the sulfurization process both kinetically and thermodynamically. The downside of constrained MD is that a reaction coordinate is chosen a priori, which may bias the system towards unnatural products with incorrect reaction mechanisms 33 and higher reaction barriers 34 .
Other methods for sampling of rare events such as nucleation of a new structure include umbrella sampling, transition interface sampling, and metadynamics. The umbrella sampling technique was introduced by Torrie and Valleau 37 , to improve the sampling of systems with energy landscapes containing high energy barriers. The weighted histogram analysis method 38 can be used to analyze a series of umbrella sampling methods. In the transition interface sampling 39 , the transition region is divided into subregions of intermediate states. The rate constant of a reaction in this method is the multiplication of transition probabilities between different intermediate states. The metadynamics technique was introduced in 2002 and is usually used within an atomistic modeling framework 40 . Selected collective variables of the system not only evolve with time but also periodically leave behind positive Gaussian potentials that are added to the original Table 1. Examples of applications of DFT to 2D materials, the DFT-based first-principles methodology, and corresponding experimental observables.

Aspect of growth
First-principles framework Experimental observable and impact Thermodynamic stability Thermochemistry, high-throughput screening Successful synthesis potential energy of the system, which effectively push the system out of a local minimum and into a neighboring energy well. Informally, metadynamics resembles "filling the free energy wells with computational sand." 41 Growth front advancement. Early studies focused on calculating the surplus energies of edges in the high-symmetry directions and then applying the thermodynamic Wulff construction 42 . One difficulty in calculating the edge energies is the polar and noncentrosymmetric nature of materials, where the usual ribbon calculation geometry makes the energies of dislike edges inseparable. The usual workaround is by constructing triangular flakes exposing a single type of edge 43 . Another "energy density method" was introduced in ref. 44 and a recently proposed method 42 aimed at finding a general method suitable for high-throughput calculations introduces capping groups to passivate a surface. One major criticism of further employing thermodynamic Wulff constructions is that growth is by definition out of equilibrium. Therefore, edges that dominate over others should not be the energetically favorable ones, but the slowest growing ones. The step-flow approach was formulated in ref. 45 for graphene growth and was further developed in ref. 46 for polar materials, establishing the use of kinetic Wulff construction. A similar approach involves DFT calculations and fitting to experimental grain morphologies, to construct a kinetic Monte Carlo (KMC) model 47 .
Defects. The formation of defects may occur within thermodynamic equilibrium such as thermally generated point defects or follow from growth imperfections such as dislocations and grain boundaries (GBs), or reflect the finite size of crystals such as edges and surfaces. Defects can also be deliberately introduced using methods such as electron or ion irradiation and chemical treatment 40 . The calculation of defect formation energies follows a well-established procedure detailed in ref. 48 . 2D materials also host lattice-specific defect types, such as the Stone-Wales defect in graphene 49 . Single-atom vacancies are another type of defects in 2D materials. A general strategy to ensure correct convergence behavior for charged defects in 2D systems is presented in ref. 50 . Defect formation energies in MoS 2 have been comprehensively studied in ref. 51 ; possible strategies of introducing extrinsic dopants have been examined in ref. 52 . Defect complexes are also frequently studied to identify likely combinations between simple intrinsic defects and external contamination 53 , devise possible defect-pairing strategies to neutralize harmful defects 54 , and investigate their influence on the growth behavior of 2D materials on a different 2D sheet 9,55 . Multiple vacancies may exist in 2D materials such as double vacancies in graphene, resulting in (i) two pentagons and one octagon-V 2 (5-8-5) defect; or (ii) three pentagons and three heptagons-V 2 (555-777) defect; or (iii) four pentagons, a hexagon, and four heptagons-V 2 (5555-6-7777) defect 56 . The formation of defects with an even number of vacancies in graphene is energetically favorable due to the lack of any dangling bonds 57 , whereas a large number of vacancies may bend and wrap the 2D material 58 . Table 2 lists the typical values of formation and migration energies of these point defects from atomistic calculations.
Dislocations are line defects in 2D materials that can be of inplane edge or out-of-plane screw type. The former consists of pentagon-heptagon (5-7) pairs 59,60 and the latter makes the structure to become 3D 61 . In general, dislocations in 2D materials have high gliding energy barriers and are immobile at room temperature. They hinder the plastic deformation of 2D materials 62 . The in-plane edge dislocations with the smallest Burgers vector, b = (1,0), is typically the energetically most favorable dislocation configuration. The trilayer bi-atomic MoS 2 has a complex dislocation structure with two types, i.e., Mo-rich that resembles (5-7) and S-rich that resembles  pairs. The complex nature of dislocations in MoS 2 is a source of its rich chemistry. For example, at specific chemical potentials of S, the (5-7) dislocation reacts with 2S vacancies to form a (4-6) defect.
GBs in 2D materials can be considered as an array of dislocations, which can be divided into low-angle and highangle GBs depending on the dislocation core density. Atomistic simulations revealed that the mechanical properties of 2D materials, similar to the one-dimensional nanostructures, will be affected by both the density and arrangement of the defects 63 . Notably, graphene sheets with high-angle tilt GBs can be as strong as the pristine material and much stronger than those with lowangle GBs 64 . First-principles calculations also revealed that sinuous GB structures are energetically favorable when the straight GB line cannot bisect the tilt angle 65 . GB energies of different 2D materials are also reported based on atomistic calculations. For MoS 2 , the energy of GB with a 20.6°misorientation angle is 0.05 eV/Å 62 . The GB energy for phosphorene varies between 0.05 and 1.5 eV/Å, depending on its tilt angle 66 .

Molecular dynamics
The MD approach is one of the main tools that have been used to study the growth of 2D materials at the atomistic scale. The quality of potential function governing the atomic interactions is the key for the accuracy of this technique in predicting the structures and properties of 2D materials. In general, two classes of interatomic potentials are used in MD simulations to study the growth of 2D materials, i.e., empirical and reactive potentials, which we elaborate on here.
Empirical potentials. Empirical potentials are faster and more computationally efficient, but less accurate than reactive potentials. Common empirical potentials for modeling 2D materials are reviewed below.
Lennard-Jones potential: In 1924, Lennard-Jones (LJ) published the 12-6 pairwise potential to describe the vdW interaction between two neutral particles (i.e., rare gases) 67 , where r is the distance between two interacting atoms, ϵ is the potential well depth and σ is the distance at which the potential is zero. The r −12 term captures the strong repulsion occurring when two atoms get closer. The repulsive effect originates from the Pauli exclusion principle, which penalizes the overlap of electron orbitals. The r −6 term describes the vdW interaction, which is due to the coupling between the instantaneous polar charges induced in two neutral molecules and is significantly weaker than the electron overlap effects. The popularity of the LJ potential is due to its simplicity and it has been widely applied to simulate the interlayer interactions in atomically layered materials. While being used to capture vdW interactions, a limitation of the pairwise interaction is the inability to capture frictional processes, such as the sliding of atomic layers.  [68][69][70] . The two potential parameters ϵ and σ can be determined by using one energy quantity like the binding energy and a structural quantity like the density. The LJ parameters for most elements can be determined systematically, as demonstrated by Rappe et al. 71 . The LJ parameters for typical 2D materials and substrates in Table 3 are listed in in Table 4 (for σ) and Table 5 (for ϵ), respectively. Stillinger-Weber potential: Stillinger and Weber (SW) proposed this potential in 1985, to model bulk silicon 72 . The SW potential includes both pair (two-body) and three-body interactions, (2) where V 2 describes the pairwise bond-stretching energy and V 3 captures the bending energy associated with an angle with initial value θ 0 .
The cut-offs r max , r max12 , and r max13 can be determined by fitting to the average of the firstand second-neighboring bond lengths. The seven unknown potential parameters, including five unknown geometrical parameters, i.e., ρ and B in V 2 and ρ 1 , ρ 2 , and θ 0 in V 3 , and two energy parameters A and K, can be determined either by a least-squares fit to targeted quantities such as the elastic constants, Poisson's ratio, and the phonon spectrum, or by analytical derivations from the linear valence force-field potential. It is noteworthy that the nonlinear parameter B could not be determined by the latter approach, which should be determined by nonlinear quantities such as the third-order elastic constants 73 .
The SW potential can describe nonlinear processes including large deformations, as higher-order terms for the variation of bond length and angle are included in its formulation. Although the SW potential is efficient and suitable for simulating thermal transport and other nonlinear phenomena 74 , it is not able to describe chemical bond formation or breaking. It should also be noted that the SW potential cannot provide bending energy for atomically layered materials without out-of-plane bonds such as graphene 75 . Due to its efficiency, the SW potential has been used in the simulation of many atomically layered materials such as 2D nanostructures of Si and their thermal properties 76 2 77 , and black phosphorus 73,78,79 . The SW potential has been parametrized for 156 emerging atomic layered materials, along with other available empirical potentials for these atomically layered materials 80 .
Force-field potential: The force-field (FF) potential describes the variation in the potential energy of a deformed structure with respect to its equilibrium configuration. It is essentially a Taylor expansion of the total energy in terms of the variation in the bond length and angle, where Δb, Δθ, Δϕ, and Δγ are variations of the bond length, bond angle, twisting angle, and inversion angle, respectively; θ , k ϕ , k γ , k bb′ , k θθ′ , and k bθ are force constants. V b is the bond-stretching energy, V θ is the bond angle-bending energy, V ϕ is the torsional energy for a bond, V γ is the inversion energy among four atoms, and V c is the off-diagonal coupling interaction. It is noteworthy that the cubic and quartic anharmonic terms in the bond-stretching and angle-bending interactions enable the FF model to capture nonlinear phenomena such as thermal expansion and thermal transport.
The FF potential has two significant features: (i) the number of interaction terms and (ii) the force-constant parameters and the equilibrium structural parameters (bond lengths and angles). MM3 81 and COMPASS 82 use all of the interaction terms and can simulate large deformations and nonlinear phenomena. In contrast, the Keating model 83 is the simplest model including only the harmonic terms. Parameters of the potential are usually fitted to available data and, thus, the volume and accuracy of the fitted data set determine the accuracy of the FF model. Several generations of FF models have been developed as more data become available, including four versions of the molecular mechanics FF: MM1 84 , MM2 85 , MM3 81 , and MM4 86 . The UFF 71 provides the FF parameters for all the elements in the periodic table.
As the energies in the FF model are additive, this potential has been widely used to calculate physical or mechanical properties for various materials via analytic expressions, including many 2D materials. For example, an FF model, with bond stretching and angle-bending terms, has been used to derive the phonon dispersion of graphene layers 86 . The bond-stretching and anglebending terms were utilized to extract the bending properties of graphene 87 . A mapping between the FF model and the elastic beam model has been developed, which is used to derive analytic expressions for Young's modulus and Poisson's ratio of graphene and carbon nanotubes 88 .   92 , and Reactive Force Field (ReaxFF) 93 . These different reactive potentials are reviewed in detail in ref. 94 . Performances of several reactive and non-reactive potentials in graphene-based materials are compared in ref. 95 . These potentials allow the realization of the thermodynamic stability of graphene-based materials under different conditions, environmental chemistry, substrate, interlayer interactions, and structural defects. However, only a few of them have been employed in the study of computational synthesis of other 2D materials to reveal thermodynamic stability, simulate the process of the growth of layers, model the top-down synthesis of 2D materials such as exfoliation techniques 38 , and evaluate the properties of 2D nanostructures. This is due to several challenges including the following: (1) the long time required to develop a reactive potential, (2) the nontransferability of potential parameters from one material to another, and (3) the lack of universal functional forms for the reactive potentials. The Tersoff potentials are mainly used to explore various structures and properties of 2D hexagonal boron nitride (h-BN) materials 96 . The REBO and AIREBO potentials are used to simulate the CVD process and growth of 2D amorphous carbon on substrates, mechanical properties of MoS 2 , and even study the thermal stability of C 60 2D nanostructures 97 . The COMB3 potentials have been used to simulate the experimental CVD deposition and growth of graphene and metal on substrates 98 . The ReaxFF potentials have been used to explore the structures of several 2D materials under intercalation, study the defects and groups on the surfaces of MXene materials, simulate CVD growth of MoS 2 layers, and explore the synthesis of 2D polymeric materials in experiments 7,99 . In general, COMB3 and ReaxFF can further explore 2D materials, as their functional forms can describe heterogeneous systems. It remains crucial to find accelerated methods to develop reactive potentials with high accuracy and transferability for 2D materials.
Tersoff potential: The Tersoff potential was first published in 1986 100 and further developed in 1988 101 . The second version of the Tersoff potential has the following form, where f R ðr ij Þ and f A ðr ij Þ are the repulsive and attractive components, respectively. The cutoff function is: A distinct feature of the Tersoff potential is the absence of an explicit three-body interaction term. Instead, the concept of bond order is introduced, in which the strength of a given bond depends strongly on its local environment. Bond order is accounted for via the function b ij , which has the following form a ij is an alternative parameter to improve the accuracy of the potential, with a similar form The 13 unknown parameters are typically determined by fitting to data from first-principles calculations or experiments.
The Tersoff potential was used for modeling of graphene, which gives a better description of the phonon spectrum 102 . It was also used to model germanene 103 and silicene 103 . By introducing one scaling parameter for the quantity b ij with atoms i and j of different elements 104 , the Tersoff potential was generalized to describe the interaction for multicomponent systems including C-Si and Si-Ge 104 , and boron nitride 105 . Using the Tersoff potential, the thermal conductivity of silicene was calculated to be 61.7 and 68.5 W mK −1 for the armchair and zigzag configurations, which is reduced when it is doped with heavier silicon isotopes 106 . This potential has also been used for modeling the misfit dislocations in 2D materials, where a critical thickness was found for the formation of an interface misfit dislocation 107 .
ReaxFF potential: The ReaxFF 93 is a bond-order-based force field, which allows formation and breaking of bonds. It has a general form of where the energy terms on the right-hand side in the order of appearance are bond, over-coordination penalty, undercoordination stability, lone pair, valence, torsion, non-bonded vdW interactions, and Coulomb energies, respectively. The bondorder parameter determines the bonded interactions between the where BO σ ij ; BO π ij ; and BO ππ ij are the partial contributions of associated σ, π, and ππ bonds between the i and j atoms; r ij is the distance between the i and j atoms; r σ o , r π o , and r ππ o are the associated bond radii; and p bo are empirical constants that will be determined from ab-initio simulations or experiments. The energy terms in Eq. (14) are then formulated as a function of the bondorder parameter (see ref. 93

for details).
ReaxFF studies on 2D materials were initiated from studies of various aspects of graphene 108 and inspired investigations for a number of other 2D materials, which are elaborated below.
Graphene-Effect of the environment's chemistry on properties and responses of graphene has been investigated using ReaxFF. For example, effect of oxygen on the mechanical response of graphene 109 , effect of the addition of fluorine and hydrogen 110 on the structural quality of graphene, as well as the Ni-catalyzed growth process of single-layer graphene 111 , are investigated. The growth mechanism involves atomic deposition and segregation of C atoms on the Ni substrate, forming dome-like carbon nanotube cap and polygonal ring patterns via atoms→chains→graphenelike layer formation (Fig. 1a). The ReaxFF has also been applied to study reactive events at the graphene-water interface (Fig. 1b) 112 and selective desalination via the bare and hydrogenated graphene nanopore 113 .
Silicene and phosphorene-Silicene and black phosphorene (BP) are the silicon and phosphorus analogs of graphene, respectively. A favorable route for silicene formation has been proposed by employing the graphene bilayer as a template (Fig.  1c) 114 . It was found that vacancy defects reduce the thermal stability of silicene, where the critical temperature reduces by more than 30% 115 . A P/H parameter set has been developed for monolayer BP 114 to examine the mechanical, thermal, and chemical stability of both pristine and defective BP.
Molybdenum disulfide-The ReaxFF potential for 2D MoS 2 systems was introduced in ref. 6. This potential could precisely calculate the formation energy of ripplocations in MoS 2 (Fig.  1d) 116 , which result from the slippage of the upper layer against the lower layer to relax the strain without breaking covalent bonds. This potential was also trained over the vacancy formation energies, diffusion barriers, bending rigidity, and kinetics. It has also been applied to design S vacancy defects on MoS 2 surface 117 . Reactive MD simulations have also been used to model the synthesis of MoS 2 layers via sulfidation of MoO 3 surfaces 7 . These simulations identified a three-step reaction pathway as follows: (i) evolution of O 2 and reduction of MoO 3 surface; (ii) S 2 -assisted reduction and formation of SO/SO 2 ; and (iii) sulfidation of the surface via Mo-S bond formation. The ReaxFF MD simulations have also been employed within a multiscale framework to capture thermal properties of MoS 2 118 . The predicted thermal conductivity of MoS 2 structures with infinite length (37 ± 3 W/mK) falls within the range of experimental results (34.5-110 W/mK) 119 .
MXenes-MD simulations with ReaxFF of MXenes revealed that intercalating with potassium cations drastically improves water stability and homogeneity of MXenes, and reduce the selfdiffusion coefficient of water by two orders of magnitude 120 . The dynamical response of MXenes with different surface terminations to intercalating ions was studied with DFT calculations and ReaxFF MD simulations 121 . Figure 1e compares the interlayer distance of Ti 3 C 2 (OH) 2 MXene bilayer with and without K + ions. Influence of metal ions intercalation on vibrational properties of water molecules trapped between MXene layers was studied using ReaxFF MD simulations 122 . Widening of the interlayer gap may enable the penetration of molecular reactants such as urea, which decomposes readily and leads to the intercalation of ammonium cations 123 . ReaxFF simulations were performed for defect formation and homoepitaxial growth of TiC layer on Ti3C2 MXene structures 124 .
REBO/AIREBO potentials: The bond-order interaction ratio, b ij , is dependent on the local coordination in the atomic environment. The REBO potential can define the conjugation of bonds using b ij , and therefore the potential has been widely used for modeling hydrocarbon systems 125 . The repulsive term is expressed by the Brenner equation 90 , Here, the Q ij , r ij , and α ij parameters depend on i and j. The bond weighting parameter w ij (r ij ) depends on switching function S(t) as The attractive interaction is expressed as the bond-order interaction ratio b ij is Here, p σπ ij and p σπ ji are not necessarily equal as they depend on the penalty function (g i ) of bond angle θ jik between the vector r ij , and vector r ki , i.e., The AIREBO potential is the modification over the REBO potential, which also considers the torsional ðE TORSION kijl Þ and LJ ðE lj ij Þ interactions. The integration of pairwise interactions in AIREBO potential can be represented by the following equation The REBO potential can describe the breaking, formation, and hybridization of covalent bonds. The REBO and AIREBO potentials have been used to study the growth of 2D amorphous carbon, mechanical properties of MoS 2 , and even thermal stabilities of C 60 2D nanostructures using atomic coordination, internal stress, and mass density analysis 126 . The subplantation phenomenon during CVD, i.e., penetration of the substrate bulk carbon by carbon atoms, can be explained by this REBO potential. REBO potential is also suitable for modeling and analysis of TMD materials such as MoS 2 . The potential results in stabilized structure and good agreement of the structural and mechanical properties 97,127,128 . It was also used for modeling the synthesis of diamane from graphene and its stability under different conditions 129 .
COMB3 potential: The general formation of the COMB3 potential can be expressed as 92 The total potential energy consists of electrostatic energy , and the correction terms. Here, q and r represent the charge and coordinate array of the elements. The electrostatic energies U es q f g; r f g ½ ð Þ can be represented as The self-interaction energy U self q f g; r f g ½ À Á consists of fieldeffect and ionization energy, or the affinity energy. Charge-tocharge interaction U qq q1s f g; r f g ½ ð Þ is associated with the charge density distribution functions. Nuclear-charge interaction U qZ q f g; r f g ½ depends on the coulombic interaction between charge density distributions in the system. Polar interaction energy U polar q f g; frg ½ depends on the charge interactions and dipole distributions in response to the external electric field.
The charge-dependent short-range interaction in Eq. (24) is associated with bond distance functions of the atom and the associated charge functions. This energy term can be expressed by Here, V bond ij is the bond energy that is associated with pairwise attraction V A ðr ij ; q i ; q j Þ À Á and pairwise repulsion energy ðV R ðr ij ; q i ; q j ÞÞ. They exponentially decrease with the interatomic distance r ij . F c (r ij ) is the cutoff function.
The long-range interaction energy is defined by the vdW interaction, which can be expressed by the LJ interaction formula 130 K. Momeni et al.
Here, σ vdW ij and ε vdW ij are the equilibrium distance and interaction energy, respectively; r ij is the cutoff radii where the interaction assumes zero.
The correction factor in Eq. (24) adjusts the energy penalties associated with specified angles. It can be calculated by a sixthorder Legendre polynomial and bond bending term.
Here, K θ ijk indicates the specific bond angles and θ ijk is the bond angle.
The COMB3 potential has been introduced to study the CVD process for various 2D materials. The charge-optimized COMB3 potentials have been used to simulate the experimental CVD deposition and growth of graphene on metal substrates. Especially, the introduction of dihedral interaction terms in the COMB3 potential can capture the delocalized bonding and bond bending during the CVD process on a particular substrate. The wrinkle formation due to the size mismatch of graphene and Cu-substrate has been studied using COMB3 potential, which was comparable to experimental results 98 . This potential was also used to understand the role of growth parameters such as absorption energy, migration barrier, and temperature in the Cu deposition on ZnO substrate 131 . COMB3 potential has also been used to investigate the effect of surface hydroxylation of amorphous SiO 2 substrate on heat conduction of supported graphene 132 .

MESOSCALE SIMULATIONS
Phase-field model At mesoscale, a simple and effective model for 2D materials growth is the classical Burton-Cabrera-Frank (BCF) theory 133 , where the island growth is realized by the deposition of adatoms from the supersaturated gas atmosphere and the incorporation of diffusing atoms at the step edge. The direct numerical implementation of the BCF model requires tracking the evolving step edges and applying the boundary conditions due to the sharpinterface nature of the model; in contrast, phase-field diffuseinterface models 134 avoid the explicit tracking of the interfaces.
The BCF-based phase-field models have been developed since the 1990s for a series of problems related to crystal growth, including collective step motions in a 1-D step train 135 , spiral surface growth during thin-film epitaxy 136 , and step flow under different kinetic regimes 137 . More recently, the combined effects of edge diffusion, the Ehrlich-Schwoebel barrier, deposition, and desorption for epitaxial growth have been investigated in 2-D 138 . Based on these studies, Meca et al. 139,140 applied the model to the growth of 2D materials with anisotropic diffusion on substrates, using the epitaxial graphene growth on copper foil during CVD as an example, with experimental validations. The typical BCF-based phase-field model for 2D materials growth uses a smoothly varying order parameter ϕ to distinguish the phases, e.g., ϕ = 1 in the 2D island, ϕ = −1 in the substrate, and −1 < ϕ < 1 at step edges. The free energy functional F of the system has the general form of: where the first term f(ϕ) is a double-well (or multi-well, e.g., see refs. 136,141 ) function with global minima at ϕ = 1 and ϕ = −1; the second term describes the coupling with the reduced saturation field u, with g(ϕ) an interpolation function and λ the coupling coefficient; the third term is the gradient energy due to the inhomogeneous distribution of ϕ at interfaces, where W(θ) is the angle-dependent interface thickness with θ ¼ arctanðϕ y =ϕ x Þ. u can be directly related to the adatom concentration c through u ¼ Ωðc À c 0 eq Þ, where c 0 eq is the equilibrium concentration for a flat interface and Ω is the atomic area of the solid phase. The evolution of ϕ is governed by where τ ϕ is the characteristic time. The evolution of the saturation field u is governed by where D is the diffusion coefficient, F d is the effective deposition rate, and τ s is the characteristic time scale for atom desorption. The phase-field equations (Eqs. (30) and (31)) can recover the sharp-interface BCF model in the limit of W(θ) → 0 140 . The applications of the phase-field models in 2D materials growth can be categorized into the following three aspects. (1) Reproducing and explaining experimentally observed growth morphologies. Examples include the reconstruction of the spiral growth of SnSe 2 141 , the verification of impurity-induced bilayered graphene growth 142 , and the explanation to island shape change from quasi-hexagons to triangles and dendrite-like morphologies 143 . (2) Investigating the effect of experimental parameters on growth morphologies. Taking the CVD growth of graphene as an example, the effects of flux of carbon species 144 , deposition rates 145 , and equilibrium saturation 145 , and substrate orientation 139 on graphene island morphologies have been investigated using phase-field simulations, which could provide useful guidance to the experimental control of the growth qualities. (3) Further development of the model to include the effect of additional physical/chemical processes during growth. In addition to the deposition, edge diffusion, and edge anisotropies, phasefield simulations have been extended to model multi-island interactions 146 , tilt GB topology of 2D materials on substrates with topological cones 147 , and to include the chemical reaction kinetics in the gas phase using a microkinetic model 148 . These works improve the fidelity and applicable range of the phase-field model for 2D materials growth.
Nevertheless, a primary obstacle to applying phase-field models for computational synthesis of 2D materials is the lack of connection between the model parameters and the experimentally controllable CVD processing parameters, which is largely due to the difference in size scales of interest. Most of the existing phase-field simulations only focus on the growth of a few 2D islands (up to micrometer scale), which can hardly represent the effect of the macroscopic CVD parameters. One possible solution is to enable the "multiscale" modeling approach 8 , integrating the FEM (see the section on "Macroscale models") to account for the effect of macroscopic CVD parameters on transport phenomena and the mesoscale phase-field model for 2D crystal growth within selected zones of the whole substrate. In particular, the steadystate velocity, concentration, and temperature profiles from FEM can be incorporated in the phase-field model, as it has been proposed and applied by Momeni et al. 8 (see Fig. 2). The incorporation of chemical reaction kinetics in the phase-field model can improve the accuracy of the predictions and expand its applicability. Phase-field models have also been applied to phase changes in 2D materials. Notable examples include the multidomain microstructure during H-T' structure transformation in MoTe 2 and the domain switching behaviors under external stimuli 149 . A fast Fourier transform algorithm has been implemented to study the domain switching during bending of MoTe 2 (see Fig. 3 149 ).
Phase-field-crystal model Different from the continuous mesoscale phase-field approach, the phase-field-crystal (PFC) approach 150 describes the thermodynamics and dynamics of phase transformations through an atomically varying order parameter related to the atomic density field. With the application of the classical dynamic DFT, PFC can capture the atomistic-scale morphology and evolution dynamics, within diffusional time scales, which is usually challenging for conventional atomistic methods. Recently, PFC has been actively applied to investigate the defect formation and the atomic structure of interfaces, GBs, and triple junctions during the growth of 2D materials 151 . The unique GB structure and collective domain dynamics in binary 2D materials, e.g., h-BN, was investigated (see Fig. 4 152 ).
Compared with the continuous phase-field approach, the advantage of PFC lies in resolving atomistic-scale interface structures and growth dynamics of 2D islands; however, due to the limitation of simulation size, interface dynamics under the effect of realistic CVD parameters are difficult to realize. Future directions of PFC simulations for 2D materials growth include the consideration of atomistic-scale island-substrate interactions and the integration with the larger-scale phase-field or FEM simulations.
Kinetic Monte Carlo KMC is another mesoscale technique for predicting the growth morphology and kinetic mechanisms of 2D materials, where all possible kinetic events are listed in an event catalog, and their stochastic sequence are randomly selected based on their activation energies 153 . Activation energies of the related kinetic processes are usually obtained from atomistic-scale calculations. KMC identifies the governing kinetic pathways and simulates the kinetic process. Details of its implementation can be found in ref. 154 .
KMC has been extensively applied to simulate the growth of graphene, e.g., the evolution of vacancy complexes and formation of vacancy chains 155 , the formation and kinetic effect of multimember carbon ring complexes 156 , etching of graphene GBs due to oxygen migration and reaction 157 , and GB evolution following the Stone-Wales mechanism 158 . Growth of nanoscale graphene islands   from carbon monomer nuclei or pre-existing growth fronts has been captured, showing catalytic growth behaviors 159 , anisotropic morphological patterns 160 , temperature-and deposition-flux-ratedependent size and shapes 161 , inhomogeneous and nonlinear growth kinetics due to lattice mismatch 162 , and geometrydetermined growth mechanisms 163 . The "coastline" graphene morphology during sublimation 164 and the step-flow growth of epitaxial graphene have also been reported 165 . Based on KMC simulations for graphene growth on Cu (111) with and without hydrogen, the growth protocol was designed for bilayer growth and N-doped graphene growth 166 .
Regarding TMDs, the KMC simulations for the growth of WSe 2 monolayer on graphene have been used to develop the phase diagram of domain morphologies as a function of flux and precursor stoichiometry (Fig. 5) 154,167 . It was found that the fast kink nucleation and propagation, rather than edge attachment and diffusion, could lead to ultrafast growth of monolayer WSe 2 168 . KMC simulations also guided CVD growth of large-scale WSe 2 grains by controlling the three-stage adsorption-diffusion-attachment mechanisms 169 . KMC model is also a key component in a more generalized mechanistic model for growth morphology predictions of 2D materials 47 .
In summary, the KMC model can be a useful tool for investigating the kinetic pathways and morphologies during the growth of 2D materials. However, the probability-based nature of KMC makes it most suitable for cases where atomic fluctuations are high, i.e., the atomistic and nanoscale morphology and kinetics of 2D materials. For larger-scale simulations, due to the significantly increased system size and the disparate rates of different KMC events, a full-KMC is computationally expensive. As a compromise, multiscale KMC has been developed 163 ; simplifications should also be made to account for the key events that are most relevant to the large-scale kinetics 163 . Under such situations, the phase-field approach could be a more efficient option. Meanwhile, the KMC model can provide the governing kinetic mechanisms for phase-field simulations to improve the validity and accuracy of the simulation.

MACROSCALE MODELS
Although the growth of 2D materials occurs at the nano-and mesoscale, it is controlled by physics and parameters that have a macroscopic nature, e.g., heat and mass transfer, furnace configuration, and gas-phase reactions. Thus, having a thorough understanding of the macroscale physics and processes is essential for controlling the growth of 2D materials and their synthesis by design. We can classify the macroscale models of the growth chamber into four groups: (i) experiment-based models, where rate equations and their constants used to describe the growth 170 are determined from experiments; (ii) analytical models, where the governing equations are simplified and solved analytically 171 ; (iii) adaptive models, where a set of experiments are used to train the model 172 ; and (iv) multiphysics models, where the coupled system of governing equations at different length and temporal scales are solved numerically 173 . Among these methods, the last group of models have key advantages, providing a profound understanding to the growth process, flexibility to apply to different growth conditions, and the ability to optimize the process.
A practical macroscale model of the growth chamber should capture the critical governing physics, e.g., the heat and mass transports and chemical reactions. Setting up these models requires several key information that can be obtained from lower scale simulations 7 or experiments 174 . Identifying the gas-phase reactions, we may approximate some of the reaction parameters using classical theories, e.g., the collision rate to estimate the chemisorption rate of species 175 or the group contribution methods to determine the diffusion coefficient 176 . The other approximation is for low concentrations of reactive species, where the change in pressure and heat of reaction can be neglected. In the latter case, we may decouple the fluid and heat transfer of gaseous materials from the mass transfer and kinetics. In contrast, for high concentrations of reactive species in the gas phase, such as in metalorganic chemical vapor deposition (MOCVD), the coupled system of equations must be solved 176 . A summary of main equations in macroscale models are presented below.

Gas flow
The Navier-Stokes equation governs the flow rates in the growth chamber, where u is the velocity field, p is pressure, I is the unit matrix, μ is the dynamic viscosity, and F is the volumetric applied force, i.e., Adapted from ref. 167 .
weight. Furthermore, the density ρ is a function of temperature and precursor concentration. In general, the concentration of precursor in the gas phase is low and ρ can be assumed to be only a function of temperature.

Heat transfer
The two main heat transfer mechanisms in the growth chamber are convection and conduction that respectively are where Q is the heat source, C P is the heat capacity at constant pressure, T is the temperature, q is the heat flux vector, and k is the thermal conductivity coefficient.

Mass transfer
Commonly, the concentration of precursor materials in the growth chamber is negligible. Thus, the equation governing the mass transfer can be formulated with the flow-assisted diffusion of dilute species where R is the source term for precursor, c is the precursor concentration, D is the diffusion coefficient, and N is the flux of the precursor. Equations (32)-(34) must be solved for a set of initial and boundary conditions, as shown in Table 6.
Chemical reaction Consider a set of j reactions involving i species of the form aA þ bB þ Ð xX þ yY þ with forward and reverse reaction rate constants k f j and k r j , respectively. These reaction rates depend on the temperature via an Arrhenius expression, e.g., Here, A f j denotes the frequency factor, n the temperature exponent, E f j the activation energy, and R g the gas constant. According to the mass action law, the reaction rates r j can be described as, where c i is the concentration of i-th species and v ij is the stoichiometric coefficient (negative for reactants and positive for products). The rate equations provide the information regarding the concentration of different species, which in combination with mass and heat balance equations form the complete system of equations.
Macroscale models provide insight into the growth mechanisms and flow regimes, and play a crucial role in determining proper furnace design and growth parameters. Figure 6 shows the application of a macroscopic model of the flow and concentration diffusion in a vertical furnace with a side inlet to study the effect of the gap between the susceptor and substrate. It shows that the presence of the gap results in the formation of eddy currents close to the substrate edge, which disturbs precursor concentration and hinders the formation of 2D materials that is consistent with the experimental observations. These models have been successfully used in determining the mechanisms governing the growth of 2D materials. For example, it was revealed that the concentration gradient of the precursor is a critical factor in determining the growth morphology of 2D materials, where a low concentration gradient along the substrate leads to the formation of multigrain 2D films. In contrast, moderate planar concentration gradients lead to isolated islands of varying morphologies, and an out-ofplane gradient leads to the formation of standing 2D materials 173,177 . Furthermore, macroscale models can be combined with lower scale models to form a multiscale platform with high fidelity, e.g., ref. 8 .

MATERIALS GENOME APPROACHES
Data-driven methods have been employed to satisfy the urgent need for theoretical development of 2D materials. The improvement of computational power has resulted in extensive collections of data and scientifically organized databases, which can be used by methods such as data mining and machine learning, to accelerate the computational design of new 2D materials.
Database of 2D materials Compared with the databases of 3D bulk materials such as the Crystallographic Open Database 178 , the Inorganic Crystal Structure Database 179 , and the Cambridge Crystallographic Data Centre 180 , the database focused on 2D materials is still under development. Several 2D materials databases have been attempted primarily based on a high-throughput approach coupled with first-principles and quantum theory to address this shortcoming. Besides, data-mining methods have also been applied to screen those bulk materials databases and identify possible 2D materials from 3D structures. These approaches can provide various information about 2D materials that have been or not been synthesized in experiments. For example, using a criterion to identify 2D materials, a database has been developed on the basis of experimental and DFT data 181 . This database provides information including exfoliation energy, bandgap, work function, and elastic-constant values of various 2D materials, and also predicted hundreds of new 2D materials for the experimental synthesis. DFT and many-body perturbation theory have been used to generate the Computational 2D Materials Database (C2DB) 182 . The C2DB contains structural, thermodynamic, elastic, electronic, magnetic, and optical information of about two thousand 2D materials in more than 30 crystal structures. A few hundred new monolayer structures Table 6. Common boundary conditions for the heat and mass transfer in macroscale models.

Navier-Stokes
No slip on walls u = 0 Normal inflow velocity a u = −U 0 ·n that were previously unknown, but potentially synthesizable were also identified. A data-mining method has also been used to screen the bulk materials from the open database and identified a few hundred 2D weakly bonded solid materials 183 . The band gaps and point groups that determine piezoelectric and nonlinear optical properties in these materials were also provided. A search in the database of bulk materials was performed and DFT calculations were used to explore these materials for possible exfoliation, where more than one thousand materials have been identified 14 . Information including electronic, vibrational, magnetic and topological properties of promising materials was also provided. Similarly, the crystal database was screened for layered motifs, and hundreds of stable layered materials were found that can be considered as candidates for the formation of 2D monolayers via exfoliation 184 .
Machine learning method Machine learning models trained using 2D material databases can be parameterized, validated, and therefore used for predictive exploration of novel 2D materials with desired properties. Machine learning methods can also aid the understanding of the complex correlations between structures and properties in 2D materials. A combination of machine learning with MD simulations and in-situ high-resolution transmission electron microscopy to explore TMDs has been performed 185 , where machine learning provided information for structural optimization and evolution of defects to help understand the structural transition in 2D TMDs. A practical method to explore hybrid 2D materials was developed by coupling machine learning with DFT calculations 186 . The structural and electronic properties of different hybrid 2D materials were provided and various parameters for vdW heterostructures were screened. A machine learning model with force-filed-inspired descriptors in material screening for complex systems has been introduced and used to discover exfoliated 2D-layered materials 187 . An artificial neural network for titanium dioxide systems was trained based on the a DFT calculated database, where a novel quasi 2D titanium dioxide structure was revealed 188 . Similarly, new 2D materials with high magnetic moments were found using a machine learning model trained by first-principles data 189 . Machine learning methods have also been used to aid the development of force fields for classical simulations of materials 190 . For instance, a force field for classical simulations of stanene was developed using a machine learning method trained by data sets from ab-initio results to calculate the mechanical and thermal properties of stanene 190 . 2D materials without bulk layered counterparts are also being discovered, using genetic algorithms 191,192 or particle swarm optimization, both using energy-based merit criteria and additional biases towards 2D sheets (see refs. 17,193,194 ).

OUTLOOK
This overview summarizes the state-of-the-art of modeling efforts on the growth of 2D-layered materials and outlines the advantages and limitations of computational models and methods at different length and time scales. The emphasis is on modeling, understanding, and predicting the thermodynamics and kinetics of 2D materials growth. The eventual goal is to develop computational tools and techniques to realize the "synthesis by design" of 2D materials. An ideal computational model should not only capture accurately the correct physics to allow robust designs but also be computationally efficient to be employed for active control of the growth process. The complex nature of the growth of 2D materials requires applications of a range of models, which operate at multiple length and temporal scales to capture the full range of the growth phenomena from macroscale flow and concentration profiles in the furnace to atomistic reactions leading to formation and growth of 2D materials at mesoscale. Modeling the full spectrum of mechanisms involved in the growth of 2D materials is challenging due to the limitations associated with each of the computational models at different length and temporal scales.
A high-fidelity prediction of 2D materials growth often requires an integration of methods applicable across different scales, ranging from electronic and atomistic, to mesoscale phenomenological, and large-scale continuum models. For example, chemical potentials of different phases and the migration energies of atomic species in different structural and chemical environment from electronic and atomistic models can be employed to parametrize mesoscale phase-field models for predicting the morphology and thus characteristics of as-grown 2D materials with the boundary conditions specified by the information obtained from large-scale continuum models such as temperature distribution in the growth chamber as presented in ref. 8 . Furthermore, technical obstacles to in-situ experimental monitoring of growth hinder the validation of computational predictions. The current mathematical and numerical models are computationally expensive and commonly do not capture simultaneously all the physics involved in a growth process of 2D materials, limiting their widespread adoption in active control of 2D materials growth. One promising way to circumvent this challenge is to use machine learning models trained using experimental and computational databases of the 2D materials structure and growth conditions such as the ones available at the 2D Crystal Consortium-Materials Innovation Platform. Finally, the models and approaches applied in understanding the growth of 2D materials can be adapted to the growth of other materials that use the same synthesis techniques, e.g., CVD growth of thin films.
As a future perspective, we can envision two main strategies for the development of computational models and their applications to the design and synthesis of 2D materials, which are long-term challenging goals of the field. First is the open-loop design approach, where the growth chamber design and synthesis conditions are determined by performing a series of simulations mimicking the thermophysical conditions and reaction and growth kinetics that would be expected in experiments. In this approach, the models must have high fidelity and a comprehensive sensitivity analysis should be performed to obtain a robust design. The second approach is the closed-loop design that the developed models will be utilized to determine the states of the system to be passed to a controller which adjusts the growth conditions in real-time. In this approach, the computational efficiency of the models is key to the success of the design.

DATA AVAILABILITY
All data that was obtained during this project are available from the authors.