Raman Spectra of Titanium Carbide MXene from Machine-Learning Force Field Molecular Dynamics

MXenes represent one of the largest class of 2D materials with promising applications in many fields and their properties tunable by the surface group composition. Raman spectroscopy is expected to yield rich information about the surface composition, but the interpretation of measured spectra has proven challenging. The interpretation is usually done via comparison to simulated spectra, but there are large discrepancies between the experimental and earlier simulated spectra. In this work, we develop a computational approach to simulate Raman spectra of complex materials that combines machine-learning force-field molecular dynamics and reconstruction of Raman tensors via projection to pristine system modes. The approach can account for the effects of finite temperature, mixed surfaces, and disorder. We apply our approach to simulate Raman spectra of titanium carbide MXene and show that all these effects must be included in order to properly reproduce the experimental spectra, in particular the broad features. We discuss the origin of the peaks and how they evolve with surface composition, which can then be used to interpret experimental results.

Two-dimensional (2D) materials have attracted a lot of attention from the scientific community in the last two decades. Among them, MXenes represent one of the largest class of 2D materials with more than 40 synthesized compositions and many more still left to discover. 1 Due to their high electrical conductivity, MXenes show promise for applications in various fields, including energy storage, gas sensors, and electromagnetic interference shielding. [2][3][4] MXenes take the form M n+1 X n T x where M is an early transition metal (e.g. titanium), X can be carbon or nitrogen, T are the surface terminations which depend on the synthesis methods, and n = 1-4 defines the thickness of the layer. The first synthesized and the most studied MXene Ti 3 C 2 T x is obtained by selective etching of the Al layer from its precursor bulk phase Ti 3 AlC 2 . 5 Due to this process, the surfaces are passivated by functional groups from the etching solution, with the most common ones being -O, -OH and -F. Balls-and-sticks representation of pure -O and -OH surfaces of Ti 3 C 2 are represented in Fig. 1(a,b), respectively. Under the common wet-etching synthesis conditions, the surface is never purely O or OH but contains a randomly arranged mixture of these functional groups [6][7][8][9][10] Unfortunately, it has turned out to be challenging to not only control the ratio of these terminations, but also to quantify them.
One commonly used method to study the composition of 2D materials is Raman spectroscopy, which is a non-destructive method for characterizing the vibrational properties of molecules or solids. Raman spectra of Ti 3 C 2 T x has been measured and carefully studied experimentally. [11][12][13] To explain the experimental observations, Raman spectra have been calculated using density functional theory (DFT) for pure surfaces, 14,15 but the comparison is not straightforward. To illustrate these problems, in Fig. 1(c) we reproduce typical experimental Raman spectrum and the frequencies of Raman active modes of pure surfaces (-O and -OH) from DFT (description of the synthesis and Raman measurement for the experimental results are available in the Supporting Information). For the modes around 120 cm −1 , 200 cm −1 , and 700 cm −1 , the experimental frequencies appear to be close to, or in between, those of the calculations for pure surfaces, and thus these peaks can  2 , respectively. Titanium is in pink, carbon in blue, oxygen in red and hydrogen in white. (c) Experimental Raman spectrum of Ti 3 C 2 T x compared with the frequencies of the Raman active modes calculated using DFT. The red and blue lines correspond to calculations of pure Ti 3 C 2 O 2 and pure Ti 3 C 2 (OH) 2 , respectively. Full lines represent the in-plane modes while the dashed ones represent the out-of-plane modes.
be quite safely assigned to the calculated modes. More importantly, other peaks around 300-400 cm −1 and 600 cm −1 both show extremely broad features which cannot be explained by static calculations of pure surfaces. There are few possible origins for this discrepancy. First, static calculations give no information about the width of the Raman peaks. Calculation of anharmonic interatomic force constants 16 or molecular dynamics (MD) simulations are usually necessary to account for temperature effects and obtain realistic line widths. In the latter, Raman intensities can be obtained from the Fourier transform of the polarizability χ(t) autocorrelation function, 17 but requires tens of thousands of polarizability calculations. It has successfully been used to compute Raman spectra of solid ice, 18 liquid water, 19 and molecules, 20 but such calculations are very expensive and usually restricted to a small number of atoms. Secondly, wide peaks of the MXene Raman spectra might come from the structural disorder, i.e., the mixed surfaces. In order to simulate the different surface terminations, it is necessary to use large supercell, making the calculation of polarizabilities at every time step of an MD run computationally intractable.
With the recent developments in machinelearning force fields (MLFF) it has now become possible to carry out MD simulations for long times and large number of atoms efficiently, yet still nearly matching the accuracy of firstprinciples calculations. 21-26 Moreover, by using MLFF, it becomes possible to run MD over a large configurational space, such as different compositions and distributions of surface terminations in MXenes. Unfortunately, MLFF usually cannot be used to evaluate polarizabilities as they only model energies and forces. To compute the polarizabilities of these large supercells, Hashemi et al. recently developed a computational scheme (denoted RGDOS), which allows to obtain the Raman tensors of a supercell from a projection onto those of the unit cell. 27 The scheme has already been successfully applied to, e.g., transition metal dichalcogenides alloys (Mo x W 1 -x S 2 and ZrS x Se 1 -x ), 27,28 defects in MoS 2 , 29 and SnS multilayer films. 30 Adapting RGDOS for evaluating time-dependent polarizability χ(t) and combining it with MLFF MD trajectories would lead to a highly efficient scheme to obtain the Raman spectra at finite temperature and in large supercells.
In this work, we present a scheme for combining MLFF trajectories with RGDOS (or DFT calculated Raman tensors) to evaluate finitetemperature Raman spectra of complex systems. The method leads to highly efficient MD runs and calculations of the polarizabilities, allowing us to easily probe large configurational space. We apply the method to 2D titanium carbide MXene to study the effect of heterogeneous surfaces on the Raman spectra and attempt to give a better description of the experimental results. We first investigate the effects of temperature and mixed surface terminations on the Raman spectrum, and finally study the effect of additional disorder by considering the vibrational modes outside Γ-point. We discuss the origin of the peaks and broad features and how they evolve with the surface composition.
When using MD, the Raman spectrum is obtained as the Fourier transform of the autocorrelation function of polarizability χ: 17,31 where x(τ )x(t + τ ) τ denotes the autocorrelation function. The challenge therefore is how to obtain the polarizability of a large supercell for every configuration visited during MD at a reasonable computational cost. In RGDOS, 27 the Raman tensors of a large supercell are obtained by first projecting its vibrational modes onto those of the pristine unit cell and then combining the Raman tensors of the unit cell weighted with these projections. This means that only the Raman tensors of the unit cell have to be calculated, resulting in an immense reduction of computational effort. The scheme has already been successfully applied to study the effect of alloys, defects, and film thickness on the Raman spectra, all involving large supercells. [27][28][29][30] Since this method does not rely on the electronic structure of the supercell, it is perfectly suited to be used with empirical or MLFF MD simulations. Moreover, when combined with MD simulations, such method would also account for temperature effects and anharmonicity.
In this work, we adopt an approach conceptu-ally similar to the RGDOS. We start by writing the polarizability as a Taylor expansion around the relaxed positions r 0 , which leads to the following equation where χ 0 is the polarizability at r 0 and u represent the atomic displacements One can notice the resemblance between the partial derivative ∂χ ∂u and the definition of Raman tensors R m = ∂χ ∂e * m (in the Placzek approximation 32,33 ), where e * m are the mass-scaled eigenvectors of the unit cell. This suggests that we could use the set of eigenvectors {e * m } as a basis for the expansion, with Raman tensors giving the expansion coefficients and the "distance" given by the projection of u to e * m , written out as where P m (u) is the projection of u onto e * m . We first note, that since e * m do not form an orthogonal basis, the projection coefficients P m (u) have to be found by solving the following system of linear equations i e * m |e * i P i (u) = e * m |u .
Moreover, since only Γ-point modes contribute to first-order Raman scattering, the projections are here done from the supercell displacements u to the unit cell eigenmodes at Γ-point e * m . The final equation for the time-dependent polarizability is then When the eigenmodes and the Raman tensors of the unit cell and the MD trajectories are known, χ(t) can easily be obtained from equations 4 and 5. This method is therefore highly efficient and compatible with MLFF MD. We also note that in equation 5, the second order term could be added, essentially corresponding to the second-order Raman scattering. In this case, the sum should also be carried over different points of the first Brillouin zone (BZ) to account for the modes outside of Gamma-point. In the case of MXenes, there is no indication of dominant second order contribution and thus we have decided not to include the second order term.
All calculations are done using the Vienna ab initio simulation package (VASP). 34,35 We adopt the Perdew-Burke-Ernzerhof exchangecorrelation functional for solids (PBEsol) 36 and set the plane-wave cutoff to 550 eV. Phonon eigenmodes are obtained using the Phonopy software. 37 For determining the force constants, we used 4×4 supercells with the Brillouin zone sampled using a 4×4 k-point mesh. Once the eigenvectors are known, the Raman tensors were obtained from a centered finite difference scheme. Since MXenes are metallic, the polarizabilities have to be evaluated at finite excitation energy ω, i.e., using χ(ω) instead of ω → 0 limit normally used in nonresonant Raman. Although Placzek approximation is derived under non-resonant conditions, it has been found to work well also in resonant conditions. 38 Here, for the projections we used the eigenvectors and Raman tensors of -OH unit cell, the latter evaluated at ω = 516 nm. More information regarding resonant Raman and the choice of Raman tensors are available in the Supporting Information (see in particular Fig. S1). The frequencydependent polarizability (or dielectric function) is computed using a summation over empty bands. 39 To guarantee accurate evaluation of small changes in χ(ω), a denser k-point mesh of 48×48 (in the unit cell) was adopted and the number of orbitals was increased to 120 (roughly 4 times the recommended amount). For phonon modes and Raman tensors calculations, the electronic structure was relaxed with a precision of 10 −7 eV.
The machine learning force field (MLFF) was trained using on-the-fly machine learning, as recently implemented in VASP. 25,26 This method has already been successfully applied to different materials including different phases of zirconium 40 and hybrid perovskites MAPbX 3 . 41,42 Fig. 2(a) shows a schematic of the workflow and how MLFF, MD, DFT and RGDOS are all combined to obtain the Raman spectrum. In the training, we again used 4×4 supercells with a 4×4 k-points mesh, but the electronic structure convergence criterion was loosened to 10 −4 eV. The model was first trained for pure surfaces before be-  ing trained for mixed -O and -OH surfaces as well. F-terminated surfaces were not included in order to reduce the complexity of the MLFF model and since they are expected to behave fairly similar to OH-groups. 10,43,44 The vibrational frequencies are generally bit lower than for OH-surface 45 due to the slightly larger mass of F. Each step consisted of MD simulation in the NVT ensemble at 300 K for 20 hours (40000-70000 steps), and the MLFF configurations after each step carried on to the next step. In total, the training set contains 1364 configurations. The root mean square error (RMSE) between DFT and MLFF forces is 45.8 meV/Å. To benchmark the model, we compared the phonon dispersion curves from DFT and the MLFF model. Fig. 2(b) and (c) show the phonon bands and density of states (pDOS) for pure Ti 3 C 2 O 2 and Ti 3 C 2 (OH) 2 , respectively. Results are overall in good agreement for the dispersion curves as well as for the pDOS. Fortunately, the few modes where errors at Γ-point are somewhat larger, such as the 350 and 700 cm −1 modes in Ti 3 C 2 O 2 , happen to be Raman-inactive modes. We also note, that since MLFF is fitted into finite-temperature trajectories, the 0 K force constants used in the phonon dispersion curves may differ slightly from the DFT results (the effect of temperature is discussed below). Those benchmarks indicate that the model is sufficiently trained and that we can use it for reliable MD production runs.
Molecular dynamics were performed using the MLFF model in the NVT ensemble using a Nosé-Hoover thermostat 46,47 for 8x8 supercells. The time step was set to 1 fs and 10 5 steps were performed, resulting in runs of 100 ps. To reduce noise in the final spectra, five runs with different random distributions of the surface terminations are performed for each concentration. Using different distribution of the surface termination does not impact the main features of the resulting Raman spectra, see Fig. S2 in the Supporting Informations for more details.
We first investigate the effect of temperature on the width of the peaks. Fig. 3 Raman intensity Raman intensity DFT results presented in Fig. 1(c) for both systems, showing that our approach leads to reasonable spectra. Overall, the frequencies do not shift markedly with temperature, but there is now a realistic description of the peaks widths. For some peaks, such as the two modes at 100 and 200 cm −1 localized to titanium atoms, show weak dependence on temperature. While our results agree with experimental observations for the first, the width seems underestimated for the latter. For other peaks, the width increases with temperature. The A 1g carbon mode at around 700 cm −1 shows a fairly narrow peak at low temperature but it gets significantly wider with temperature, resulting in a better agreement with experiment at 300 K. Similar broadening with increasing temperature is observed for the in-plane modes of surface terminations and carbon at 300-400 and 500-600 cm −1 , respectively. Since surface terminations (-OH groups especially) should be strongly anharmonic, significant broadening could be expected for modes involving these atoms. In those cases however, the widths at 300 K are still underestimated and does not compare well to experiments. Overall, using finite temperature leads to a better comparison with experiments, but it is clear that something is still missing in order to correctly reproduce experimental observations. We next look at the effect of mixed surfaces. Fig.  4(a) shows the Raman spectra for different concentrations of surface terminations and Fig. 4(b) illustrates the distribution of surface terminations for different concentrations. The top and bottom lines show results for pure Ti 3 C 2 O 2 and Ti 3 C 2 (OH) 2 , respectively. The titanium A 1g mode at 200 cm −1 is largely unaffected by the composition changes. The other two sharp peaks from the titanium E g mode at around 120 cm −1 and the carbon A 1g mode at around 700 cm −1 both show a straightforward linear shift with O/OH ratio (a so-called onemode behavior). This suggests that these peaks could be used to find the surface termination composition from experimentally measured spectrum. To this end, we extracted the frequencies by isolating the peaks and fitted them using a Lorentzian function (see Fig. S3 in the SI for details), and collected the results in Table 1 together with experimental values. In particular, we include the results from our experimental spectrum shown in   Fig. 4(a) below the x=0.625 spectrum) and those from Ref. 11, where these peaks were found to shift depending on the sample preparation conditions. First, if we use the Ti A 1g mode as a probe for the accuracy of our calculations, we find that the calculated frequencies are 0-10 cm −1 lower than the experimental values. Second, by comparing the frequencies of the other two modes, best agreement with experiments is found for x around 0.5-0.625. Even accounting for the inaccuracies in the calculated frequencies, our results strongly point to all these samples having surfaces with mixture of surface terminations. Interestingly, the composition range found from this comparison also agrees well with that predicted by first-principles calculations in Ref. 9: O 0.75 OH 0.25 (x=0.75) at higher pH and O 0.5 OH 0.25 F 0.25 at low pH (corresponding to a total concentration of 0.5 for -OH and -F terminations). For these two peaks, calculations also reproduce well the widths of the experimental peaks.
In addition, there are the two in-plane modes at 300-400 cm −1 and 500-650 cm −1 , which show distinct two-mode behavior with changing composition. For the former, it seems that the two broad features in the experimental spectrum arise from these two modes. For the latter, it is difficult to distinguish the two peaks from the experimental spectrum, but as this feature is very broad it is consistent with being comprised of two peaks as found in our calculations. Also note that at concentrations close to x=0.5, the out-of-plane mode of carbon also seem to contribute at a similar frequency. Even though Raman spectra simulated using mixed surfaces agrees much better with experiments and would appear to contain all the main features, the agreement is clearly still far from perfect. In particular the width of some peaks, mainly those arising from the surface terminations and the carbon E g modes, are still greatly underestimated and by more than one might expect from instrumental broadening.
The peak broadening could arise from the interaction of surface terminations with water or other surface adsorbates, from the interaction with neighboring MXene layers, or from disorder arising, e.g., from point defects, grain boundaries, etc. We simulated the Raman spectra of multilayer MXene flakes, but this did not lead to marked broadening of these features as shown in Fig. S4(a) in the Supporting Information. Moreover, experimentally measured spectra of monolayer flakes, multilayer flakes, and colloidal dispersions all show largely similar broadening. [11][12][13] We also investigated the presence of (small amount of) water on the surface, but also did not yield the desired effect (see Fig S4(b) in the SI).
One might expect that the wet chemistry and strong acids used in the synthesis would lead to a large number of defects, in addition to the defects inherited from the precursor MAX phases. [48][49][50] Consequently there is also a large variety of possible defects: in addition to flake edges and grain boundaries, 2,48,51 the reported point defects include Ti vacancies and Ti adatoms, 52,53 C vacancies, 54 and substitutional O in C site. 55 Unfortunately, the dominant types of defects and their density is largely unknown, which makes simulating the role of defects on the Raman spectra challenging.
To circumvent this problem, we simulate the effect of disorder by considering the modes outside Γ-point. In ideal crystals, first order Raman spectra only comes from modes at Γ-point because of the condition k = 0. 56 However, in the presence of defects, the symmetry is lowered and the modes from the rest of the first Brillouin zone can contribute. 57,58 In fact, this effect is already present in our calculations for randomly distributed surface terminations, but clearly the degree of disorder from surface terminations is insufficient. To explicitly include the contributions from modes outside the Γ-point, we study here the phonon density of state using phonon-projected velocityautocorrelation function (VACF). 42,59-61 By projecting only onto the Raman active modes, we can isolate the corresponding part of the vibrational spectrum and in this way approximate Raman spectrum. Note however that the peak intensities do not correspond to those obtained from Raman and only the frequencies and lineshapes are relevant. Comparison to the total pDOS (i.e. containing both infrared and Raman active modes) is presented in the Supporting Information (see Fig.  S5).
We focus here on x=0.5 surfaces since we already determined this concentration as a good candidate to reproduce experimental results. Every Table 1: Frequencies (in cm −1 ) of the E g (Ti), A 1g (Ti), and A 1g (C) modes for different ratios of O/OH surface terminations, and comparison to experimental values.  point of the reciprocal space are not expected to contribute to the Raman spectra with the same weight. To find the weights, we study the pDOS for each wave vector along the Γ − M path, which are shown in Fig. 5(a). pDOS at and close to Γ mostly show low frequencies modes, coninciding with the Raman spectrum in Fig. 4. In experiment, these modes show narrow peaks which tends to indicate a small contribution from modes outside Γ. On the other hand, modes close to the edge of the first Brillouin zone show wide peaks around 300-400 and 600 cm −1 . Based on the experimental observations, we therefore choose larger weights for modes close to the BZ edges. The selected weights are represented in Fig. 5(b) and built from a sum of two Gaussians such that the weight approaches 0 at Γ-point and 1 at the BZ edge, except pDOS at Γ-point is also added with a weight of 1.
Resulting pDOS using these weights is represented in Fig. 5(c). When summing the weighted pDOS with pDOS at Γ, we find great agreement with experiment. The very wide peaks in the 300-400 and 600 cm −1 regions are well reproduced. The peak at 200 cm −1 also gets wider with the additional contributions from outside Γ, further improving the agreement with experiment. As previously discussed, pDOS does not reproduce well the intensities of peaks, which is especially true for the peak at 700 cm −1 . Also note that the spectrum now shows an additional peak around 500 cm −1 which also appears in the experimental results, even though the agreement in frequency is not perfect. In Fig. 4(a), this peak shows up weakly in the OH-surfaces, but disappear in the mixed surfaces. However, since this mode is relatively flat [cf. the phonon dispersion in 2(c)], the inclusion of modes outside Γ-point could lead to increased intensity. Thus, this peak seems to arise from the A 1g (OH) mode with its intensity enhanced by disorder, which we could confirm by inspecting the eigenmode-projected VACF in Fig. S5(b) in the Supporting Information.
Similar pDOS for other concentrations is available in the Supporting Information (Fig. S6). The results in Fig. S6 and those in Fig. 4 suggest that the relative intensity of the two broad peaks at 300 cm −1 [E g (OH)] and 400 cm −1 [E g (O)] and of the two peaks that form the broad peak at 600 cm  Figure S6 also illustrate that the surface concentration has a relatively minor impact on these peaks. This finding, together with previous computational predictions of a limited range of accessible surface compositions under typical synthesis conditions, 9 explains why the Raman spectra of Ti 3 C 2 T x appears to depend weakly on the adopted synthesis procedure and on the sample morphology.
In conclusion, we have developed a computational approach to efficiently compute the Raman spectra of complex materials at finite temperature, which is based on a combination of machinelearning force fields and a method to efficiently evaluate the time-dependent susceptibility. Our approach can be applied to many other materials with high anharmonicity and pronounced temperature effects. Alternatively, it can be used for systems requiring large supercells when modeling, e.g.,alloys and/or defects.
Here, we applied the approach to 2D titanium carbide MXene to study the effects of temperature and surface terminations on the Raman spectra. We find that the temperature plays a relatively small role and does not explain the wide peaks observed experimentally. On the other hand, mixed surfaces play an important role in the reproducing experimental Raman spectra. We identified three peaks that undergo straightforward evolution with the surface composition and thus could also be used to extract the composition from experimental spectra. Mixed surfaces lead to better agreement in the peak position for some peaks and, importantly, they also reproduce the two broad peaks observed at 300 and 400 cm −1 . Additionally, the modes outside Γ-point are found to play an important role, activated by the disorder such as defects. In particular, modes close to the Brillouin zone edge are found to contribute significantly to the widening of peaks around 300-400 and 600 cm −1 , as well as the appearance of an additional peak around 500 cm −1 . Only by including the effect of heterogeneous surfaces and contribution from modes outside Γ-point, a good agreement with the experimentally measured spectra can be achieved. and its application to liquid S-methyloxirane.

S1
Synthesis and Raman measurement of Ti 3 C 2 T x An aqueous dispersion of Ti 3 C 2 T x was synthesized using MILD method. 1 In a typical recipe, 2 g Ti 3 AlC 2 (325 mesh, Carbon-Ukraine) was gradually added to a mixture of 40 mL 9 M HCl and 2 g LiF at 35°C. After stirring for 24 h, the product was centrifuged and washed with deionized water until pH was above 5. Then, 40 mL of water was added to the sediment and vortexed for 30 min. The mixture was centrifuged at 3500 rpm for 30 min to obtain the supernatant containing few-layered Ti 3 C 2 T x and then stored at 4°C in N 2 atmosphere. The Ti 3 C 2 T x concentration was determined by weighting the dried self-standing film filtrated from a certain volume of the Ti 3 C 2 T x dispersion.
The Raman spectrum of Ti 3 C 2 T x was performed on the Horiba Jobin-Yvon confocal Raman system. For Raman sample preparation, the Ti 3 C 2 T x dispersion was first diluted to 0.1 mg mL −1 .
Then 5 µL of the diluted dispersion was drop-casted on a 1×1 cm 2 highly ordered pyrolytic graphite (HOPG) substrate pre-treated with 5 min of O 2 plasma. The Raman spectrum was acquired from few-layer Ti 3 C 2 T x flakes under 60 s exposure time of 15 mW 785 nm semiconductor laser with ×100 objective.

Resonant Raman spectra
Since Ti 3 C 2 T x is metallic, the measured Raman spectra is always under resonant conditions. The S3 which is in good agreement with experimental measurements (see Fig. 1

Comparison of distributions
To reduce the noise in the spectra shown in Figs. 4 and 5 of the main paper, we performed five MD runs for each concentration, with a different distribution of surface terminations in each run.

Lorentzian fitting of the peaks
In order to compare the frequencies of the peaks with experiments, we have to fit the peaks. We use a lorentzian function, as written in equation S1, where a is the intensity, ω is the frequency, and Raman intensity (arb. units) Figure S2: Raman spectra for x=0.5 concentration but different distributions of the surface terminations.

S5
G is the width of the peak. One of the interesting features of RGDOS is that the peaks can easily be isolated by only considering spectra projected to a particular eigenmode. Figs. S3 (a)-(c) show the E g (Ti) and A 1g (Ti) modes as well as the A 1g (C) mode for different concentrations of surface terminations. The Lorentzian function nicely fits the spectra and the evolution of the frequencies can easily be followed. Numerical values of the frequencies are reported in Table 1 of the main article.
(S1) Effect of multilayer and surface water.
To explain the wide peaks observed at 300-400 and 600 cm −1 , we investigated the effect of multilayer MXene and interaction with water molecules at the surface on the Raman spectra.

S6
MLFF model is first extended to multilayer MXene by adding bulk and bilayer configurations.
Bulk and bilayer structures were constructed from the monolayer and using a O-H-O chain to link the layers, which should be energetically the most favorable configuration at x=0.5. 3 Note that two different stacking configurations were studied for the bulk phase. Both converged to the ideal configuration where O-H-O are aligned between layers during the training MD. Production runs were therefore started from this stacking for bulk phase and multilayers. Results for different number of layers and bulk MXene are shown in Fig. S4(a). While the transition from monolayer to bulk shows clear changes, there are no sign of wide peaks. We conclude that multilayer is not responsible for the widening of peaks.
MLFF model is then extended to describe water molecules at the MXene surface. The model is trained by adding small amount of molecules (2 in the beginning) at the surface and progressively increasing this number (up to 10 molecules). Note that all of the water molecules are located on the same surface for both training and production MDs. Results for various number of molecules are shown in Fig. S4(b). Overall, water molecules seem to have fairly little impact on the Raman spectra and again do not explain the wide features observed in experiments.

Phonon density of states
To observe the effect of modes outside Γ-point on the Raman spectra, we study the phonon density of states obtained from velocity auto-correlation function (VACF).
The total phonon DOS from VACF is shown in Fig. S5(a), but cannot be directly compared to experiments since it also contains the Raman-inactive modes. Fig. S5(b) shows the VACF projected on the Raman active modes only and decomposed into contributions from each mode.
To better reproduce the Raman spectra, one can keep only the contribution from the Raman active modes and sum over the first Brillouin zone using the weights shown in Fig. 5