Determining by Raman spectroscopy the average thickness and N-layer-specific surface coverages of MoS₂ thin films with domains much smaller than the laser spot size

Introduction

The advent of two-dimensional (2D) layered materials beyond graphene has initiated a new field of research [1-3]. In the family of 2D layered structures, transition metal dichalcogenides (TMDs) have attracted considerable attention from academia and regarding potential applications [4-9] because of a number of remarkable properties [10-12]. Particularly, it was found that the properties of layered TMDs drastically change when their thickness is reduced to a monolayer [13,14]. Layered TMD structures have a graphite-like structure with each graphene sheet replaced with an X–M–X or MX₂ triatomic layer, where X is a chalcogen atom (e.g., sulfur, selenium, or tellurium) and M is a transition metal atom (e.g., molybdenum or tungsten) [10].

Among the layered TMD materials, molybdenum disulfide, MoS₂, is of particular interest in optoelectronic applications because of its transition to a direct bandgap semiconductor with very high photoluminescence quantum yield when thinned down to a monolayer [13-17]. Its unique electronic and optical properties could provide an edge in many future applications.

The multilayers MoS₂ structures are of the most common 2Hc type, where atomic layers are arranged in such way that the stacking between two adjacent layers corresponds to a twist angle of θ = 60°, and any Mo atom is sitting on top of two S atoms of the adjacent layers [18,19]. However, during the synthesis process (e.g., chemical vapor deposition (CVD) synthesis) or when using precise transfer or AFM tip manipulation techniques [20], twisted MoS₂ can be formed with two adjacent layers stacked with a relative twist angle (θ) varying from 0 to 60°. Such twisted-layered MoS₂ structures can exhibit a variety of interesting physical properties including unconventional super conductivity [21,22], non-linear optics [23,24], and moiré excitons [25].

Because the properties of MoS₂ flakes are first a function of their thickness, or layer number (N), it is of a primary importance to determine the N of MoS₂ flakes, including twisted MoS₂ flakes and defective MoS₂ flakes, synthesized by different ways. Independently of the structural organization between adjacent layers, a MoS₂ flake is usually named NL-MoS₂, or simply NL, with N being the number of MoS₂ triatomic layers, which defines the thickness of the flake.

Several optical techniques have been developed to identify the N of MoS₂ flakes produced by different methods. Among these techniques, Raman spectroscopy is widely used thanks to its convenience, non-destructiveness, and sensitivity to materials change, including strain, temperature, doping, and defects [26]. Concerning the characterization of MoS₂ flakes, different information can be derived from the measurement of the Raman features (frequencies, linewidths, and intensities) of intralayer phonon modes as well as those of the interlayer modes, the so-called layer breathing (LB) modes and shear (S) modes.

Recently, we have developed the reproducible direct growth of wafer-scale MoS₂ thin films on SiO₂/Si substrates by direct liquid injection pulsed-pressure chemical vapor deposition (DLI-PP-CVD) using low-toxicity precursors [27]. Such MoS₂ thin films showed good stoichiometry (Mo/S = 1.94–1.95) and the potential for high photoluminescence quantum yield. However, atomic force microscopy revealed that they are constituted of nanoflakes (with a lateral size of typically 50 nm) with possibly a distribution of thicknesses. Furthermore, depending on the synthesis conditions, the MoS₂ surface coverage can be incomplete, and the thin film average thickness can vary. These samples thus have characteristics, especially thickness inhomogeneities smaller than the laser spot size, that differ from the ones used to establish Raman spectroscopy-based MoS₂ layer counting methods [26,28-33]. In this context, the primary purpose of this work is to develop and validate an approach for determining the average thickness of such sub-laser spot size inhomogeneous MoS₂ thin films using Raman spectroscopy.

First, we reassess here as a ground work the information that can be derived from the Raman spectra of MoS₂ flakes for the evaluation of their thickness, N. Different Raman criteria for the determination of the thicknesses of MoS₂ flakes are first recollected; after the specification of the experimental protocol, domains and limits of application of these criteria are precisely defined from measurements performed on reference samples. These samples are well-characterized, either mechanically exfoliated or standard CVD MoS₂ large flakes deposited on 90 ± 6 nm SiO₂ on Si substrates. Then, we determine which Raman information is relevant to estimate the average thickness of MoS₂ samples produced by the DLI-PP-CVD method, which are constituted of nanoflakes and, thus, significantly different from the reference samples. Finally, an original procedure based on the layer breathing mode intensities is proposed to evaluate the surface coverage for each N, that is, the ratio between the surface covered by exactly N layers and the total surface, in DLI-PP-CVD samples.

Results and Discussion

Experimental procedure

To define a robust experimental Raman protocol to evaluate the thickness of a MoS₂ flake (i.e., its number of layers, N), it is first necessary to specify some parameters that can have a direct influence on the quality and accuracy of the results. The first parameter is the wavelength of the incident laser light used in the Raman experiments. As it will be detailed in the following, the measurements of the frequencies, linewidths, and intensity of first-order Raman active phonon modes of MoS₂ have to be obtained with good accuracy in order to evaluate the thickness of a MoS₂ flake. These phonons modes are (i) the in-plane phonon mode involving relative motion of Mo and S atoms with E′ symmetry for a monolayer (E¹_2g for bulk) and (ii) the out-of-plane phonon mode involving only out-of-plane motions of S atoms with A′₁ symmetry for a monolayer (A_1g for bulk). These modes are located around 385 and 405 cm⁻¹, respectively, in neutral and defect-free MoS₂ monolayers [33,34]. More precisely, in MoS₂ multilayers, the symmetries of these phonon modes are E′ and A′₁ for an odd number of layers, and E_g and A_1g for an even number of layers. For simplicity, hereafter when we will discuss the dependence on N of the features of these phonon modes, they will be simply referred to as in the bulk, E¹_2g and A_1g, independently of the number of layers.

A drastic change of the Raman spectra, especially in the frequency range of the A_1g and E¹_2g modes, occurs when the spectra are excited at an energy close to those of the A and B excitons located around 655 nm (1.89 eV) and 601 nm (2.06 eV), respectively, in MoS₂ monolayers [35,36]. When the incident laser energy is in the range of the A and B exciton energies (the so-called resonance conditions), other bands associated to different second-order processes are observed in the Raman spectra with a strong intensity, their frequencies, widths, and intensity depending on the excitation energy [36]. In addition, resonance conditions alter the symmetry selection rules of phonons of MoS₂ [35]. Some of the second-order bands overlap with the A_1g and E¹_2g modes, complicating the exact determination of the parameters of these modes recorded under resonance conditions. Furthermore, since the MoS₂ exciton characteristics (energy, width, and spectral weight) can be changed by several factors (e.g., stacking, strain, doping, and defects), the Raman intensities measured with a single laser wavelength close to exciton energies can be affected by external factors and differ for samples elaborated by different methods. For these reasons and in the aim to use Raman spectroscopy to count the number of MoS₂ layers, one must necessarily work under off-resonance conditions, that is, by using incident laser energy far from both exciton resonance energies. In this work, we chose to perform Raman experiments using 532 nm (2.33 eV) laser excitation, because this is sufficiently far from the energy range of A and B excitons [35].

All Raman spectra reported in this paper were recorded on different samples deposited on SiO₂/Si(100) substrate. Hence, the second parameter essential to define is the SiO₂ thickness. Indeed, the multiple reflection interferences that occur in the air/MoS₂/SiO₂/Si structure influence significantly the intensity of the phonon modes [28,37]. In this work, we chose to focus on substrates with a SiO₂ thickness around 90 nm, which corresponds to the first optimum value for MoS₂ monolayer (N = 1) Raman enhancement with a 532 nm excitation energy and also amplifies the signal in the wavelength range of photoluminescence emission (around 650 nm).

The third parameter to define is the power of the 532 nm light, P_λ, impinging the sample. Much of the Raman information available to evaluate the thickness of MoS₂ flakes is based on the following parameters: (i) on precise measurements of frequency of the A_1g and E¹_2g phonon modes of MoS₂. These lead to a precise knowledge of the frequency difference Δω_A−E. It was established that Δω_A−E depends monotonously on the number of layers, and Δω_A−E is largely used as criterion to evaluate the thickness of MoS₂ flakes [26,29,30]; (ii) on the precise evaluation of the integrated intensities of the phonon modes of MoS₂, namely A(A_1g) and A(E¹_2g), with respect to the integrated intensity of the 521 cm⁻¹ mode from a bare area of the oxidized silicon substrate, A₀(Si), used as an intensity reference [31], or from the silicon substrate underneath the MoS₂ flake, A_2D(Si) [28]; (iii) on the precise measurement of the A_2D(Si)/A₀(Si) intensity ratio [31]; and (iv) on the measurement of ultralow-frequency modes, the so-called breathing modes and shear modes. The frequencies and the number of LB and S modes allow one to identify the number of layers [32,33] and the presence of twist between adjacent layers from the vanishing of the S modes in twisted MoS₂ flakes [20,38-41].

Then, it is essential to determine the limit value of the laser power so that the above measurements are not affected by laser irradiation. Figure 1 shows the evaluation of the temperature of MoS₂ flakes prepared in different ways and that of the Si substrate as functions of the laser power impinging on the sample through a 100× objective (N.A. 0.9). The power was cycled between ≈5 μW and ≈2 mW. The temperature of MoS₂ flakes is evaluated from the Stokes/anti-Stokes intensity ratio of A_1g phonon modes (similar results are obtained using E¹_2g) and that of silicon from the Stokes/anti-Stokes intensity ratio of the 521 cm⁻¹ Si mode (see [42] for method details). While the silicon temperature is quasi-insensitive to P_λ, the temperature of MoS₂ flakes changes monotonically, reversibly, and quasi-linearly with P_λ (see inset of Figure 1). For MoS₂, we found an increase rate of 25–30 °C/mW for monolayers (1L-MoS₂) and 40–45 °C/mW for bilayers (2L-MoS₂). Usual effects of sample heating are the frequency shift of the phonon modes and their concomitant broadening. In Supporting Information File 1, the frequency and the linewidth of the Si mode as functions of the laser power are displayed (Figure S2). These two parameters are found to be insensitive to P_λ below 0.5 mW. More intriguing is the evolution of the frequency (Figure 2a) and width (Figure 2b) of the phonon modes of 1L-MoS₂ as a function of P_λ. For a simple thermal effect [43] and given the 25–30 °C/mW temperature increase rate determined previously, the frequencies of A′₁ and E′ modes should both downshift by 0.3–0.4 cm⁻¹/mW and the width of A′₁ should increase by ≈0.2 cm⁻¹/mW (the width of E′ should remain constant) contrary to what is observed in Figure 2a,b. To clarify this point, we present in Figure 2c (filled dots) the relative shift of the frequency of the A′₁ mode versus that of the E′ mode measured on CVD 1L-MoS₂ at different laser powers. In the same plot the expected shifts of these modes are reported (i) as functions of a pure thermal effect (Figure 2c, red line) [43] and (ii) as functions of the doping state (Figure 2c, magenta curve) [44]. Clearly the relative shift of the A′₁ mode frequency versus that of the E′ mode frequency as a function of P_λ significantly differs from the behavior expected by considering a simple thermal effect. Consequently, the results reported in Figure 2c clearly evidence photo-doping of 1L-MoS₂ concomitant with a thermal effect, as already observed for MoS₂ on SiO₂/Si [45] as well as for graphene [42]. Furthermore, the evolution of the A′₁ and E′ widths with P_λ (Figure 2b), that is, the weak change of the E′ width and the significant increase of the A′₁ width concomitant with the A′₁ frequency decrease, support this interpretation [44]. For a laser power smaller than 0.3 mW, photo-doping remains rather low, but it is the dominant contribution to the shift of the modes. Similar results were obtained on other samples including exfoliated 1L-MoS₂. The effects were found irreversible in some cases when P_λ exceeded 1 mW but remained always reversible if P_λ was kept below 1 mW.

Based on the above information, all Raman results reported and discussed in this paper were obtained by using a P_λ around 0.1 mW chosen as a good compromise between mitigating laser effects and maintaining measurement efficiency to ensure the accuracy of the Raman criteria discussed in the next part of this paper.

In summary, unless specified otherwise, all Raman spectra reported and discussed in this paper were recorded at an excitation wavelength of 532 nm, with a laser working-power close to 0.1 mW, and using a 100× objective (N.A. 0.9), on MoS₂ flakes or thin films deposited on SiO₂/Si substrate with a SiO₂ thickness of 90 ± 6 nm.

Application of Raman criteria to characterize MoS₂ flakes

In this part, we report and discuss the advantages and limits of some Raman criteria that were found to be efficient to derive the thickness (i.e., the number of layers N) of large MoS₂ flakes prepared by different ways, namely mechanical exfoliation and standard CVD (including twisted CVD 2L-MoS₂). Then, we discuss the application of Raman spectroscopy to characterize samples synthesized by DLI-PP-CVD. In contrast to the first two kinds of MoS₂ samples, the latter are constituted of nanoflakes with possibly a distribution of thicknesses and twist angles between adjacent layers of multilayer domains as well as a higher number of defects.

Exfoliated MoS₂ flakes as reference samples

We performed Raman experiments on mechanically exfoliated MoS₂ [1] that will serve as reference samples. The stacking sequence in exfoliated MoS₂ flakes is of the 2Hc-type [34]. The common feature of all these samples is to have a limited number of defects. Note also that all exfoliated flakes have a lateral size (few micrometers at minimum) significantly larger than the diameter of the laser spot. In such flakes, the exact number of layers, N, is determined by combining optical microscopy, spectral reflectivity, and the measurement of the breathing modes and shear modes in the ultralow frequency (ULF) range of the spectra [32-34].

One of the most popular criteria to determine the number of layers of MoS₂ flakes is the measurement of Δω_A−E, that is, the frequency difference between the A_1g and E¹_2g phonon modes [26,29,30]. Figure 3a shows the dependence of Δω_A−E on the number of layers measured on exfoliated MoS₂ flakes deposited on Si/SiO₂ substrates with four different SiO₂ thicknesses. As previously well documented in the literature, we confirm that Δω_A−E depends monotonously on the number of layers and does not depend on the SiO₂ thickness (Figure 3a). The separation between N and N + 1 values are larger than the experimental uncertainties (error bars in the graph) up to N = 3. The error bars start to overlap between N = 4 and N = 5. Comparison with data from the literature (see inset in Figure 3a) shows that this overlap occurs even between N = 3 and N = 4 when additional variability due to setup and samples is taken into account. Above N = 4, the separation becomes too small compared to the uncertainty. Thus, the measurement of Δω_A−E in exfoliated MoS₂ flakes allows one to evaluate with good accuracy the number of layers for N ≤ 3. It is then necessary to supplement the Δω_A−E criterion with others to reliably count thick multilayers. In addition, we will establish in the following that this criterion has to be taken with care to derive N in MoS₂ samples other than reference exfoliated MoS₂, because the A_1g and E¹_2g frequencies, and thus the value of Δω_A−E, can be affected by different factors such as stacking order, strain, doping, and defects which can be present in MoS₂ flakes prepared by other ways [44,46-49].

To evaluate the number of layers, we can also use information associated with the integrated intensity of MoS₂ phonon modes. Figure 3b and Figure 3d show, respectively, the dependences of the normalized integrated intensities of the A_1g and the E¹_2g mode as functions of N for four values of the SiO₂ thickness. For normalization, we use here an external reference, which is a bare Si(111) wafer with only native oxide. In the following, A(Si₁₁₁) stands for the integrated intensity of the Si(111) 521 cm⁻¹ mode. This reference is preferred to the Si(100) substrate with 90 ± 6 nm SiO₂ to avoid the effects of the SiO₂ thickness variations and crystal orientation. For comparison with other setups or references, the polarization ratio of our setup and the relative values measured on Si(100) with native oxide and 90 nm SiO₂ are given in Supporting Information File 1. As noted by several authors and predicted by the optical interference model, the normalized integrated intensities of MoS₂ modes, namely A(A_1g)/A(Si₁₁₁) and A(E¹_2g)/A(Si₁₁₁), increase first with N and then decrease showing a maximum for N = 4–5 for all SiO₂ thicknesses. Obviously, this non-monotonous dependence prevents using these measurands alone to evaluate the number of layers for N > 4. Moreover, a significant dependence of the MoS₂ Raman intensity on the SiO₂ thickness occurs for N > 2, pointing out the importance to determine precisely this latter parameter.

Another criterion to derive the thickness of MoS₂ flakes is the A_2D(Si)/A₀(Si) intensity ratio [31]. For the evaluation of this ratio, it is of great practical advantage to use the same silicon (the silicon below the oxide, which is Si(100) in the present work) in the measurement of A_2D(Si) and A₀(Si). A necessary precaution is that the Si(100) substrate orientation has to be kept the same for both measurements. Another advantage is to give a common origin to the plots of A_2D(Si)/A₀(Si) as a function of N (A_2D(Si)/A₀(Si) = 1 for N = 0) for any SiO₂ thickness.

Figure 3c displays the A_2D(Si)/A₀(Si) ratio measured on exfoliated MoS₂ flakes deposited on SiO₂/Si(100) substrates with four different SiO₂ thicknesses as a function of N. We confirm the monotonous decrease of the A_2D(Si)/A₀(Si) ratio with increasing N [31]. For each N, the value of this ratio depends on the SiO₂ thickness (Figure 3c; black, blue, green, and red symbols correspond to a SiO₂ thickness of 84, 87, 89, and 96 nm, respectively). Despite the monotonous dependence of the A_2D(Si)/A₀(Si) ratio, its dependence on SiO₂ thickness combined with experimental errors lead to the conclusion that the measured values for N and N + 1 can overlap for any N if the SiO₂ thickness is not known with good accuracy. For a given SiO₂ thickness, the gap between the A_2D(Si)/A₀(Si) ratio for N and N + 1 is sufficient to ensure a rather good reliability only for N ≤ 5.

In summary, for exfoliated MoS₂, considering jointly the three Raman criteria (i) value of Δω_A−E, (ii) value of the normalized integrated intensities of the A_1g and E¹_2g modes, and (iii) value of the A_2D(Si)/A₀(Si) ratio, one can unambiguously derive the number of layers as long as N ≤ 4 and the SiO₂ thickness is precisely known. It has also been suggested in the literature to use the intensity ratio A(A_1g)/A_2D(Si) (or equivalently A(E¹_2g)/A_2D(Si)). As it will be discussed in the following, we see two major problems with this approach. The first is the dependence of the Si signal on the crystal orientation and the SiO₂ thickness. The second relates to the fact that using this ratio, rather than using each measurand independently and contrasting them, even if more practical, can hide some information.

Finally, we compare in Supporting Information File 1, Figure S3, the dependence on N of A(A_1g)/A(Si₁₁₁) and A_2D(Si)/A₀(Si) for three SiO₂ thicknesses and two microscope objectives with different numerical apertures, N.A. = 0.9 (blue symbols in Figure S3), and N.A. = 0.5 (red symbols in Figure S3). We observe that the normalized integrated intensity of the A_1g mode, A(A_1g)/A(Si₁₁₁) (Figure S3a–c), is independent of the value of N.A. Concerning the dependence on N of A_2D(Si)/A₀(Si) (Figure S3d–f), we found that this ratio is slightly smaller for N.A. = 0.5 than for N.A. = 0.9 and is in a better agreement with the model of Li and coworkers [31] (black solid line in Figure S3d–f). (Note that in this latter work the experimental data on which the model has been adjusted were recorded using a numerical aperture N.A. ≈ 0.45).

MoS₂ flakes prepared by CVD

In this part, we analyze the pertinence of the previous criteria to derive the thickness of large MoS₂ flakes synthesized by CVD. In a first part, we probe the effectiveness of these criteria to evaluate the thickness of large standard CVD MoS₂ flakes. Such flakes have a limited number of defects and, like in exfoliated MoS₂, the stacking sequence is of the 2Hc type. In the second part, we examine the relevance of these criteria to evaluate the thickness of twisted CVD MoS₂ flakes.

Standard CVD MoS₂ flakes: As derived from the features of the LB and S ultralow-frequency modes, these samples do not show any twist between adjacent layers (presence of S modes for all flakes with N ≥ 2). The flakes are thus characterized by 2Hc stacking (or close to 2Hc stacking) and a low number of defects, and are named standard CVD MoS₂ flakes. On the basis of the latter features, the structure of these flakes is close to that of exfoliated MoS₂ flakes. However, the high temperature used in the CVD synthesis and interaction with the substrate can lead to lattice distortion and the presence of vacancies and doping. In the following, we limit our study to a number of layers N ≤ 4.

Figure 4a compares the values of Δω_A−E measured on exfoliated (Figure 4a, black symbols) and standard CVD MoS₂ flakes (Figure 4a, red symbols) for N ≤ 4. As previously, the exact number of layers is obtained by combining optical microscopy, spectral reflectivity, and number and frequencies of LB and S modes. For both kinds of MoS₂ flakes, Δω_A−E increases monotonously with N, but for the same N, the values of Δω_A−E are systematically larger in standard CVD MoS₂ flakes. We attribute this discrepancy mainly to a difference of strain states between the two kinds of samples. Exfoliated samples are mostly found with low or slight compressive strain, while CVD samples are under tension. Other factors such as doping, defects, or stacking were shown to lead to large changes of Δω_A−E [44,46,49]. In summary, the value of Δω_A−E is clearly and significantly sample-dependent. Consequently, Δω_A−E cannot be considered as a definitive criterion to derive the number of layers in any of MoS₂ flakes prepared in different ways. In other words, one cannot define a single master curve, Δω_A−E vs N, which would be valid for all the MoS₂ flakes independently of their preparation method or environment.

The dependencies on N of the normalized integrated intensities of A_1g and E¹_2g modes and the A_2D(Si)/A₀(Si) ratio measured on standard CVD flakes are compared with the average values of exfoliated samples with the same substrate SiO₂ thickness (Figure 4b–d). In contrast to Δω_A−E, the N dependencies of these intensities are very close in exfoliated and standard CVD MoS₂ flakes. Only A_2D(Si)/A₀(Si) and the normalized integrated intensity of the E¹_2g phonon modes slightly differ for N = 4. However, this may be due to the fact that the statistics is rather poor on this measurement, because this flake is rather small compared to the others. With regards to these results, these measurands give important information to evaluate the number of layers of 2Hc-stacked MoS₂ flakes independently of the elaboration procedure as long as N ≤ 4.

Twisted CVD MoS₂ flakes: Other interesting samples are large CVD MoS₂ flakes that present a twist angle, θ, between adjacent layers. We exemplify here the complexity to characterize such samples from the previous Raman criteria with the case of twisted MoS₂ bilayers. The identification of the bilayer character of the investigated flakes was unambiguously obtained independently from spectral reflectivity and optical contrast.

Figure 5a shows the low-frequency range of spectra recorded on three types of MoS₂ bilayer (named 2L-MoS₂ in the following), namely a bilayer with θ ≈ 30° (this sample belongs to the so-called twisted-bilayer family for which 0 < θ < 60° and is named in the following as θ-2L-MoS₂), the so-called 2Hc-2L-MoS₂ and the so-called 3R-2L-MoS₂. In the latter structure, the stacking between two adjacent layers corresponds to a twist angle of θ = 0°, and it is such that the S atoms of the top monolayer are superimposed on the Mo atoms of the bottom monolayer, and the Mo atoms of the top monolayer are above the hexagon centers of the bottom monolayer [50]. The spectrum in the low-frequency range is dominated by the contributions of the LB and S modes, the frequencies of these modes depending on the twist angle [20,39]. The LB mode emerges in the Raman spectra of all 2L-MoS₂ samples. In line with previous results [20], the Raman shift of the peak position of the LB mode in 3R-2L-MoS₂ is smaller than that of the 2Hc-2L-MoS₂, and the LB mode Raman shift in 30°-2L-MoS₂ is even smaller. Also, in agreement with the literature [20,38-41], the S mode vanishes in 30°-2L-MoS₂.

In Figure 5b, the 130–240 cm⁻¹ range of the Raman spectra recorded on monolayer, 2Hc-2L-MoS₂, and two θ-2L-MoS₂ is displayed. As well documented in the literature, this frequency range is dominated by the contributions of second-order Raman processes [20,26]. The general profile of the spectra is similar in 1L-MoS₂, 2Hc-2L-MoS₂, and θ-2L-MoS₂ with the exception that in the latter flakes, new bands, named FLA and FTA, are superimposed to the second-order Raman spectra. The FLA and FTA modes in θ-2L-MoS₂ are attributed, respectively, to folded longitudinal acoustic phonons and folded transverse acoustic phonons of the monolayer due to the presence of a moiré superlattice [20]. As shown in the literature [20], the frequencies of these modes depend on the twist angle (see Figure 2e in [20]). Unfortunately, these dependencies show a mirror behavior with respect to θ = 30°. This means that from given FLA and FTA positions, two values are possible: θ ∈ [0,30]° or its mirror 60° − θ. As a consequence, θ will be given in the range of 0–30° in all plots in Figure 6 with the possibility that the values attributed to θ-2L-MoS₂ could be 60° − θ instead. For instance, the data from both 2Hc-2L-MoS₂ and 3R-2L-MoS₂ are reported at θ = 0° in these plots. From the positions of the FTA and FLA, we claim that the spectra of the two θ-2L-MoS₂ displayed in Figure 5 correspond to 20°-2L-MoS₂ (Figure 5b–d, solid green line) and 30°-2L-MoS₂ (Figure 5a–d, solid red line), respectively.

The dependence of A_1g and E¹_2g modes on the twist angle (derived from the positions of FTA and FLA modes) is reported in Figure 5c. The frequency of the E¹_2g mode in 3R-, 2Hc-, and θ-2L-MoS₂ is downshifted with respect to its frequency in 1L-MoS₂, and it does not show a clear dependence on the twist angle. In contrast, the profile of the A_1g mode significantly depends on the twist angle, and a new mode, named FA′₁, appears on the high-frequency side of the A_1g mode. The FA′₁ mode is identified as Raman scattering from moiré phonons associated with the A′₁ dispersion curve of 1L-MoS₂. It is folded onto the zone center and, consequently, becomes Raman active [20]. Obviously, its frequency depends on the twist angle and the θ-dependence of the FA′₁ frequency was recently established both theoretically and experimentally (see Figure 3e in [20]). On the basis of these previous results, we have been able to evaluate the value of θ for each 2L-MoS₂ investigated from the position of the FA′₁ mode. The values of the angles derived from the position of FTA/FLA and FA′₁ are in close agreement.

The objective of this work is to characterize the thickness of all MoS₂ flakes. The relevance of the criteria based on the frequency (Δω_A−E) and normalized integrated intensity (A(A_1g)/A(Si₁₁₁)) of the A_1g mode has to be reevaluated in twisted 2L-MoS₂ flakes. As shown in Figure 5d, the normalized intensity of the 521 cm⁻¹ Si mode from the substrate underneath MoS₂ flakes is close in all the 2L-MoS₂ and independent of the twist angle.

Figure 6 summarizes and details the dependence on the twist angle of the four Raman criteria defined above for 2L-MoS₂. In all plots of Figure 6, the values of angles were derived from the positions of FTA, FLA and FA′₁. The values of the different criteria measured for θ-2L-MoS₂ are compared with the average values of the same criteria measured on exfoliated 1L-, 2L-, and 3L-MoS₂ flakes. In θ-2L-MoS₂, the value of Δω_A−E strongly depends on the twist angle and significantly differs from the average value measured in 2Hc-2L-MoS₂ (Figure 6a). For θ = 30°, the value of Δω_A−E is close to the one found in CVD 1L-MoS₂ [51]. In consequence, using Δω_A−E alone could lead to a wrong evaluation of the thickness of twisted 2L-MoS₂.

The normalized integrated intensity A(A_1g)/A(Si₁₁₁) significantly decreases when the twist angle increases (Figure 6b, red dots), and in 30°-2L-MoS₂, the value of A(A_1g)/A(Si₁₁₁) is close to the average value found in 1L-MoS₂ (Figure 6b, blue solid line). The behavior of A(A_1g)/A(Si₁₁₁) is opposite to the one of the normalized integrated intensity of the FA′₁ mode, A(FA′₁)/A(Si₁₁₁), the latter increasing with the twist angle (Figure 6b, gray squares). These results are in qualitative agreement with those reported in [40]. It can be emphasized that the integrated intensity of A_1g and FA′₁ bands taken together (Figure 6b, orange triangles) is close to the average value found for 2Hc-2L-MoS₂ (Figure 6b, red solid line). The reason for this compensation between A(A_1g) and A(FA′₁) is not clear yet, but it could present a practical advantage in the use of the global integrated intensity of the spectral band, located around the position of the A_1g mode for the evaluation of the thickness of twisted MoS₂ flakes.

We also observed a tendency for A(E¹_2g)/A(Si₁₁₁) to be slightly lower for θ-2L-MoS₂ than for 2Hc-2L-MoS₂ (or 3R-2L-MoS₂, which is similar), but to a lesser extent compared to A(A_1g)/A(Si₁₁₁), that is ca. 20% vs ca. 50% at maximum, respectively (Figure 6d). These results are also in qualitative agreement with those reported in [40]. Finally, only the value of the A_2D(Si)/A₀(Si) ratio seems to provide a robust/reliable information to characterize the thickness of MoS₂ flakes, since it is found largely independent of θ in all measured 2L-MoS₂ samples (Figure 6c). Even if further work is needed to complete the data presented here with other values of θ and twisted MoS₂ samples with N > 2, we anticipate that the value of A_2D(Si)/A₀(Si) ratio would be close in twisted and 2Hc-stacked MoS₂ multilayers. However, as previously recalled, the sensitivity of this ratio to the SiO₂ thickness and the gap between the A_2D(Si)/A₀(Si) ratios for N and N + 1 permit to ensure the determination of N with a rather good reliability only for N ≤ 5.

DLI-PP-CVD MoS₂ nanoflakes

The aim of this part is to define which Raman information is relevant to estimate the thickness of MoS₂ samples produced by DLI-PP-CVD. These samples are significantly different from the previous ones (exfoliated and standard CVD). Indeed, they are constituted of nanoflakes (with a lateral size of typically 50 nm, i.e., well below the laser spot size) with possibly a distribution of thicknesses and twist angles between adjacent layers of multilayer domains and a higher number of defects (the average inter-defect distance ranges from 3 to 6 nm as estimated from the LA and A_1g intensity ratio [52]). In addition, the MoS₂ surface coverage is a priori unknown and can be incomplete. It is then necessary to implement a first check criterion that ensures that the thickness estimation method based on the comparison with results obtained on exfoliated samples is still valid. More generally, this point is critical for the characterization of samples synthesized using new methods or new precursors that can lead to the co-deposition of several by-products (such as carbon, oxides, and metals), which can significantly change the measured Raman intensities. Based on the results presented in the previous sections, we have shown that the value of the A_2D(Si)/A₀(Si) ratio provides a robust/reliable Raman information to characterize the thickness of MoS₂ flakes for N ≤ 5. However, this parameter does not rely unambiguously on the presence of MoS₂. The deposition of any other material would influence its value and could lead to a wrong estimation. In the most general case, the sample characteristics are not perfectly known and can be significantly different from the reference characteristics. As a consequence, it seems mandatory to compare the thickness estimated from the A_2D(Si)/A₀(Si) ratio with other measurands directly related to the presence of MoS₂. To this aim, we propose to use jointly the normalized integrated intensity of the MoS₂ phonon modes, namely A(A_1g)/A(Si₁₁₁) and/or A(E¹_2g)/A(Si₁₁₁).

In Figure 7, the values of A(A_1g)/A(Si₁₁₁) (Figure 7a) and A(E¹_2g)/A(Si₁₁₁) (Figure 7b) are plotted as functions of the value of A_2D(Si)/A₀(Si). In these graphs, the data obtained on DLI-PP-CVD samples are compared with the average reference measurements established previously on exfoliated MoS₂ deposited on Si/SiO₂ substrate with the same SiO₂ thicknesses, namely 96 nm (red open dots in Figure 7) and 87 nm (blue open dots in Figure 7). Note that in the exfoliated samples, the exact number of layers N is perfectly known and given on the plots of Figure 7 close to corresponding open dots. The idea behind this representation comes from the expectation that the presence of contaminations or deposition of others species would have a different impact on the Raman intensity coming from MoS₂ in the film and on the one coming from the Si substrate underneath the deposited thin film. It is, thus, expected that the measurements on contaminated or highly defective MoS₂ thin films will fall off the reference curve. Indeed, data obtained on poorly crystalline MoS₂ films synthesized by DLI atomic layer deposition (not shown) are found systematically and significantly below the corresponding reference curve. Concerning the DLI-PP-CVD samples presented in Figure 7, the A(A_1g)/A(Si₁₁₁) vs A_2D(Si)/A₀(Si) dependence is found fully compatible with the respective reference exfoliated curves (Figure 7a). The A(E¹_2g)/A(Si₁₁₁) vs A_2D(Si)/A₀(Si) data points mainly agree for thin layers (A_2D(Si)/A₀(Si) > 0.8) but fall systematically below the corresponding reference exfoliated curves for thicker layers (0.8 > A_2D(Si)/A₀(Si) > 0.6) as shown in Figure 7b.

Another way to compare the results is estimating the thickness of DLI-PP-CVD samples by interpolation from exfoliated data of the measured values for A_2D(Si)/A₀(Si), A(A_1g)/A(Si₁₁₁), and A(E¹_2g)/A(Si₁₁₁). In Figure 8a (respectively 8b), the average number of layers () obtained using A(A_1g)/A(Si₁₁₁) (respectively A(E¹_2g)/A(Si₁₁₁)) are plotted as a function of the number derived from A_2D(Si)/A₀(Si). It can be emphasized that non-integer values are found for , indicating the presence of a mix with unknown proportions of bare substrate (0L), 1L-MoS₂, 2L-MoS₂, 3L-MoS₂, and so on in the investigated DLI-PP-CVD films. It is also noticeable that the errors of estimated from A(A_1g)/A(Si₁₁₁) become larger when is close to 3 as a consequence of the smoother dependence of this parameter with . In agreement with the conclusion drawn above from Figure 7, Figure 8a illustrates the coherence between the values of derived from A(A_1g)/A(Si₁₁₁) and A_2D(Si)/A₀(Si). All data remain close to the red solid line that represents the ideal relation y[ via A(A_1g)/A(Si₁₁₁)] = x[ via A_2D(Si)/A₀(Si)]. Figure 8b as well confirms that the values of derived from A(E¹_2g)/A(Si₁₁₁) and A_2D(Si)/A₀(Si) agree well for < 1.5, but the values of from A(E¹_2g)/A(Si₁₁₁) are systematically lower than those obtained from A_2D(Si)/A₀(Si) when > 1.5. One explanation could the presence of a larger proportion of multilayer regions in the thicker samples, for which, if they are twisted, A(E¹_2g)/A(Si₁₁₁) has been shown to be attenuated in the previous section. If so, the question then arises why the same behavior is not observed for A(A_1g)/A(Si₁₁₁) contrary to what would be expected. A possibility could be that because of the observed broadening of the A_1g mode in DLI-PP-CVD samples (presumably due to local heterogeneities in terms of doping, strain, defects, or thickness), the FA′₁ mode becomes indistinguishable from the A_1g mode. As a consequence, the intensity of FA′₁ would merge with A(A_1g) and compensate its attenuation. Other explanations relying on the presence of defects or strain cannot be disregarded, and further works are needed to fully clarify this point.

In order to further confirm the validity of the estimations of for DLI-PP-CVD samples, we compare in Figure 8c the values derived from A_2D(Si)/A₀(Si) with the ones obtained independently from spectral microreflectivity. A very good agreement is found between the two series of data, establishing definitively the relevance of the A_2D(Si)/A₀(Si) ratio to give with good accuracy the average thickness of DLI-PP-CVD MoS₂ samples for ≤ 3. This agreement justifies the use of the values of derived from A_2D(Si)/A₀(Si) as abscissa axis in the previous plots.

Finally, in Figure 8d the frequency difference between the A_1g and E¹_2g phonons is plotted as a function of estimated from A_2D(Si)/A₀(Si) for DLI-PP-CVD samples and compared to the data obtained on exfoliated and CVD MoS₂. DLI-PP-CVD data are distributed between the two curves obtained from the reference samples. This further confirms that this measurand cannot be used to evaluate with good accuracy their average thicknesses. Nevertheless this comparison can be informative, showing that samples with < 1 are most certainly mainly composed of 1L-MoS₂ and suggesting that the proportions of 2L-MoS₂, 3L-MoS₂, or more gradually increase with , which is compatible with AFM observations (not shown).

To get further insight on the number of layer distributions in DLI-PP-CVD samples, we have measured their ULF modes. Representative ULF spectra are shown in Figure 9a for samples with average thicknesses ranging from 0.6 up to 2.8 MoS₂ layers as estimated from A_2D(Si)/A₀(Si). Up to = 1.3, only the LB mode of 2L-MoS₂ is observed around 40 cm⁻¹ [20,38-41], showing that these samples can only be composed of 1L-MoS₂ and twisted 2L-MoS₂ plus possibly uncovered (bare substrate) regions. For thicker samples, the S mode of 2L-MoS₂ around 24 cm⁻¹ is additionally visible, as well as a signal between 25 and 30 cm⁻¹, corresponding to the LB and S modes of 3L-MoS₂ [32,33]. For ≥ 2.4, the LB mode of 4L-MoS₂ is also present around 21 cm⁻¹; there may also be a weak signal around 17 cm⁻¹ (corresponding to the LB mode of 5L-MoS₂) reflecting the presence of 5L-MoS₂. The S mode of 4L-MoS₂ could be present as well around 28 cm⁻¹, but it is hardly distinguishable from the LB and S modes of 3L-MoS₂. Thus, ULF Raman spectra give valuable qualitative information on the different N present in each sample. Quantitative information relies on the determination of the surface coverages for each N (σ_N), that is, the ratio between the surface covered by exactly N layers and the total surface. With N = 0 standing for the bare substrate and N_max being the largest number of layers present in the sample, the definition of the average number of layers can be written as

and the total coverage (including bare substrate areas) is obviously 100%:

AFM imaging (see Supporting Information File 1, Figure S4) reveals that for > 1.25, the surface is fully covered by MoS₂, that is, σ₀ = 0, which removes an unknown. In addition, for < 1.3, there is no signature of more than two layers, and we can set σ_N≥3 = 0 with confidence. Hence, for 1.25 < < 1.3, the set of Equations 1 and 2 simplifies to

This allows one to readily determine the two remaining unknowns σ₁ and σ₂, since is known from A_2D(Si)/A₀(Si).

Hereafter, a linear relationship between the Raman signal , which is the LB mode peak intensity of 2L (N = 2) areas (the broad but well-identified 40 cm⁻¹ peak), and the surface coverage is assumed, namely σ₂ = . The ratio α₂ = is determined from five samples (1.25 < < 1.3) for which we now have both the coverage σ₂ and the Raman signal .

Because α₂ is now known and assuming that the linearity between σ₂ and the Raman signal holds (which should be a good approximation for the thin multilayers considered here), σ₂ = can be obtained directly for all samples from the Raman spectra, and is thus no longer an unknown.

We now turn to the samples with < 1.25, which may present some bare substrate areas, so σ₀ and σ₁ are a priori unknown. Both and σ₂ are determined as explained above, and σ_N≥3 = 0 is again a safe estimate. Hence, Equations 1 and 2 reduce to the system

which may be solved trivially for σ₀ and σ₁.

A similar approach can be used for samples with 1.3 < < 2 that are fully covered (σ₀ = 0) and might present some trilayers but show no trace of thicker layers. We set σ_N≥4 = 0, and the system of equations reduces to

Both and σ₂ are determined as explained above. Thus, the system can be trivially solved for the two remaining unknowns σ₁ and σ₃.

At this point it would be natural to get the proportionality between σ₃ and a Raman signal attributed to 3L areas (N = 3), and proceed recursively to obtain σ₄ in slightly thicker layers, and so on and so forth. In practice this becomes challenging because of the uncertainty on the 3L (N = 3) Raman signal, which is less clear than the 2L (N = 2) peak. Another approach gave better results.

Three samples with between 2.75 and 2.85 are thick enough to neglect σ₀ and σ₁, yet thin enough for σ₅ to also be negligible as a first approximation. Equations 1 and 2 then reduce to

where and σ₂ are known, so σ₃ and σ₄ can be determined readily.

The LB mode of 4L-MoS₂, located around 21 cm⁻¹, is sufficiently separated from other modes to be identified (which was not the case for ), so that can be extracted from the spectra. From the three 2.75 < < 2.85 samples α₄ = is determined. Then, assuming again a linear relationship σ₄ = , the coverage by 4L (N = 4) layers can be determined for all samples. This removes another unknown.

Now the last remaining case of 2 < < 2.75 samples can be solved, as

give σ₁ and σ₃ directly, since , σ₂ and σ₄ are known.

The results obtained using this procedure are shown in Figure 9b where σ_N (with N from 0 to 4) is plotted as a function of , the average sample thickness. On this graph, all values of σ₂ (respectively σ₄) are calculated using (respectively ) even for the samples used to derive the proportionality coefficient α₂ (respectively α₄). As shown in Figure 9b for samples with 1.25 < < 1.3 (respectively 2.75 < <2.85), we find by this way −0.03 < σ₀ < 0.04 (respectively −0.01 < σ₁ < 0.01) with little fluctuations around the expected value of 0.

Just below the full coverage of the sample surface by MoS₂ (σ₀ > 0), both σ₁ and σ₂ increase with a slight tendency of σ₂ to increase faster than σ₁. Indeed, 1L-MoS₂ represents 80–90% of the deposited MoS₂ for = 0.5 and 70–80% for = 1.3. The maximum of σ₁ is reached around = 1.3 when the sample surface is totally covered by MoS₂ (σ₀ = 0), and σ₁ starts to decrease above this value. Around = 1.6, 1L-MoS₂ only represents 50% of the MoS₂. σ₂ continues to increase and reaches a maximum value of ≈50% around = 2 and then decreases for thicker samples. 3L-MoS₂ starts to appear after the substrate surface is completely covered by MoS₂ and increases continuously, representing about 50% of the thickest samples ( ≈ 2.8).

In order to verify our approach, we implemented a 2D growth toy model (see Supporting Information File 1 for details). The model results are shown in Figure 9b as full lines and give a good agreement with the experimental results. It should be noted that within this representation (σ_N = f()), the results of the model are remarkably robust to any parameter changes (the curves are almost insensitive to either doubling or halving the cell size and, thus, the advance rate, or to multiplying or dividing the growth rate by 5). In other words, this means that this comparison with the experiment cannot be used to validate any model parameters but demonstrates the relevance of the proposed procedure to estimate the σ_N from the experiments. Nevertheless, it has to be noted that while for < 1.3 the ULF Raman signature of 2L-MoS₂ remains very similar, it is not the case for thicker samples with the notable appearance of the S mode of 2L-MoS₂ around 24 cm⁻¹ [32,33]. This could mean that the stacking order distribution changes. As a consequence, the hypothesis based on the proportionality between and σ₂ would probably be less valid above = 1.3, and an error on the absolute values deduced can be expected. However, the appearance of the S mode of 2L-MoS₂ around 24 cm⁻¹ could also be related to NL-MoS₂ (with N ≥ 3) constituted of a stacking sequence where 2L are not twisted, for example, the so-called t(1+2)L, t(2+2)L, … structures [26,53]. In this case, our hypothesis would remain more appropriate. Despite this unknown as well as the other approximations made, we believe that the main tendencies can still be captured by the proposed analysis. Further works are needed to determine and improve the accuracy of the method.

Conclusion

In this work we have reviewed all Raman information leading to the evaluation of the thickness of MoS₂ flakes, that is, the layer number N. First, we have analyzed in detail the effects of some experimental parameters, namely the wavelength of the incident laser light used in the experiments, the power of the incident light, and the oxide thickness of the SiO₂/Si substrate on which the flakes are deposited, on the quality and accuracy of Raman results. Based on this analysis, an experimental protocol has been defined and systematically applied to large MoS₂ flakes (i.e., single-domain flakes much larger than the laser spot), including twisted MoS₂ flakes, prepared by different methods on the one hand and to MoS₂ thin films composed of nanoflakes prepared by the DLI-PP-CVD method on the other hand. Special attention was paid to the measurement statistics.

The limits of different Raman criteria which allow one to determine the thicknesses of MoS₂ flakes, namely (i) the value of Δω_A−E, (ii) the value of the normalized integrated intensity of A_1g and E¹_2g MoS₂ modes, and (iii) the value of the A_2D(Si)/A₀(Si) ratio, have been precisely studied in the different types of MoS₂ samples. We definitely confirm that Δω_A−E cannot be considered a robust criterion to derive the number of layers in MoS₂ samples. We found that the value of the A_2D(Si)/A₀(Si) ratio provides the most robust/reliable information to characterize the thickness of MoS₂ large flakes, especially since it is found largely independent of the twist angle. The limit of application of this criterion is N ≤ 5, under the condition that the SiO₂ thickness is precisely known.

We then apply this analysis procedure to DLI-PP-CVD samples constituted of nanoflakes with a lateral size of typically 50 nm (well below the laser spot size) with possibly a distribution of thicknesses and twist angles between adjacent layers of multilayer domains and a higher number of defects. Our results definitively establish the relevance of the A_2D(Si)/A₀(Si) ratio to give with good accuracy their average thickness , for ≤ 3. Nevertheless, we emphasize that this criterion is not only related to the presence of MoS₂ and can be influenced by several factors, such as the co-deposition of by-products or the presence of defects, leading to a wrong estimation of . We propose to combine A_2D(Si)/A₀(Si) with the normalized integrated intensity of the MoS₂ phonon modes, namely A(A_1g) and/or A(E¹_2g). Although limiting the application to ≤ 3, this approach enables the validation of the A_2D(Si)/A₀(Si) ratio to determine in the presented case, and we anticipate that it would avoid possible errors in unfavorable situations.

Finally, to get further insight on the number of layer distributions in DLI-PP-CVD samples, we have measured their ULF modes. An original procedure based of the measurement of the intensity of the layer breathing modes allows one to evaluate the surface coverage (σ_N) for each N. A 2D growth toy model gives a good agreement with the experimental results supporting the proposed procedure to estimate the σ_N from the ULF spectra.

Experimental

Samples preparation

Mechanical exfoliation

MoS₂ flakes were obtained by micromechanical cleavage of a MoS₂ crystal (HQ graphene) using scotch tape (Nitto) and PDMS slabs (Gel-pak). They were then transferred onto Si substrates with SiO₂ layers of different thicknesses, namely 84, 87, 90, and 96 nm. Flakes were selected by optical microscopy and their thicknesses were determined by optical contrast.

Standard CVD process

MoS₂ was grown by CVD on 87 nm SiO₂ on Si substrates using MoO₃ (Sigma-Aldrich, 25 mg) and sulfur (Sigma-Aldrich, 250 mg) powders as solid precursors using a 1 inch quartz tube furnace. MoO₃ powder was placed in the center of the heating zone of the furnace, while sulfur was placed upstream at the furnace inlet. Prior to growth, air was evacuated by flowing Ar (ultrahigh purity, Linde) for 15 min at 200 sccm, after which the tube was heated to 200 °C for 10 min. The temperature was then increased to 750 °C under Ar (100 sccm), and it was held at this value for 15 min before cooling naturally to room temperature.

Direct-liquid injection pulsed-pressure chemical vapor deposition (DLI-PP-CVD)

The 12 × 11 mm SiO₂/Si (with 87 nm or 96 nm SiO₂ thicknesses) substrates were cleaned in acetone (C₃H₆O, technical, Acros Organics), isopropanol (C₃H₇OH, 99.8%, Höfer Chemie GmbH) and deionized water (H₂O, Acros Organics) under ultrasonic agitation for 10 min each, before being blown dry with nitrogen. They were then immediately loaded on the susceptor of the reaction chamber (Annealsys MC-050) for deposition. Solutions of 0.001 M molybdenum hexacarbonyl (Mo(CO)₆, 98%, Strem Chemicals) and 0.002 M sulfur (S, 99.999%, Acros Organics) in anhydrous toluene (C₆H₅CH₃, 99.8%, Sigma-Aldrich) were used as precursors. The process is as follows: Following sample installation, the chamber is closed and brought to about 0.02 mbar. For monolayer depositions, it is imperative that the substrate is thoroughly cleaned and free of adsorbates. Therefore, to ensure complete desorption of remaining contaminants, the samples were kept for 30 min under vacuum at room temperature inside the deposition chamber. For the first part of the process, the pumping direction is reversed so that all species are pumped from the deposition chamber to the back of the reactor.

Nitrogen (800 sccm) is flowed through the chamber (200 sccm through the gas line, and 600 sccm through two injection heads) and the substrates are brought to 750 °C at a ramp of 2 °C/s. The reactor is kept in this state for 5 min for homogenization purposes. While still in reverse direction pumping, 0.3 g/min of both precursors are injected and vaporized to prepare the evaporation system for deposition. Then, the Mo(CO)₆ injection is stopped, the pumping direction is switched back to the deposition direction and hydrogen (40 sccm) is added to the gas mix. For 1 min, sulfur is injected to clean any remaining contaminants, and to prepare the surface of the substrate for MoS₂ deposition, then the deposition works in 20 s cycles. During one cycle, a single pulse of 3 to 10 ms of Mo(CO)₆ is injected while the S injection is set to 0.3 g/min. This 20 s cycle is repeated 80 to 160 times. The quantity of MoS₂ deposited is controlled by the quantity of Mo(CO)₆ injected, that is, the pulse duration and the number of cycles.

Raman spectroscopy

Raman spectra and maps were recorded using an Acton spectrometer fitted with a Pylon CCD detector and a 1800 grooves/mm grating (≈0.6 cm⁻¹ between each CCD pixel). The samples were excited with a 532 nm (2.33 eV) laser (Newport Millennia Prime or Cobolt Samba) throughan Olympus microscope objective either 100× (numerical aperture 0.9) or 50× (numerical aperture 0.5). The full width at half-maximum (FWHM) of the focused laser spot with the 100× objective is about 380 nm. Optimized focus conditions were checked for each measurement. The samples were mounted on a three-axis piezoelectric stage (Physik Instrumente) to ensure the precise positioning and focusing of the laser spot. A Si(111) wafer with only native oxide sample was used as a daily reference for the system. The laser power was continuously measured during acquisitions allowing for intensity normalizations of the Raman spectra at each point of the maps. All data presented in this paper, unless specified otherwise, are extracted from Raman maps constituted by hundreds to thousands points (see Supporting Information File 1 for an example), which were analyzed using a custom-made software. All reported points are the average values obtained by Gaussian fitting of the data distribution extracted from Raman maps (corresponding to hundreds to several thousands of spectra), and the error bars correspond to 99.7% confidence intervals (±3 standard deviations).

2D growth toy model

The model used the DynamicGrids.jl package, which was part of the Dispersal.jl framework [54], see Supporting Information File 1 for more details.