Dynamic Growth/Etching Model for the Synthesis of Two-Dimensional Transition Metal Dichalcogenides via Chemical Vapour Deposition

The preparation of two-dimensional transition metal dichalcogenides on an industrially relevant scale will rely heavily on bottom-up methods such as chemical vapour deposition. In order to obtain sufficiently large quantities of high-quality material, a knowledge-based optimization strategy for the synthesis process must be developed. A major problem that has not yet been considered is the degradation of materials by etching during synthesis due to the high growth temperatures. To address this problem, we introduce a mathematical model that accounts for both growth and, for the first time, etching to describe the synthesis of two-dimensional transition metal dichalcogenides. We consider several experimental observations that lead to a differential equation based on several terms corresponding to different supply mechanisms, describing the time-dependent change in flake size. By solving this equation and fitting two independently obtained experimental data sets, we find that the flake area is the leading term in our model. We show that the differential equation can be solved analytically when only this term is considered, and that this solution provides a general description of complex growth and shrinkage phenomena. Physically, the dominance suggests that the supply of material via the flake itself contributes most to its net growth. This finding also implies a predominant interplay between insertion and release of atoms and their motion in the form of a highly dynamic process within the flake. In contrast to previous assumptions, we show that the flake edges do not play an important role in the actual size change of the two-dimensional transition metal dichalcogenide flakes during chemical vapour deposition.

By now, for the fabrication of 2D TMDCs in large scales, an access via top-down methods exists [16][17][18]. However, in particular bottom-up methods such as chemical vapour deposition (CVD) are promising due to their potential compatibility with processes for thin film fabrication established in the semiconductor industry. With bottom-up methods, the desired materials are formed by self-assembly of the corresponding precursor atoms. Requirements for this are, for instance, suitable high temperatures and the supply of sufficient precursor material.
After Lee et al. first reported the successful growth of the 2D TMDC molybdenum disulphide (MoS 2 ) with CVD in 2012 [19], various 2D TMDCs have been synthesized on different substrates, as well as on other van der Waals and 2D materials . One of the most rudimentary realization of CVD for 2D TMDCs is based on the use of two solid precursor sources containing either the chalcogen (e.g. elemental sulphur powder) or the transition metal (e.g. transition metal oxides or chlorides). The similarities of the process design for many 2D TMDC species and of their resulting morphology imply identical atomic kinetics during CVD.
In the last decade, deeper insights into the growth mechanisms of 2D TMDCs (mostly 2D MoS 2 ) have been collected by refining the process systems and recipes as well as by developing models. For example, the initial nucleation has been studied extensively [24,30,31,33,34,[36][37][38]. As a result, concepts for controlling nucleation were presented, e.g. by the use of seeding promoters [24] or artificial defects in the substrate [36,37]. Furthermore, different growth rates for the 2D TMDC crystal facets have been identified as the reason for the typically (equilateral) triangular flake shape from CVD. These growth rates differ in their dependencies of the ratio of the precursor atom species (transition metal or chalcogen), allowing potentially to control of the edge termination and even the shape of the resulting flakes by the precursor atom concentrations in the gas phase [27,34]. Recently, concepts have been proposed in order to describe the dynamics and the stability of the orientation of 2D TMDC flakes growing on crystalline van der Waals materials [39]. Here, the flake orientation might be controlled by substrate defect engineering [37,41].
However, one phenomenon has yet only been insufficiently elucidated. It is experimentally found that 2D TMDC flakes first grow and then shrink again as the process duration increases, see the publication by Chen et al. [32] or Fig. 3. Because degradation of 2D TMDCs is also facilitated by increased temperatures [10,[42][43][44][45][46][47][48][49][50][51] [and see SI 1 ( Fig. S1)], a dynamic process between growth and etching during a CVD process must obviously exists. Chen et al. assume -without providing any theoretical model -that insertion and release of atoms only takes place at the edges of the grown 2D TMDC flake (labelled and discussed below as growth rate G 1D and etching rate E 1D , respectively). This assumption seems to be straightforward and intuitive for describing synthesis of 2D TMDCs as the basal plane is chemically rather inert, while edges represent the active sites, at least at low temperatures. In addition, Chen et al. explain their observations by the absence of adsorbed material on the basal plane of growing flakes. That material supply takes place exclusively via the substrate is also implied by Wang et al. [27]. But are these assumptions really adequate to describe the growth process?
In this paper, we introduce an advanced concept explaining the experimentally observed fact, that the size of 2D TMDC flakes at first increases and then decreases again during CVD. Our mathematical concept is based on considerations taking material supply and transport into account as well as its change over time and is complemented by thermal degradation/etching mechanisms concluded from experiments. Finally, we apply the resulting equations to data for MoS 2 on sapphire by Chen et al. [32] and to our own data for tungsten disulphide (WS 2 ) on sapphire. Contrary to intuitive assumption we find, both mathematically and by the best fits on both data sets, that the change in area of the 2D TMDC flakes is largely proportional to the flake area itself.

Results and Discussion
To describe the experimental results correctly, our model must account for growth and shrinking mechanisms which will depend on material transport and supply. Within the model, the time-dependent change of the area of a two-dimensional TMDC flake dA dt depends on its own, current size (defined by its area A or its edge length L, respectively). In this context, we distinguish between rates that result in an increase in the flake size and those that lead to a decrease. Accordingly, they are called growth rate G and etching rate E, respectively. At first, we will discuss these rates in order to develop our basic differential equation shown later in Eq. (1).
The growth rate of a single 2D TMDC flake G depends on the size of the total supply area from which precursor material can agglomerate to form a new flake (nucleation) or diffuse to appropriate sites in pre-existing flakes in order to increase its size (direct growth) or to compensate etching (indirect growth). The supply area itself is related to the flake size -or simplified: larger flakes can "catch" more precursor material. Because, corresponding to Fig. 1 a), the total supply area is composed of different areas with different dependencies on the flake size (and on the different surfaces, see Fig. 1 b)), we distinguish between three supply areas and thus also separate the growth rate into three individual growth rates G nD (n = 0, 1, 2). The illustrations for various flake sizes in Fig. 1 a.i)-a.iv) clarify the relationship between the separate supply area sizes and flake sizes. In detail, the orange supply area remains constant (independent of spatial dimensions of the flake: 0D), the green supply area is proportional to the edge length of the flake (dependent on one spatial dimension: 1D), and the blue supply area corresponds to the flake area itself (dependent on two spatial dimensions: 2D). These three supply areas leading to the corresponding growth rates G nD will be discussed individually in the following paragraphs.
At the very beginning, when no flake is present, only the flake size independent 0D supply area contributes to the growth or -in this very particular case -to the formation of the first 2D TMDC flake by its corresponding growth rate G 0D . This case is visualized in Fig. 1 a.i), in which only the circle-shaped orange supply area is present. The size of this area depends on the adsorption and desorption rates as well as on the diffusion constant of precursor material on the substrate and thus mirrors the probability of the event of randomly agglomerating precursor atoms adsorbed on the substrate surface.
Once a 2D TMDC flake is formed (corresponding to the blue triangle in Fig. 1 a.ii)a.iv)), the additional 1D and 2D supply areas emerge. The background of the 1D supply area (green) and its corresponding growth rate G 1D is, that precursor material adsorbs on the substrate surface near the flake, where it can diffuse to the flake and be built in at its edge (1D line of reactive sites) before desorbing. As the distance from which  The size of the supply areas are either independent of spatial dimensions (orange, 0D), or dependent on one (green, 1D) or two dimensions (blue, 2D), respectively. b) Precursor atoms either diffuse on the surface of the substrate (0D, 1D) or of the flake (2D) to be inserted at the flake edge. The constant c micros. summarises microscopic constants such as the adsorption and desorption rate, the diffusion constant as well as the reaction probability and differs for the supply areas. c) The schematic of a basic CVD system illustrates the dependency of the adsorbed material on the amount of source material M (t) and its specific vaporisation rate v. The constant c macros. takes into account losses caused due to the transport from the source to the target substrate.
precursor material can diffuse to the flake edges is constant, the 1D supply area increases proportionally to the edge length L. Given that precursor material from this supply area is built in only at the edges of the flake, its contribution to the increase of the flake area A is direct (direct growth). Therefore, this mechanism might be the most intuitive one.
Because the third supply area (blue) is the flake area A itself, in this case -in contrast to the previously described 0D and 1D supply areas -the precursor material is not adsorbed on the substrate, but on the already grown 2D TMDC flake, see Fig. 1 b). It is very likely, that adsorption/desorption rates as well as diffusion constant differ for the precursor material on the flake itself from those on the substrate. If a flake is very small (Fig. 1 a.ii)), material adsorbed on the flake is very likely to diffuse to the flake edge to contribute to the increase of the flake size by direct growth. The degree of material supply via the flake itself is expected to depend on the actual flake area A for small flakes. Once a flake becomes rather large ( Fig. 1 a.iv)), the material supply via the flake and thus the growth would also become proportional to the edge length L due to a limited diffusion range. We propose that the material supply via the flake area A itself, even for very large flake areas, depends (approximately) on A, resulting in the growth rate G 2D within the 2D term of Eq. (1) valid for a wide range of flake sizes. This hypothesis is derived from the fact, that the precursor material consists of the same atoms as the 2D TMDC flake. Therefore, we anticipate complex, dynamic mechanisms taking place on/within the flake itself. The dynamics are discussed in more detail below, once the mechanisms underlying the etching rate are described.
For the etching rate E, we again introduce individual rates, E 1D (proportional to the flake edge length) and E 2D (proportional to the flake area). A flake size independent etching rate E 0D is not considered because no etching takes place without a flake being present. As soon as a few atoms agglomerate, already the smallest resulting agglomerate (nanoflake) has a spatial extension, so its decrease in size can be described by the size dependent etching rates E 1D and E 2D .
The etching rate E 1D mirrors that 2D TMDC flakes have an increased chemical reactivity at their open edges with respect to their pristine basal planes. This manifests, for example, in an increased reactivity with oxygen [52] or in an increased catalytical activity at the edges [6,7,14]. Therefore, it seems intuitively reasonable, that etching also occurs preferentially at 2D TMDC edges as reported by Lv et al. for pristine 2D TMDC nanoflakes [45]. Obviously, the release of built-in atoms at edges (enhanced by oxygen) directly contributes to the decrease of the flake size.
Less intuitive are the mechanisms for the change of the flake size that may account for the etching rate E 2D as well as for the growth rate G 2D . In the following, we will discuss various dynamically interacting mechanisms on and within the 2D TMDC flake at elevated temperatures.
Firstly, we start with 2D TMDC flakes heated up to temperatures significantly lower than their typical CVD temperatures (> 650°C) and at or close to ambient pressure. In this case, it is experimentally observed that the flakes begin to degrade and finally decompose completely [10,[42][43][44][45][46][47][48][49] [and see SI 1 ( Fig. S1)]. However, during this kind of degradation, the flakes do not become smaller from the edges. Instead, the atoms are released also from the basal planes of the 2D TMDC flakes. In some of these studies, annealing was intentionally performed with oxygen being present in the atmosphere. These include a comprehensive study by Cullen et al. showing for ten of the most common TMDCs degradation under ambient conditions [49]. For all of these TMDCs the degradation temperature is spectroscopically determined to be (far) below 400°C. That oxygen plays important role in the etching process is experimentally evident from the study by Yamamoto et al.: While no etching takes place in 2D MoS 2 under Ar/H 2 atmosphere at 350°C, etching is observed under Ar/O 2 atmosphere already at temperatures around 300°C [44]. However, we observe such an etching effect in 2D WS 2 under completely inert Ar atmosphere at temperatures above 300°C [48] [and SI 1 (Fig. S1)]. The reason might be small leakages that still let small amounts of air (oxygen) into the annealing system. On the other hand, because hydrogen binds oxygen, and because basal plane etching occurs in 2D MoS 2 even under Ar/H 2 atmospheres at temperatures in the range of 400-500°C [10], oxygen apparently only has a promotive role but is not necessary. If oxygen is present, the formation energy of S vacancies in pristine 2D MoS 2 basal planes becomes indeed negative [50]. However, the calculated oxygen dissociative adsorption barrier on pristine MoS 2 is rather large (1.59 eV) [53]. At sites of S vacancies in the basal plane of 2D MoS 2 , the oxygen dissociative adsorption barrier is halved [53]. Hence, it is more likely to extend pre-existing defects (a certain number of intrinsic defects is always present) in the basal plane of 2D TMDCs, than to create new ones. This is supported by experiments giving evidence for grain boundaries and induced vacancies to be the preferred sites for the release of built-in atoms [45][46][47][48] [and SI 1 (Fig. S1)] and by studies showing the creation of defect clusters in form of triangular pits in 2D TMDCs [10,[42][43][44][45][46]48]. The latter mechanism is often referred to as anisotropic (oxidative) etching. The increased chemical reactivity of defect sites is consistent with experimental studies reporting a high catalytic activity of defect sites in 2D TMDCs [9,11,12,15], rendering these sites to be chemically more like 2D TMDC edges than the pristine basal planes. Density-functional theory (DFT) calculations further confirm the increased chemical reactivity (catalytic activity as well as oxidation) at defect sites in the basal plane of 2D TMDC flakes [8,9,12,13,53].
From the previous paragraph we conclude, that etching does not only apply to the edges of 2D TMDC flakes (L-dependent/1D component), it also has an A-dependent/2D component (E 2D ). However, the etching process at the basal plane does not result in a reduction of the 2D TMDC flake size as reported by Chen et al. [32] (relevant data points in Fig. 3 a)) and as shown by our own data in Fig. 3 b). The major differences between this experimental observation of shrinking flakes and the studies mentioned in the previous paragraph are the conditions, under which the experiments are performed: (i) the latter are performed at much lower temperatures and (ii) in the absence of (at least one) precursor atom species.
If higher temperatures are applied, i.e. temperatures typically used in CVD and thus in the studies showing shrinking flakes (Fig. 3), diffusion is also facilitated. For example, it has been shown by transmission electron microscopy (TEM), that even already built-in atoms are able to diffuse within the 2D TMDC lattice, if a neighbouring atomic site is empty (vacancy) [50,51,[54][55][56]. This effect could be called also defect/vacancy diffusion. Once a vacancy reaches a 2D TMDC flake edge due to diffusion, the defect vanishes by reducing the flake size. This mechanism corresponds to an A-dependent etching component.
The vacancy diffusion barrier has been calculated by DFT for 2D MoS 2 and MoSe 2 to be between 0.6 and 2.9 eV [8,50,51,54,55,57,58]. The actual calculated barrier value depends on the type of vacancy (transition metal vacancy, single chalcogen vacancy, double chalcogen vacancy) and on the environment of the diffusing vacancy. For instance, the diffusion barrier of a single S vacancy is strongly reduced once an additional vacancy exists on a neighbouring site [8,54]. Therefore, pairs of single S vacancies would migrate faster through the 2D TMDC lattice. Because existing single S vacancies facilitate oxidation [53] as discussed above and thus improve the release of neighbouring S atoms, the basal plane etching may initially enhance the vacancy diffusion velocity. When S vacancies agglomerate, they tend to form vacancy rows as experimentally observed by TEM even at room temperature [51,54]. At high temperatures, the number of rows decreases, but their length increases [51]. DFT calculations confirm that these vacancy rows are energetically favored [8,54]. Due to a large diffusion barrier, S atoms at the edge of the vacancy row are unlikely to diffuse into the row [8]. Instead, atoms within the vacancy rows (both S and Mo atoms for 2D MoS 2 ) are able to migrate rapidly through the lattice [51]. Because in this way a lot of material can be moved, such "channels" are important for the formation of triangular pits within the basal planes of 2D TMDCs [51]. In this way, at 800°C, a triangular pit with a diameter of a few nm can be opened at the end of a vacancy row within one minute.
It is very likely, that the vacancy diffusion observed experimentally is triggered by the high kinetic energy of the electrons during the TEM measurements mostly performed at room temperature. However, Lin et al. report similar morphological structures formed by defect diffusion within a 2D MoS 2 lattice after annealing at 700°C in high vacuum as within the 2D MoSe 2 lattice after extensive defect diffusion triggered due to the electron beam [50]. Hence, such high temperatures, which are also typical for CVD of 2D TMDCs, may also be sufficient for a reasonable high diffusion of defects.
In addition, during CVD, simultaneous to defect creation and diffusion, the growth still takes place, i.e. new precursor atoms adsorb, diffuse, and are built in, if they reach an appropriate site. Not only edges are appropriate sites, but also the diffusing vacancies. Hence, adsorbed precursor atoms, which statistically would not be able to reach the edge (diffusion range), would at least compensate the reduction of the flake size due to the release of built-in atoms and subsequent defect diffusion to the edges (indirect growth). In SI 1 Fig. S1 we demonstrate the influence of precursor atoms existing in the gas phase on the degradation velocity. We compare the degradation of 2D WS 2 flakes under Ar, and sulphur containing Ar atmosphere and find a reduced degradation velocity, if S atoms are present.
Summarising up to here, a CVD process is not only about the growth of the flakes. Rather, CVD is a highly dynamic process including adsorption, etching, diffusion (of adsorbed atoms and vacancies), agglomeration, healing, and growth. In this respect, the most dynamic region during growth is the flake surface A itself.
Taking all contributions to growth and etching into account, we arrive at the following differential equation to describe the dynamic behaviour of a 2D TMDC flake during synthesis: This is the basic equation of our dynamic growth/etch model. In general, and in accordance with the previous consideration of the dynamic mechanisms during CVD, this equation is composed of three terms with different dependencies on the flake size: one 2D term and one 1D term reflecting either the dependence of the changing flake size on the flake area A or the edge length L, and one flake size-independent 0D term.
Actually, besides the growth rates G nD , also the etching rates E nD in Eq. (1) would be time-dependent. While for the latter the time dependence stems from their dependencies on temperature T (t) and pressure p(t), the growth rates additionally depend on the amount of precursor source material M (t). For our following discussion and application of Eq. (1), we assume a constant temperature and a constant pressure during the entire growth process with the duration t. Therefore, only the growth rates G nD depend on t, or more precisely, on M (t).
Further, we reduce our model to only one solid precursor source resulting in the schematic process configuration shown in Fig. 1 c). The mathematical description of the growth rates G nD for only one solid precursor source, which is introduced below in the solution (3), are a good approximation for many cases under following conditions: either, when most of the time one of the two precursor atom species is abundant relative to the other one, and/or, when the specific evaporation rates v of both precursor sources are approximately equal (see SI 2 for more details). The reduction to one solid precursor source is also experimentally confirmed by the typically (equilateral) triangular shape of 2D TMDC flakes grown by CVD [10, 12, 14, 20, 22-29, 31-34, 37-41, 45, 48] as this type of shape occurs, when one of the two precursor atom species is abundant relative to the other one [27,34]. Therefore, it is justified to take only one solid precursor source into account for the reaction rate and thus the growth rates G nD in CVD.
From the typically (equilateral) triangular shape of grown 2D TMDC flakes, we derive the relationship A = Furthermore, as the flake size is often expressed by the edge length L (or by the lateral size in one spatial dimension) in literature [10, 12, 22-28, 30, 32, 38, 39, 45, 48], the solutions and results are presented as a function of the edge length L below. Nevertheless, we decided to use dA dt in Eq. (1) because A is directly proportional to the mass of the flake m using the two-dimensional density of one TMDC layer ρ 2D and, thus, also to the number of built-in atoms. Therefore, dA dt is proportional to the mass change dm dt . We believe that this convention renders our basic differential equation (1) to be more intuitively understandable.
Next, we want to comment on the t-dependency of the growth rates G nD . As this t-dependency stems from the depletion of the solid precursor source during a running CVD process -either by being consumed or by forming a passivation layer on its surface (so-called poisoning) [32] -, we describe the growth rates by the differential equation with the solution (initial condition: Equation (2) expresses that the growth rates G nD (n = 0, 1, 2) are proportional to the temporal change of the precursor source mass − dM dt . In other words, the more material from the material source moves into the gas phase, the more material can adsorb on the substrate surface (including the surface of the already grown 2D TMDC flakes) and contribute to the growth of 2D TMDC flakes. On the other hand, at constant temperatures, − dM dt is also proportional to the mass of material available in the precursor source M . v is the specific vaporisation rate and c nD is a proportionality factor which takes into account microscopic as well as macroscopic factors. The latter include transport losses (see c macros. in Fig. 1 c)) and the fact, that one precursor source supplies hundreds of flakes simultaneously. The microscopic factors are different for the three supply areas because of varying conditions for e.g. adsorption, desorption, diffusion and reaction probabilities (see c micros. in Fig. 1 b)). Hence, c nD in general is specific for each of the terms in Eq. (1).
Solution (3) for G nD (t) is a strictly monotonically decreasing function. Hence, with adequate constants c nD and M 0 , the growth rates G nD initially dominate over the etching rates E nD in Eq. (1). The specific vaporisation rate v finally leads to a dominance of the etching rate E. Therefore, solutions of Eq. (1) can describe 2D TMDC flakes, which first become larger and later on shrink again with time, and can thus in principle explain the experimental observations of Chen et al. [32] and our own data in Fig. 3.
Next, we will present and discuss actual solutions of the total as well as parts of  (1) and of the partial differential equations (4a)-(4c). a) Numerical solution of the total differential equation and the 2D solution 5c. b) Derivation of the total solution and the contribution of the three terms to the total differential equation 1. c) Analytical 0D, 1D, and 2D solutions (5a)-(5c).
the differential equation (1). Unfortunately, the total differential equation (1) is not analytically solvable with the growth rates from Eq. (3). We therefore begin with its numerical solution (see methods). The typical shape, namely, first growth rates G and then the etching rates E become dominant, is shown in Fig. 2 a), red curve. Initially, the curve rapidly increases to a maximum of the flake size (here: flake edge length L). Thereafter, it drops somewhat less rapidly, but still rather quickly. Figure 2 b) shows the derivative of the numerical solution (red) including the contributions of the three terms of Eq. (1). Obviously, after a short time, the 2D term (blue) becomes predominant in growth and etching. The inset of Fig. 2 b) illustrates that the 0D term (orange) and then the 1D term (green) dominate in the early stages of the flake growth.
In order to identify the leading term(s) and thus the dominant physical mechanism(s) we split the total equation into the following three partial equations, each for one of term: The analytical solutions -again in terms of the solution for the respective growth rate [solution (3)] -are as follows (with L(t = 0) = L 0 ): Note, that a closed form solution also exists for the 1D2D differential equation, i.e. equation 1 without the 0D term. This solution is in fact to complex to be of practical value. However, we show and discuss it in SI 3.
Exemplary quantitative curves of these solutions are shown in Fig. 2 c) and their parameters are listed in SI 4 (Tab. S2). The parameters were chosen so that for the 1D solution (5b) and the 2D solution (5c) the local maxima are congruent. Without an etching rate, the 0D solution (5a) has no local maximum of course, but it is a monotonically increasing function with the largest change existing for t → 0. The solution of the 1D and 2D partial equations differ in such a way that the curve of the 1D solution is more convex near the local maximum, while the curve of the 2D solution is comparatively sharp. Mathematically, this behaviour is related to the nature of the two solutions: the 2D solution corresponds in its form to the exponential function of the 1D solution.
It becomes evident that the solution of the 2D partial equation (5c) is very similar in shape to the numerical solution of the total differential equation (1). In order to elucidate this, the 2D solution normalized to its maximum has been added to the plot of the numerical solution in Fig. 2 a). The 2D solution mainly diverges from the total solution for small or large times (both corresponding to small flakes), which is consistent with the expectation from Fig. 2 b) that the total differential equation (1) and its solution is largely dominated by its 2D term. Hence, the 2D solution (5c) is widely applicable as an analytically derived approximation for the actual solution of the total differential equation. In this context, the parameter L 0 must be chosen in order to compensate the neglected nucleation and early growth stages.
In order to test our model for plausibility, we apply it to experimental data of Chen et al. for 2D MoS 2 grown on sapphire [32] and on our own data for 2D WS 2 grown on sapphire. Both data sets are acquired by analysing flake size distributions from several growth processes for a varying growth duration t at maximum (growth) temperature and are shown in Fig. 3 a) and Fig. 3 b), respectively. Because the 2D solution (5c) (blue) can be fitted to both data sets obviously better than the 1D solution (5b) (green), our model and in particular our hypothesis, that the 2D term is largely predominant for the 2D TMDC growth, is confirmed. Now, we will discuss the data point at 60 min in the data set of Chen et al. (red in brackets in Fig. 3 a)). Both, the 1D and 2D solutions are fitted to the data set except for this last data point. However, even if it is included, this data point significantly differs from the best fit of the 2D solution and, thus, increases the fitting error, see SI 5.
Considering the history of this data point from Chen et al. [32], it is possible not only to justify why the data point can be excluded in Fig. 3 a), but it even confirms the premises for our model. For all other data points, our model of initially growing and later shrinking flakes can be applied. At the growth duration of 60 min, in fact, Chen et al. found that the grown 2D MoS 2 flakes behave differently. At this time anisotropic oxidative etching is observed: Instead of shrinking, triangular pits are formed within the basal plane of the flake resulting in a fragmentation into many smaller triangular flakes. This fragmentation/anisotropic etching effect is beyond our model, which does not take into account degradation of the basal plane occurring once the growth rate is significantly smaller than the etching rate. Apparently, Chen et al. evaluated the size of the small flake fragments for their data set, so the data point is shifted towards lower values as expected from our fit of the 2D solution.
Moreover, the fragmentation of the flakes proves that etching occurs not only at the edges, but also and in particular on the basal plane of the flakes. Obviously, this etching mechanism is only observable at later stages of the growth process, once the supply of precursor atoms is significantly diminished. The reduced quantity of new atoms can no longer sufficiently compensate the etching of the basal plane, resulting in the formation of defect clusters turning into pits with a low mobility that prevent the diffusion to the edges. Subsequently, the pits continue to enlarge anisotropically in the preferred directions of the TMDC lattice, becoming triangular, which is typically found for the flake degradation in absence of precursor atoms [10, 42-46, 48, 51].

Conclusion
In summary, we have presented a mathematical concept to describe the growth as well as the etching processes of 2D TMDCs during their synthesis at high temperatures. The individual mechanisms are represented by a differential equation that is made up of three parts. By a detailed analysis we found: (i) the numerical solution of the total differential equation differs only very slightly from the 2D solution in wide range, (ii) the growth and shrinking phenomenon experimentally found can be best analytically approximated by the 2D solution. These results imply, that both the material supply for growth and the material loss/etching are largely determined by the size of the flake area. This is in contrast to the common assumption of the importance and dominance of the flake edges, also suggested (but not discussed) by Chen et al. Our findings are corroborated by the fact, that the 1D solution corresponds to growth dominated by edges and the fits to both experimental data sets, including the one of Chen et al., are rather poor.
The predominance of the 2D solution of our model alters and advances the understanding of how the synthesis of 2D TMDCs takes place. It is based on the interplay between highly dynamic mechanisms at the atomic level within the flake itself and is thus consistent with the expectations at the high temperatures typically used for synthesis.
Our model provides an explanation for the rapid growth of 2D TMDC flakes by CVD (typical average: order of 100 nm/s). In order to exceed the etching term for a sufficient period of time in solid precursor source based CVD to actually obtain 2D TMDCs as product, the initial growth rates have to be chosen unphysically large. This clearly reveals that metalorganic CVD has the key advantage of a continuous supply of material over time. Within our model, the growth rate of metalorganic CVD would no longer be time dependent because of non-depleting precursor sources. This allows the etching rate to be precisely compensated and a constant net growth rate to be set. Note however, that the method suffers from other disadvantages such as small grain sizes, for example Kang et al. [59].
Because the actual synthesis of TMDC monolayer is based on the in-and on-flake dynamics due to the high process temperature, we finally conclude with a hypothesis for the frequently observed mutlilayer growth. Because the dynamic processes decay with decreasing temperature, we believe, that multilayer growth primarily takes place in the cooling phase. That is, once the temperature is sufficiently reduced, a transition state with a time constant depending on the cooling rate appears. Within this transition, remaining precursor atoms would continue to adsorb on the flake surface. Because their diffusion constant as well as the release of atoms (i.e. the etching rate is reduced at low temperatures, the precursor atoms might merge to additional layers on the first one. Optimizing growth during the cooling phase could thus be a successful strategy to suppress or enhance bilayer formation.

Materials and Methods
Chemical Vapour Deposition Tungsten disulphide (WS 2 ) flakes for the time-dependent study shown in Fig. 3 b) has been grown by chemical vapour deposition on c-face sapphire substrates. Therefore a custom made process system consisting of a heating belt and a tube furnace (ThermConcept ROS 38/250/12) was used, that in this way provides two heating zones in a quartz tube. A mass flow controller is used to adjust an argon (Air Liquide, 99.999 %) flow through the tube. The sapphire substrates are cleaned by an ultra sonic bath in ethanol and prepared by homogeneously spreading of individual WO 3 powder grains (Alfa Aesar, 99.8 %) on their surfaces. The substrates are deposited in a ceramic crucible in the downstream heating zone (tube furnace). In the upstream heating zone (heating belt) a second crucible with 160 mg sulphur (Sigma-Aldrich, 99.98 %) is deposited. After sealing the quartz tube, it is flushed by argon gas. The upstream heating zone is heated up to 150°C and the downstream heating zone to 800°C. The maximum temperatures are hold for 15 to 45 min. During the whole process a constant argon flow of 10 Ncm 3 /min is used. The pressure in the tube was close to ambient pressure.
Numerical Solving Differential Equation 1 is solved numerically by Mathematica (Wolfram Research, Inc.).