A computational framework for guiding the MOCVD-growth of wafer-scale 2D materials

Reproducible wafer-scale growth of two-dimensional (2D) materials using the Chemical Vapor Deposition (CVD) process with precise control over their properties is challenging due to a lack of understanding of the growth mechanisms spanning over several length scales and sensitivity of the synthesis to subtle changes in growth conditions. A multiscale computational framework coupling Computational Fluid Dynamics (CFD), Phase-Field (PF), and reactive Molecular Dynamics (MD) was developed – called the CPM model – and experimentally verified. Correlation between theoretical predictions and thorough experimental measurements for a Metal-Organic CVD (MOCVD)-grown WSe2 model material revealed the full power of this computational approach. Large-area uniform 2D materials are synthesized via MOCVD, guided by computational analyses. The developed computational framework provides the foundation for guiding the synthesis of wafer-scale 2D materials with precise control over the coverage, morphology, and properties, a critical capability for fabricating electronic, optoelectronic, and quantum computing devices.


INTRODUCTION
Two-dimensional (2D) materials, including graphene and Transition Metal Dichalcogenides (TMD), have attracted increasing interest among the materials science community due to their unique thermomechanical 1,2 and optoelectronic [3][4][5][6] properties. They have applications in a wide range of industries, such as sensors 6 , actuators 7,8 , healthcare 5,9 , electronics 10,11 , and energy storage 12,13 . A key feature of 2D materials is the dependence of their properties on the morphology and size of the as-grown islands 14 . Thus, significant effort has been devoted to synthesizing uniform 2D materials with controlled morphology. The grand challenge for implementing 2D materials in commercial device technologies is the low-cost growth of wafer-scale single-crystal continuous films. Various synthesis methods have been developed, including mechanical exfoliation 15 , chemical exfoliation 16 , powder vaporization (PV)also called CVD [17][18][19] molecular beam epitaxy (MBE) [20][21][22][23][24] , and MOCVD techniques 25,26 . Mechanical and chemical exfoliation methods lack control over the shape and are generally accompanied by impurities and defects. The size of 2D domains grown by the MBE is commonly limited. The PV technique has been successfully utilized to synthesize domains as large as 100 μm in edge length [17][18][19] for demonstrating proof-ofconcept synthesis. However, it has limited scalability and almost no control over the precursor distribution and chemistry, which prevents its application in the industrial-scale synthesis of 2D materials. Among these methods, MOCVD is the most promising synthesis technique for the wafer-scale growth of 2D materials [26][27][28] . For example, large-area growth of monolayer and fewlayer films of WSe 2 , WS 2 , and MoS 2 using MOCVD has been realized [28][29][30][31][32][33] through exhaustive trial-and-error procedures.
Precise control of the growth morphology and coverage of atomically thin structures requires understanding the underlying physical mechanisms and realizing the correlation and significance of experimentally controllable parameters 34,35 . Several attempts have been made to understand these mechanisms and guide the synthesis process for the controlled growth of 2D materials. However, they are generally focused on one aspect of the growth and do not capture the complete picture of the synthesis process. We may refer to classical Wulff construction 36 , atomistic simulations [37][38][39][40][41] , and diffuse interface methods [42][43][44][45][46][47][48][49] as examples of these efforts. However, so far, precise control of the uniformity and coverage of as-grown 2D materials is considered a formidable challenge due to intrinsic complexities associated with the growth involving multiple physics across a wide range of length and time scales.
Our investigations show that the uniformity, morphology, and coverage of MOCVD-grown 2D materials are strong functions of precursor concentration and its gradient over the substrate, which are also influenced by the growth chamber configuration and flow characteristics. Here, we demonstrate the critical role of precursor concentration and its profile on the coverage and uniformity of layered materials grown by the MOCVD technique, both experimentally and theoretically (see Supporting Information for details of the CPM model and experimental synthesis). We develop a multiscale/multiphysics model of the growth, coupling a macroscale continuum model of the flow with reactive MD and mesoscale PF simulations to describe the growth kinetics involved in the MOCVD synthesis of layered materials across several scales. We consider WSe 2 , a critical 2D material for optoelectronic applications 50,51 , as our model material; we then validate our computational framework by comparing the simulation results with experimental data. We further perform a comprehensive study about the effect of different experimentally controllable growth parameters, such as inlet flow velocities, on the coverage 1 and uniformity of the MOCVD grown WSe 2 and potentially many other layered materials. Additionally, we demonstrate the role of substrate rotation during the growth of 2D materials on their uniformity and coverage.
We considered monolayer WSe 2 , a direct band-gap semiconductor 52 with native p-type conductivity 18,53 , which has applications in optoelectronic and photonic devices, as the prototype for this study.  Table 1) associated to different growth conditions. Second, we use reactive MD simulations 54 (Supplementary Methods Sec. 1.2 in Supporting Information) to calculate the energy of W-rich and S-rich edges at different W:Se mass ratios (see Fig. 1a) to find the equilibrium island shapes using the Wulff construction (Supplementary Figs. 5 and 6). Finally, the precursor concentration and orientation-dependent edge energies calculated in the first and second steps are passed to the PF model (see Supplementary Methods Sec. 1.3 in Supporting Information) to predict the morphology and coverage of the as-grown 2D materials.
The PF model simulates the growth process of WSe 2 in four consecutive stages that have been observed experimentally 25 : (I) random nucleation; (II) ripening; (III) lateral growth; (IV) and coalescence of the islands, forming atomically thin WSe 2 films. During stage I, a series of nuclei are randomly introduced in the PF solution domainthe substratefollowing the classic nucleation theory. Stage II consists of the growth of large nuclei at the expense of the smaller ones, i.e., ripening. In stage III, individual islands grow laterally due to the flow of precursors to the substrate. In the last stage, IV, the isolated islands will merge, forming the WSe 2 film. The morphology and uniformity of the asgrown 2D materials will be determined by competition between various thermodynamic driving forces, including the energy of different edges.

RESULTS AND DISCUSSION Experimental verification of CPM model
We detected the formation of hexagonal and full/truncated triangular-shaped islands, common patterns formed in the growth of 2D materials (see Supplementary Figs. 5 and 6). To verify the proposed computational framework and assess its robustness, when utilized to study real case growth of WSe 2 monolayers, we grow atomically thin WSe 2 on 2" c-plane sapphire substrates, using a unique MOCVD furnace ( Supplementary Fig. 1), allowing control over temperature, pressure, inlet velocities, and substrate rotation. We find that our simulated results are in precise agreement with experimental measurements. We confirm the applicability of the CPM framework at various common synthesis conditions for WSe 2 (see Supplementary Table 3). As anticipated, almost identical coverage and morphology were realized by the MOCVD growth of WSe 2 .

Effect of substrate rotation
To synthesize uniform large-area WSe 2 films, we analyze the CPM model simulation results; we revealedusing the CPM modelthat when the substrate rotates during the growth, higher uniformity in growth with large area coverage is obtained, Fig. 1b-g. Following this finding, we grow WSe 2 monolayers on a rotating substrate; this enables us to synthesize uniform large-area WSe 2 monolayers over two-inch wafers (Fig. 1h, i). Comparing the size distribution of monolayers with various orientations in the CPM model designated by different colors (Fig. 1f, g), we revealed the formation of a larger number of smaller monolayer WSe 2 islands on a rotating substrate with preferred growth morphology depending on their size. Furthermore, the maximum island size on the stationary substrate is more significant than on a rotating substrate. These results can be explained by the difference between the precursor concentrations between the two substrates, where a higher precursor concentration on the stationary substrate promotes the growth of larger islands. Furthermore, rotation of the substrate has a profound effect on the velocity profile, and thus, intermixing of precursor species and their velocity-assisted diffusion on the substrate, see Fig. 1c. It further introduces a downward component to the streamlines, increasing the adsorption of precursor species on the substrate. Our simulations also revealed the presence of higher precursor concentrations at the center of the rotating substrate and the front edge of the stationary one, potentially leading to the multilayer growth, which is consistent with our experimental observations (Fig. 1j). The large precursor concentration gradient at the front edge of the stationary substrate may promote agglomeration and out-of-plane growth via Mullins-Sekerka instability. A consistent concentration distribution obtained upon rotating the substrate further results in a more even W:Se ratio, leading to a more uniform growth morphology (Supplementary Figs. 7 and 8).

Nucleation and coverage
Upon detailed analysis, we realize that the concentration of precursor species and their gradient, under the joint effect of flow profile and temperature, determine the morphology and coverage of islands in MOCVD-grown monolayer materials. During the nucleation phase, isolated islands of WSe 2 are introduced on the substrate following an explicit nucleation algorithm and Poisson's seeding. A series of random positions are chosen to introduce nuclei with the critical radius if nucleation probability P ¼ 1 À expðÀI Á ΔtÞ is larger than a random number within [0,1], where I is the nucleation rate at a given time and position and a function of precursor concentration and temperature, and Δt is the integration time step. Precursor concentration and temperature profiles are determined by the furnace configuration and flow pattern from CFD simulations.
We assumed a uniform temperature distribution of 800 o C over the substrate, which allows us to study the precursor concentration distribution's role in nucleation and thus coverage of CVDgrown 2D materials. At the same time, precursors' concentration determines the equilibrium morphology of 2D domains and the uniformity of the synthesized monolayer films. Variation in concentration ratios of W(CO) 6 and H 2 Se can stabilize different morphologies as predicted by MD simulations and the Wulff construction (see Supplementary Figs. 4-6). Furthermore, interface energy between WSe 2 /Substrate and WSe 2 /gas phase is a function of precursor concentration, where variation in its values may promote the formation of 3D islands through the Volmer-Weber mechanism, hindering the growth of 2D monolayers. The Fig. 1 Effect of the chemistry and substrate rotation on the growth coverage and velocity profile. a The edge energy as a function of Sechemical potential varies between W-rich and Si-rich local environments for the five configured edge types; Δμ Se is a difference between μ Se (bulk) and μ Se . Steady-state concentration distribution of the precursors, W(CO) 6 (bottom) and H 2 Se (top), for the stationary substrate (b) and one rotating at 150 rpm (c), are shownarrow indicates the flow direction. The corresponding PF simulation for the stationary (d) and rotating substrate (e) are depicted along with the distribution analysis of the islands (f, g). Colors show different island orientations. The abscissa shows island area (designated A) while the ordinate shows island count in each bin (labeled with C). The equivalent experimental measurements are shown (h-i). j Variation in the thickness of the as-grown films across the diameter of the substrate. Error bars indicate standard deviation of the island thicknesses in each region. Velocity streamlines in the volume above the substrate for stationary substrate (k) and one that is rotating at 150 rpm (l). The colored volume shows the gas phase on the left half of the substrate, where the color map shows magnitude of the velocity. The substrate is shown as a green circle.
presence of precursors' concentration gradient normal to the substrate will also promote the formation of out-of-plane growth via the Mullins-Sekerka mechanism 34,35 . Thus, the highest uniformity of 2D films can be achieved when the precursors' concentration is homogeneously distributed over the substrate with the minimum gradient in the direction normal to the substrate.

Effect of temperature
Considering the aforementioned micro-processes determining the growth, we hypothesize that the precursors' concentration and their distribution are vital parameters determining the coverage and uniformity of synthesized 2D materials, and they depend on several factors such as flow of precursors and their ratios, location of gas inlets, growth temperature, and growth pressure. Adjusting the growth temperature has a multifaceted effect on precursor concentration and distribution profiles; increasing temperature reduces the density of carrier gas, dominating the buoyancy forces downstream (see Fig. 2a-c), as well as raising heat capacity, heat conductivity, and dynamic viscosity (see Fig. S3 Supplementary Information). This effect is significant for stationary substrates, where a reduction in growth coverage and thickness of 2D materials downstream is revealed (Fig. 2d-g). Lowering the growth temperature also reduces the nucleation energy, which leads to an increase in the nucleation rate, but increases thermodynamic driving forces for the growth of solid monolayers and thus reduces the radius of the critical nucleus. It will further increase the thickness of the as-grown films while forming smaller nuclei (Fig. 2d-g). Reducing the growth temperature further reduces the surface diffusion length and desorption of reactive adatomsparticularly Sefrom the substrate 55 , reducing the island size. It also reduces the diffusion that is dominated by flow velocity.

Effect of pressure
In addition to growth temperature, growth pressure is another factor affecting the concentration and distribution of precursors, as it influences density, the diffusion coefficient of gas-phase precursors, and the range of stability of different morphologies. The density of carrier gas is a linear function of pressure (ρ ¼ P Á cte=T). The diffusion coefficient is inversely proportional to pressure and can be calculated using Fuller's semi-empirical method 56,57 , i.e., D ¼ D p T 1:75 ð Þ=P, where D p is diffusivity or pressure-independent diffusion coefficient. Reducing pressure results in a higher diffusion coefficient, (D / 1=P), leading to a higher concentration of precursors, while it reduces density (ρ / P) and reduces precursor concentration. Thus, it acts as a double-edged sword, indicating the presence of an optimized pressure that maximizes the concentration of precursors and their coverage. To verify this conclusion, we synthesize WSe 2 monolayers at three different growth pressures −50, 200, and 500 Torrwhile other growth conditions are held constant. Optical images and AFM scans reveal degradation in uniformity of synthesized monolayers as growth pressure deviated from the optimum 200 Torr pressure. Furthermore, lowering the growth pressure reduces the equilibrium concentration, c eq , and increases the nucleation rate, which is consistent with our simulations and experimental observations (Fig. 3). We simulated the growth within a range of D ¼ 1 À 10 3 ½ 10 À9 m 2 s À1 to understand the effect of pressure on precursor concentration and its distribution profile while holding ρ constant to investigate the role of the diffusional part of the pressure. The simulation results indicate improved growth coverage and formation of fewer yet larger islands upon reducing pressure (or increasing diffusion coefficient), neglecting changes in carrier gas density that are consistent with our experiments.

Effect of inlet velocity
Additionally, using the CPM model allows us to study the effect of carrier gas velocity on the growth of WSe 2 films. Keeping the growth conditions fixed and varying the inlet carrier gas velocity, a more extensive growth coverage is achieved when we increase the carrier gas velocity up to a certain point beyond which changes in the growth coverage are minute. A higher inlet flow rate improves velocity-assisted diffusion of precursor species, limited by the thickness of the boundary layer, δ, forming on the substrate that is inversely proportional to the inlet velocity, V 0 , and can be established from the Blasius solution for laminar flow as a   function of distance from the susceptor x, i.e., δ / ffiffiffiffiffiffiffiffiffi ffi x=V 0 p 58 . Decreasing the boundary layer thickness reduces the diffusion distance that atoms need to travel to reach the substrate. Thus, for a specific diffusion and growth time, a larger amount of precursor material reaches the substrate upon increasing inlet velocity, improving the growth coverage. The boundary layer thickness also increases downstream, reducing the concentration of the precursor reaching the substrate and, thus, the coverage and uniformity of the 2D film. Although an increase in the inlet velocity leads to higher precursor concentrationtherefore, larger area coverage and multilayer growthit adversely affects the adsorption of adatoms and increases their desorption, increasing the chance of defect and void formation. Similarly, an MOCVDgrowth experiment was performed for different inlet carrier gas velocities with a growth temperature, T = 800 o C, at a growth pressure of 200 Torr, revealing an increase in the growth coverage and multilayer growth upon increasing the inlet carrier gas velocities, perfectly matching predicted simulation results (see Fig.  4).
In this study, we developed a multiscale/multiphysics modeling framework for the growth of 2D materials using the CVD technique, called the CPM method, exploring some of the delicate micromechanisms governing the synthesis of CVD-grown 2D materials and their complex relationship with experimentally controllable parameters. Comparing the CPM simulation results with experimental measurements, we identified two key factors influencing the coverage and growth uniformity of the 2D islands, i.e., concentration and concentration gradient of precursor. We performed a comprehensive experimental study on the effect of various growth control parameters. We demonstrated a perfect match between the developed CPM model and the experimental observations for MOCVD synthesis of WSe 2 and potentially other 2D materials. Finally, we show significant improvement in uniformity and coverage of the synthesized atomically thin films of WSe 2 by rotating the substrate during growth, which has a major impact on the wafer-scale synthesis of 2D materials and thus fabrication of high-performance optoelectronic devices, particularly with applications in flexible electronics.
This novel computational framework for the analysis and design of CVD-grown 2D materials enlightens physical mechanisms governing the growth of these materials across several scales. At the atomic scale, using reactive MD simulations, it provides the relationship between the ratio of precursor concentrations and stable morphology of 2D materials; At the macroscale, it connects experimentally controllable parameters with concertation of precursor and its gradient over the substrate; It combines the information from both nanoscale and macroscale simulations to predict the coverage and uniformity of as-grown 2D materials. This computational framework provides a unique alternative to exhaustive trial-and-error experimentations and a powerful tool to develop and optimize the synthesis of new 2D materials. It can further serve as an observer for controllers of the growth process, providing the feedback loop capability, thus, precise control over the growth process, which opens new routes to design and fabricate the next generation of nanodevices for application in quantum computing and artificial intelligence.

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