Introduction

Controllable synthesis of 2D materials with desirable morphology is challenging and yet is an important prerequisite for realizing their unique functional properties and corresponding applications. Vapor-based growth methods have been successfully applied to synthesizing various 2D materials. For example, chemical vapor deposition (CVD) has been applied to the growth of monolayer1,2,3 and multilayers4 of MoS2 and related heterostructures;5,6 powder vaporization (PV) has also been applied to synthesize both vertical and planar MoS2 crystals.7 These growth methods provide several adjustable control parameters, including powder/vapor species, gas flow rates, substrate materials and positions, and temperatures, to fabricate 2D materials with different morphologies. However, in practice, even subtle changes of these parameters can lead to significant differences in growth qualities of 2D materials. The growth process will be particularly sensitive to changes of growth parameters for materials with a narrow region of stability of desired intermediate kinetic phases. As an example, we may refer to change in the MoS2 growth mode simply by changing the orientation of the substrate with respect to the trajectory of the carrier gas flow.8,9 To further improve the controllability and growth quality, as well as to facilitate the comparisons of results from different CVD chambers, it is critical to understand the underlying mechanisms controlling the growth morphology during fabrication and the effect of each control parameter on the growth morphologies of 2D materials. Since the vapor-based synthesis methods involve complicated interactions among transport phenomena including laminar gas flows, mass transfer and heat conduction, various chemical reactions involving multiple transient species, as well as crystal growth, such comprehensive understandings have not yet been achieved. As a result, the optimization of these control parameters and the control of growth qualities of 2D materials still largely relies on the trial-and-error-based experimental studies.

On the other hand, theoretical calculations and computational simulations are playing increasingly important roles in the understanding of growth mechanisms during vapor-based synthesis of 2D materials. At the atomic scale, first-principles calculations have shown that the anisotropic edge energies are critical to the island morphology evolution of MoS210 and h-BN.11 The Reactive Force Field (ReaxFF) Method12 has been applied to revealing the chemical reaction mechanisms and pathways of forming the 2D materials from the gas/powder precursors,13,14,15 and predicting the growth morphologies of a few islands. The kinetic Monte Carlo (KMC) method has also been coupled with first-principles calculations to predict the island growth morphology as a function of flux and precursor stoichiometry.8,16,17 At the mesoscale, the phase-field approach18,19,20,21 has been applied to investigating the effect of deposition flux and diffusivity on the island morphologies20,22,23 as well as the interactions among different islands.24 At the macroscale, finite element methods (FEM) focused on transport phenomena have been employed to predict the distribution of velocities, concentrations and concentration gradients of precursors during synthesis,9,25 which have provided useful guidance for experimental synthesis by correlating the concentration gradient distribution with the growth modes of 2D materials. However, since most of these existing investigations are focused on different individual aspects of the entire synthesis process, i.e., atomistic calculations on thermodynamic properties, reaction mechanisms and pathways, mesoscale models on crystal growth kinetics, and FEM on transport phenomena, the coupling and interactions among these different aspects, which are required for developing the direct connection between the CVD control parameters and 2D island growth morphologies, are highly demanding.

A comprehensive understanding of the morphological evolution of 2D materials during a PV synthesis process requires the integration of theoretical and/or computational models at different length scales. For example, the pioneering work by Govind Rajan et al.8 formulated and validated a generalized mechanistic model for CVD of 2D materials with the linkage to CVD parameters, which provide useful guidance for controlled synthesis of such materials. Enlighten by this work, we propose a multiscale modeling approach, integrating an FEM growth chamber model and a crystal growth phase-field model (PFM), to investigate the effect of precursor concentration, its spreading, and deposition rate on the island morphology evolution and distribution of 2D materials on a substrate. While Govind Rajan’s work8 captures the fundamental thermodynamic and kinetic growth mechanisms of 2D materials, the current work can capture the island growth modes, including growth instabilities and island interactions, and can directly simulate the island size and morphology distributions on the whole substrate, which are difficult for most of the other approaches. The proposed model considers the coupling between fluid flow, diffusion, heat transfer, and buoyancy forces and sensitivity of the growth to changes of these parameters. We specifically apply the multiscale model to the CVD growth of MoS2. Simulations indicate that the 2D islands initially nucleate and grow at the regions with high precursor concentration and deposition rate, which affects the subsequent nucleation and growth of other regions including the island size, shape, and inter-island spacing. It is shown that the precursor concentration distribution and deposition rate play critical roles in the morphology and spatial distribution of 2D islands on a substrate. Experimental measurements are subsequently performed to validate the simulation results indicating an excellent match between theory, simulation, and experiments.

Results and discussion

The coupled system of equations for the large-scale heat and mass transport, Eqs. (S1)–(S3), are solved simultaneously to give the distributions of temperature, precursor concentration, and velocity field within the growth chamber. Schematic of the experimental setup and corresponding computer model are presented in Fig. 1a,b. The results of normalized precursor concentration c* and velocity field of the carrier gas are plotted in Fig. 1c. The precursor concentration on the substrate calculated using the FEM simulations is plotted in Fig. 1d, which shows an excellent agreement with the experimental observations shown in Fig. 1e. Multilayer growth was observed in the areas of high concentration. The carrier gas velocity is larger at the entrance of the crucible, which leads to the maximum spreading of precursor. Moving downstream and particularly close to the substrate, where no slip boundary conditions are enforced, the velocity of carrier gas reduces and precursor spreads mainly via diffusion. There is also a change in the flow direction close to the front and end points of the substrate due to the formation of a boundary layer, where the component of the flow velocity normal to the substrate becomes comparable to the component parallel to the plane of the substrate, as shown in Fig. 2. It results in the formation of a precursor concentration gradient (\(\nabla c^ \ast\)) normal to the substrate, which may subsequently promote the out-of-plane growth of 2D materials and formation of vertically grown 2D materials through the Mullins–Sekerka mechanism.26,27 Our results indicate that the out-of-plane component of the \(\nabla c^ \ast\) is significant in the region close to the front edge of the substrate, and it reduces drastically as we move further downstream. On the other hand, the in-plane component of the \(\nabla c^ \ast\) is negligible in the front region of the substrate while it dominates toward the center and end sections of the substrate, indicating the larger driving force for the growth of in-plane monolayer materials.

Fig. 1
figure 1

Simulation results for the large-scale model of the growth chamber. a Schematic of the experimental setup; b the corresponding 3D computer model and the dimensions; c distribution of normalized concentration, c*, of precursor (volume data) and velocity profile (stream lines) within the volume enclosed between the crucible and the substrate (only half of the volume represented here for clarity); d normalized concentration of precursor over the substrate; and e experimentally observed precursor deposition on the substrate where the color variation shows the concentration of deposited material. For d and e only half of the substrate is shown for comparison between simulation and experiment, which shows an excellent agreement

Fig. 2
figure 2

Simulation results for \(\nabla {\vec{\boldsymbol c}}^ \ast\). The results of FEM simulations indicate a large gradient in concentration of precursor materials on the bottom side of the substrate normal to its plane a, specifically on the side closer to the upstream. This can result in out-of-plane growth of monolayer materials. Large in-plane precursor concentration gradients can be seen in the area further away from the front edge of the substrate b, where the out-of-plane concentration gradient is also negligible. This is the area with potential growth of large-area 2D materials

The temperature profile within the gas-phase is shown in Fig. S7 in Supplementary Information, indicating temperature gradual increase closer to the heating region. The gas phase temperature reaches a maximum at Tmax = 1050 K followed by a slight drop, about 4°, due to the heat transfer with crucible, and a subsequent increase in temperature to Tmax. It reduces again as it exits the heated region.

We performed phase-field simulations of the island morphology evolution under the non-uniform distribution of concentration and deposition rates within a larger region on the substrate with an angle-dependent anisotropy function. For simplicity, we only consider two different island orientations 60° to each other (i.e., vertices pointing up and down for triangles). With the initial concentration cini = g(r)ceq fitted from FEM simulations (Fig. 3a), we set the initial saturation uini = (ciniceq)Ω, where Ω is the molar volume of MoS2; we also assume the deposition rate follows the same spatial distribution with the initial concentration. Langevin random noises are used to trigger the nucleation of the islands.28 The 2D simulation region is discretized into 3600 × 3600 grid points, representing a 3.6 mm × 3.6 mm region in real size. We assume the simulation region has the same concentration distribution with that from FEM simulations (Fig. 3a). The message passing interface (MPI) parallelization is used for the finite difference scheme to facilitate the simulation using up to 720 computer cores. The simulation results are shown and compared with experimental observations in Fig. 3 (800 °C) and Fig. S15 (725 °C).

Fig. 3
figure 3

Simulation and experimental results for MoS2 island morphology distribution at 800 °C. a Initial precursor concentration distribution from FEM. b Phase-field simulation results for MoS2 morphology on substrate at t = 0.25 h, the enlarged images of i–iii boxes are shown in (df). c Experimental results for deposition of MoS2 after deposition for t = 0.25 h, the enlarged images of 1–3 boxes are shown in (gi). df Enlarged images for the selected microscopic regions marked by the boxes in (b). gi Experimental results on morphology of as-grown MoS2 for selected regions in c. Comparing (df) with (gi) shows an excellent match between simulation and experiment

The island morphology dissemination is largely affected by the non-uniform distribution of concentration and deposition rates, Fig. 3b. The MoS2 islands first nucleate in the regions with high concentration and deposition rates with dodecagonal or hexagonal shapes (Fig. 3d); the nuclei then sequentially develop towards the regions with low concentration and deposition rates, forming quasi-hexagonal and triangular islands (Fig. 3d,e). The first-principles predictions of island morphology change due to local chemical potentials, i.e., the S-rich atmosphere (or low concentration of MoO3 precursor) favors the triangular shape while the Mo-rich atmosphere (or high concentration of MoO3 precursor) favors the more isotropic dodecagonal island shapes8 as shown in Fig. S9 in Supplementary Information, are well captured in the phase-field simulations, as illustrated by the corresponding island morphologies in Fig. 3d–f, respectively, which are enlarged images of the selected regions i–iii in Fig. 3b.

The developed model is further verified with the experimental results by synthesizing monolayers of MoS2 using the CVD technique. The variation in morphology of the synthesized MoS2 monolayers is shown in Fig. 3c,g–i where ~10 μm-sized hexagonal MoS2 (Fig. 3g) form at the leading edge of the substrate (spot #1 in Fig. 3c) changes to quasi-hexagonal, 3(h), (spot #2 in Fig. 3c) and then to triangles, 3(i), (spot #3 in Fig. 3c). This observation is in excellent agreement with simulation results shown in Fig. 3d–f. A more quantitative comparison between simulation and experimental results in terms of island sizes, number of edges and coverages for the three selected regions (spots #1–3) is listed in Table S4. The predicted average island sizes and number of edges per island are quite consistent with experimental results.

In addition to the shape variations, other significant morphological features of MoS2 islands are also captured in our model. The coverages in Fig. 3h,i are lower than that in Fig. 3g, similar to the phase-field simulation results (Fig. 4) at low precursor concentration regions, which is likely to be due to the later nucleation and lower deposition rate. The triangles in Fig. 3i are less equilateral with concaved edges, which is generally consistent with the phase-field simulations in Fig. 5. Therefore, the experimental observations qualitatively validate the multi-scale simulations.

Fig. 4
figure 4

Phase-field simulation of MoS2 island morphology distribution and statistical analysis at different simulation time steps at 800 °C. ac Phase-field simulation results at t = 0.1, 0.5, and 5 h, respectively

Fig. 5
figure 5

Development of growth instabilities of MoS2 islands from phase-field predictions at 800 °C. a1a4 Enlarged images for the region marked by the blue box in Fig. 4b at t = 0.5, 0.75, 1, and 2 h, respectively. b1b4 Enlarged images for the region marked by the red box in Fig. 4b at t = 0.5, 0.75, 1, and 2 h, respectively. The white arrows in a2 and b2 indicate the growth directions of the triangle tips

To more clearly show the nucleation (deposition), growth, and coalesce processes of the MoS2 islands, we plot the simulated MoS2 island morphology distributions in Fig. 4a–c at three different time steps: t = 0.1, 0.5, and 5 h. We also perform statistical analysis for the island coverage, as shown in Fig. S12(a-c). The whole simulation region is uniformly divided into 10 × 10 regions, identified by the corresponding x- and y-indices. The coverages are counted within each region. From Fig. 4a,b, one prominent morphology feature is the formation of island chains and arrays with increasing spacing among each other, which develop sequentially with time. This is because there is a critical saturation value uc = 0.251 for the spontaneous nucleation of MoS2 islands under infinitesimal perturbations, i.e., the local free energy \(f_{local} = \left( {\phi ^2 - 1} \right)^2 - \lambda u\left( {\phi - \phi ^3{\mathrm{/}}3} \right)\) becomes unstable \(\left( {\frac{{\partial ^2f_{local}}}{{\partial \phi ^2}} < 0} \right)\) in the gas phase (ϕ = −1) when u is larger than uc. Moreover, due to the non-uniform distributions of precursor concentration and deposition rates, the time required for the local u to reach uc is different at different positions of the substrate. As a result, the islands nucleate and grow along the iso-concentration lines, forming island chains. As the time required for u to reach uc increases nonlinearly with the decrease of initial precursor concentration and deposition rates at the upper-right and lower-right corners on the substrate, the spacing among the island chains are increasing along the concentration gradient. These features can be explained by the 1-D simulations, Fig. S8 in Supplementary Information.

The growth and coarsening behaviors of the islands are also affected by the non-uniform distribution of precursor concentration and deposition rates. As shown in Fig. 4b,c, within the center-left region of the substrate, due to the high precursor concentration and deposition rates, the nucleation number density is high, and the inter-island spacing is rather limited, and the growth of the islands is dominated by island competition, coalescence, merging, and coarsening. On the other hand, the upper-right and lower-right corners on the substrate have lower nucleation number density and larger inter-island spacing, and the island growth behavior is featured by: (i) coarsening along the iso-concentration line; and (ii) growth and development of growth instabilities along the concentration gradient (perpendicular to the iso-concentration line). Especially, the growth instabilities are evident by the pointy-triangular morphologies along the outer island chain in Fig. 4b, pointing to the concentration gradient direction. According to the Mullins–Sekerka theory,26,27 the growth instabilities become active with the accumulation of solutes ahead of the island edges, leading to the formation and amplification of tips at the initially flat edges. The initiation of this instability, or, the accumulation of solutes to some critical value, originates from the different solute diffusivities of the two phases, and can be promoted if the edge diffusivity is poor or the deposition flux is large. In Fig. 4b,c, however, we find the inter-island spacing can be another influencing factor for the development of growth instability, especially in this complicated CVD process involving simultaneous nucleation, growth and interactions of islands. In the center-left region of the substrate with high nucleation number density and low inter-island spacing, the islands would not have enough chance to grow to the critical size with sufficient solute accumulation to trigger the growth instability; meanwhile, the continuous growth and deposition further avoid the contact of island edges with vapor phase in the 2D plane in this region while promote the coalescence of different islands, which further kills the growth instability. On the other hand, at the upper-right and lower-right region of the substrate with low nucleation number density, low inter-island spacing along the iso-concentration line, and high inter-island spacing along the concentration gradient, the island edges have enough contact with the vapor phase and sufficient space to grow larger along the concentration gradient, causing Mullins–Sekerka-type growth instabilities when the islands grow to critical sizes with sufficient solutes accumulated in front of the island edges. Preferentially, the pointy-triangular morphologies develop at the triangle vortices pointing towards the concentration gradient direction due to the larger growth spaces. The situation can become more complicated when the multi-island interactions are involved.

We select two microscopic regions on the substrate, marked by the blue and red boxes, respectively, to track the growth instabilities under island interactions, as shown in Fig. 5a1–a4,b1–b4, respectively. For the blue-boxed region, the concentration gradient is pointing to the lower-right direction as indicated by the arrows in Fig. 5a2. Therefore, the lower-right vertices of the upward triangle array first elongate towards the concentration gradient direction, as indicated by the white arrows in Fig. 5a2, causing the island shapes to deviate from ideal equilateral triangles from the Wulff construction predictions (Fig. 5a2). The two elongated edges then become curved rather than remaining straight, developing kinks at the positions with negative curvature (Fig. 5a2,a3), leading to growth instabilities. The kinks can continuously develop at the elongated island front as long as there is enough space (Fig. 5a4). Fig. 5b1–b4 shows another example, i.e., the red-boxed region in Fig. 4b. In this case, the concentration gradient is pointing downwards so the lower vertices of the downward triangles elongate (Fig. 5b1), and the growth instabilities develop for the triangle on the left. The triangle island in the middle elongates both downwards and to the right, as indicated by the white arrow in Fig. 5b2, but no significant growth instabilities (i.e., development of kinks) are observed due to the limited growth space confined by the islands on its right and at the bottom (Fig. 5b3,b4). Similarly, due to the islands developed at the bottom, the growth instability of the island on the left is killed (Fig. 5b4), and no further kinks are developed before merging with the island at the bottom. The island morphology evolutions in these two selected regions show typical growth instability features under the effect of non-uniform distributions of precursor concentration and deposition rates during CVD growth.

The islands can further grow into a polycrystal thin film as CVD process continues for an enough long time. As shown in Fig. 4c, the center-left region of the substrate will first form thin films, with the island coverage close to 100% (Fig. 4f), leaving defects such as grain boundaries between islands with different crystal orientations. Meanwhile, the upper-right and lower-right corners on the substrate will continue to grow through the combination of island planar growth, island growth instabilities along preferred directions (along concentration gradient) and nucleation of new islands, until the full coverage.

A hierarchical multiscale model is introduced here to model the growth morphology during CVD growth of monolayer MoS2. The large-scale coupled mass and heat transport system of equations are numerically solved to determine thermodynamic state of the precursor and its concentration within the growth chamber. The calculated concentration distributions are subsequently transferred to the mesoscale PFM where the nucleation and growth of 2D materials on the substrate are investigated. The PFM incorporates the anisotropic edge energies from first-principles calculations, and considers the effect of non-uniform distributions of precursor concentration and deposition rates. The developed model is validated experimentally by growing monolayer MoS2 using CVD technique at two different temperatures, 800 and 725 °C. The excellent agreement between experimental observations and numerical simulations proved the predictive capability of the developed model, including variation in the morphology, distribution, and size of synthesized 2D materials with furnace parameters. We further studied a vertical CVD furnace with a heated susceptor (Fig. S12 in Supplemental Information), where we demonstrated formation of eddy currents when there is a gap between the substrate and the susceptor. Formation of these currents will disturb precursor concentration and subsequently growth of the 2D materials (Fig. S13 in Supplemental Information). The current framework can be used to create growth diagrams specifying the range of growth parameters that result in the formation of 2D materials. It also can be used to train neural networks to harvest the capability of artificial intelligence techniques29,30 for solving the inverse problem of finding proper growth conditions based on the desired morphology and characteristics of the 2D materials to be grown.

Methods

Computational model

The multiscale model (Sec. 1 of Supplementary Information) integrates an FEM growth chamber model for predicting the continuum level transport phenomena including the laminar gas flow, heat transfer and flow-assisted diffusion, and a PFM crystal growth model for simulating the mesoscale 2D island morphology evolution during CVD growth processes. We assumed that the macroscopic transport phenomena is much less affected by the more localized island growth because the transport phenomena take place throughout the CVD chamber while the island growth only happens on a substrate inside the chamber. Therefore, in FEM we neglect the island growth and only consider transport phenomena. The calculated concentration distribution from FEM is imported into PFM to simulate the island morphological evolution on the substrate. We assumed a vapor to solid deposition process from supersaturated precursors.

Experimental method

The CVD growth of MoS2 was conducted in a horizontal hot zone furnace, where the diameter of the tube reactor is 40 mm. The monolayer MoS2 films were grown by the PV technique. The reactor set-up is shown in Fig. 1a. In this study, 300 nm SiO2 on Si (University Wafer) was used for substrate. Before the growth, the reactor remained as vacuum (~18 mTorr) for 10 min in order to remove the moisture, etc. The growth was conducted at 800 °C (and 725 °C) for 15 min with 2 mg MoO3 (Sigma Aldrich) located 0.5 cm upstream of the face-down SiO2 (University Wafer) substrate (Fig. 1a). 200 mg sulfur powder (99.95%, Sigma Aldrich) was put ~15 cm (6 in.) upstream of the substrate. During the growth, the MoO3 powder and substrate was co-heated to 800 °C (and 725 °C) for 15 min and S powder was heated by a separated heat tape at 130 °C. 100 sccm Ar was used as the carrier gas. The pressure of the growth was maintained at 710 Torr. To analyze the morphology of MoS2 domains at different locations with different precursor concentrations on the substrate, the synthetic MoS2 was characterized by optical microscopy. Optical microscopy images were taken every 0.4 cm (Fig. 3g,h at 800 °C and Fig. S15g-h at 725 °C). Atomic force microscopy measurements were carried to ensure the formation of single layer MoS2 where the results are presented in Fig. S14 of Supporting Information.

Data availability

All data that was obtained during this project is available from the authors.