The influence of crystal thickness and interlayer interactions on the properties of heavy ion irradiated MoS2

Ion irradiation is a versatile tool to introduce controlled defects into two-dimensional (2D) MoS2 on account of its unique spatial resolution and plethora of ion types and energies available. In order to fully realise the potential of this technique, a holistic understanding of ion-induced defect production in 2D MoS2 crystals of different thicknesses is mandatory. X-ray photoelectron spectroscopy, electron diffraction and Raman spectroscopy show that thinner MoS2 crystals are more susceptible to radiation damage caused by 225 keV Xe+ ions. However, the rate of defect production in quadrilayer and bulk crystals is not significantly different under our experimental conditions. The rate at which S atoms are sputtered as a function of radiation exposure is considerably higher for monolayer MoS2, compared to bulk crystals, leading to MoO3 formation. P-doping of MoS2 is observed and attributed to the acceptor states introduced by vacancies and charge transfer interactions with adsorbed species. Moreover, the out-of-plane vibrational properties of irradiated MoS2 crystals are shown to be strongly thickness-dependent: in mono- and bilayer MoS2, the confinement of phonons by defects results in a blueshift of the A1g mode. Whereas, a redshift is observed in bulk crystals due to attenuation of the effective restoring forces acting on S atoms caused by vacancies in adjacent MoS2 layers. Consequently, the A1g frequency of tri- and quadrilayer crystals is statistically invariant on account oft competition between phonon confinement effects and interlayer interactions. The A1g linewidth is observed to decrease in bi- and trilayer crystals after low dose irradiation and is attributed to layer decoupling. This work shows that there is a complex interplay between defect production, crystal thickness and interlayer interactions in MoS2. Our results demonstrate that ion irradiation is an effective tool to modulate the electronic, vibrational and structural properties of MoS2, which may prove beneficial for practical applications.


Introduction
2D transition metal dichalcogenides have attracted intensive research efforts in recent years due to their diverse range of applications [1]. In particular, monolayer (1L) MoS 2 is a direct bandgap semiconductor [2] thus making it a promising candidate for electronic and optoelectronic devices [3][4][5].
Moreover, the mechanical properties of 2D MoS 2 crystals [6] facilitates their use in flexible and transparent electronics [7].
The controlled inclusion of defects such as vacancies and nanopores into 2D MoS 2 has proved to be an effective method to engineer its physical properties. For example, nanopores with tuneable diameters have been generated via electron beam drilling using a transmission electron microscope [8,9], showing the potential of the modified material for DNA sequencing and osmotic power generation applications [10][11][12]. Similar to electron irradiation, defectengineering of 2D crystals using ion beams has also been shown to offer unprecedented control over the type, size and concentration of the desired defects [13][14][15]. For example, metal and chalcogenide vacancies can be introduced into suspended 1L MoS 2 by 30 keV Ga + ions. The concentration and effective diameter of the defects were found to positively correlate with the ion fluence showing that defective crystals can be precisely engineered to study ionic transport and water desalination [16]. The electronic [17] and optical [18] properties of 1L MoS 2 have also been modulated by heavy ion irradiation. In both cases, the p-type changes observed were attributed to charge transfer interactions. Moreover, as the edges of MoS 2 crystals possess vacant coordination sites, this material has also shown great promise as a catalyst for the hydrogen evolution reaction (HER) [19,20]. Irradiation of 1L MoS 2 with 500 keV Au + ions has been shown to significantly improve its catalytic performance. The enhanced HER kinetics were positively correlated with the ion fluence and attributed to an increase in the number of catalytically active vacancies produced by elastic collisions [18].
To our knowledge, no study has systematically investigated the influence of crystal thickness on defect production in MoS 2 due to predominantly elastic collisions between the incident ion and the lattice atoms, i.e. in the nuclear energy loss regime. Mechanistic study within this energy regime is relevant for a broad community of scientists and technologists using ion implantators and focussed ion beam accelerators for defect-engineering. The momentum transferred to the recoil atoms during elastic collisions can displace them from their lattice sites and potentially eject (sputter) them from the surface of the crystal. However, if the recoil atom is not sputtered from the crystal, vacancy defects may still be produced; providing that the kinetic energy of the recoil atom is greater than its threshold displacement energy (T d ).
For the same total energy deposited by 2 MeV H + (1.26 x 10 16 ions cm −2 ) and 390 keV He + (1 × 10 15 ions cm −2 ) irradiations, He + irradiated 2D MoS 2 field effect transistors exhibit greater degradation of device performance characteristics [21]. This is attributed to a higher concentration of vacancy defects and interface trapped charges in He + irradiated devices on account of the larger nuclear stopping power of 390 keV He + projectiles. Therefore, it is evident that the concentration of defects produced by elastic collisions in MoS 2 crystals is a key parameter which directly modulates device performance. In light of this consideration, the parameter by which we define radiation hardness in this study is the concentration of stable defects formed as a function of fluence i.e. the degree of structural disorder introduced by 225 keV Xe + irradiation.
In the specific case of nanomaterials, due to their high surface area to volume ratio, the yield of sputtered atoms can be significantly higher than in bulk crystals [13], thus implying that nanomaterials may be less radiation hard than the corresponding bulk compounds. Conversely, for ions in the keV range, the cross section for direct defect production (i.e. elastic collisions involving incident ions) in nanomaterials decreases significantly with increasing incident ion energy [22]. A consequence of this is that the stopping power of the incident ion is effectively constant throughout its trajectory within the nanomaterial as the ion deposits only a fraction of its total kinetic energy [13]. Whereas, in bulk crystals the entirety of this kinetic energy is transferred to the lattice and the stopping power of the ion varies along its trajectory or 'ion track' . Specifically, the highest concentration of defects in a bulk target will be localised around the Bragg peak where the nuclear stopping power of the heavy ion increases. The majority of these defects are created by energetic recoil atoms produced in displacement cascades. Consequently, due to the abundance of potential recoil atoms, bulk targets could be less radiation hard than nanomaterials due to their lack of dimensional confinement.
Herein represents a systematic study of defect production in single-, few-layer and bulk crystals of MoS 2 under 225 keV Xe + irradiation in order to obtain fundamental insights into the radiation hardness and defect-engineered properties of this material as a function of crystal thickness. Our results show that MoS 2 becomes more susceptible to heavy ion-induced radiation damage as the crystal thickness decreases from quadrilayer to monolayer, despite the uniform energy deposition profile of the incident Xe + ions throughout the nanostructures. However, the concentration of defects in irradiated bulk MoS 2 are not significantly lower than those of quadrilayer crystals, under our experimental conditions. This is attributed to the contribution of displacement cascades, involving Mo and S recoil atoms, towards indirect defect production in bulk MoS 2 . The out-of-plane vibrational properties of defect-engineered MoS 2 crystals are also shown to be strongly thickness-dependant. This is attributed to competition between attenuated interlayer interactions and defect-induced reductions in phonon correlation lengths i.e. phonon confinement. Our results show that heavy ion irradiation can be used to induce p-doping in MoS 2 and clearly demonstrate the influence of crystal thickness on the radiation stability of 2D MoS 2 . The established mechanisms can be extended to a wide range of layered crystals.

Results and discussion
MoS 2 crystals of various thicknesses and 1L crystals produced by both micromechanical exfoliation (MME) [23] and chemical vapour deposition (CVD) [24] were bombarded with 225 keV Xe + ions (irradiation conditions and key dosimetric parameters are detailed in the methods section and table S1, stacks.iop.org/TDM/7/035011/mmedia). Samples produced by MME contain mono-, bi-, tri-, quadrilayer and bulk crystals as evidenced by correlative Raman spectroscopy and atomic force microscopy (figure S1, supplementary information). In contrast, samples produced by CVD contain mostly 1L single crystals as confirmed by Raman spectroscopy. X-ray photoelectron spectroscopy (XPS) is used to evaluate the changes in the chemical composition and electronic properties of defective MoS 2 , whilst Raman spectroscopy and electron diffraction are utilised to gain insights into the radiation-induced changes in the vibrational properties and crystallinity of irradiated MoS 2 , respectively. Figure 1 shows the XPS analysis of pristine and irradiated MoS 2 samples produced by MME and CVD up to a fluence of~2 × 10 14 ions cm −2 . Deconvoluted XPS spectra are provided in figures S2-S5, supplementary information. The majority of the crystals deposited onto the substrate during MME are bulk [25]. Therefore, as the diameter of the x-ray beam is on the order of a few hundred micrometres, the XPS spectra of MoS 2 produced by MME and CVD is representative of bulk and 1L crystals, respectively. Figure 1(a) shows that the S:Mo atomic ratio of both 1L and bulk MoS 2 decreases upon exposure to radiation due to sputtering of lattice atoms by the impinging Xe + ions. In addition, as the MoS 2 crystals are deposited onto SiO 2 /Si wafers, vacancies can also be produced by elastic collisions between substrate recoil atoms and Mo/S lattice sites [26,27]. The observed decrease in the S:Mo atomic ratio indicates that S atoms are preferentially sputtered with respect to Mo atoms. The reason for this is threefold: firstly, the stoichiometry of MoS 2 dictates that collisions involving S atoms are statistically more probable. Secondly, S atoms have a larger cross section for elastic collisions involving 225 keV Xe + ions, compared to Mo atoms [22]. Thirdly, the T d of S atoms is around 15-25 eV lower than that of Mo atoms [22,28], due to differences in their bonding environments.
Despite the preferential sputtering of S atoms, simulations by Kretschmer et al show that Mo atoms can be sputtered from 1L MoS 2 on a SiO 2 substrate under 225 keV Ar + irradiation [27]. Hence, as sputtering yields increase with ion mass, Mo vacancy production in MoS 2 /SiO 2 samples irradiated with 225 keV Xe + is highly likely. For instance, the nuclear stopping power for Xe + ions travelling in MoS 2 exhibits a maximum value of 332.4 eV Å −1 at 225 keV (figure S7, supplementary information) which is more than six times that of 225 keV Ar + ions (53.2 eV Å −1 ). Indeed, Monte Carlo calculations confirm that Mo vacancies are produced in bulk MoS 2 under 225 keV Xe + irradiation (figures S8-S10, supplementary information).
The S:Mo atomic ratios of irradiated 1L MoS 2 are significantly lower than those of bulk crystals at similar fluences ( figure 1(a)). Hence, the rate at which S atoms are sputtered as a function of radiation exposure is considerably higher for 1L MoS 2 . The significantly smaller decrease in the S:Mo atomic ratio of bulk MoS 2 is in agreement with previous ion irradiation studies of this material [29] and could be attributed to the smaller fraction of surface atoms in bulk crystals relative to 2D MoS 2 . More specifically, molecular dynamics (MD) and Monte Carlo simulations have shown that the sputtering of S atoms from the bottom face of 1L MoS 2 on a SiO 2 substrate is significantly suppressed, compared to suspended crystals [22,27]. The substrate is modelled as an external repulsive potential acting on the Mo/S atoms which effectively 'stops' the recoil atom from being sputtered from the crystal and could promote vacancy-interstitial recombination [27]. Drawing upon these considerations, one might argue that in bulk MoS 2 , when a vacancy is created, the lattice atoms in adjacent layers could act analogously to the modelled substrate i.e. the recoil atom experiences a repulsive potential and is 'stopped' . This effect of neighbouring MoS 2 layers could suppress sputtering yields of non-surface atoms and account for the slower decay of the S:Mo atomic ratio in bulk crystals.
A significant increase in the Mo VI :Mo IV atomic ratio is observed for 1L MoS 2 upon irradiation (figure 1(b)). The increase in the intensity of the Mo VI 3d doublet is attributed to the formation of MoO 3 as previously reported for 1L MoS 2 under swift heavy ion irradiation at grazing incidence [30]. However, it is possible that S vacancy point defects could yield Mo VI centres with a mixed bonding environment involving both S and O atoms. The observation of Mo IV oxidation is in agreement with density functional theory calculations showing that the activation energy barrier for O 2 dissociation is lower when the molecule is adsorbed to a vacancy in the MoS 2 lattice as opposed to physisorbed on the pristine basal plane of the crystal [31]. We propose that Mo IV oxidation may proceed via adsorption of O 2 to radiationinduced chalcogen vacancies followed by dissociation and Mo-O bond formation to yield MoO 3 . Moreover, XPS spectra of highly defective bulk MoS 2 indicates the presence of MoO 2 or Mo IV S y O x species, which could be attributed to amorphisation of the layered crystal structure (figures S4 and S6, supplementary information).
In order to evaluate the effect of heavy ion irradiation on the electronic properties of MoS 2 , the binding energies of the Mo IV 3d 5/2 (figure 1(c)) and S 2p 3/2 (figure 1(d)) core levels are plotted as a function of fluence. The binding energies of both photoelectron lines are observed to decrease as defects are introduced into 1L and bulk MoS 2 , in agreement with a previous study on bulk MoS 2 using 2 keV Ar + ions [29]. During XPS measurements, the samples are in Ohmic contact with the spectrometer ground, hence, the Fermi levels of the sample and the spectrometer are aligned. The binding energies of core level electrons are measured relative to the Fermi level of the system and the work function is defined by the vacuum level of the spectrometer as opposed to that of the sample.
The reduction in the binding energies of the Mo IV 3d 5/2 and S 2p 3/2 photoelectrons in MoS 2 infers that the Fermi level has been lowered i.e. MoS 2 becomes p-doped. In regard to the physical origin of the binding energy decease, when a p-type material is in contact with the spectrometer ground, electrons are transferred to the material until their electrochemical potentials are equalised. This charging of the p-type material increases the kinetic energy of the outgoing photoelectrons i.e. the binding energy of the photoelectron lines of p-type materials are lower than their undoped counterparts [32].
This observation of p-type doping in irradiated MoS 2 can be attributed to the presence of Mo vacancies, which produce acceptor states in the band structure of MoS 2 [33]. Indeed, the binding energies of Mo 3d and S 2p photoelectrons lines are observed to decrease in p-type crystals exhibiting S:Mo ratios of 2.3:1 which could be attributed to Mo vacancies [34]. Furthermore, the interaction of O 2 can also play an important role: despite the activation energy barrier for O 2 dissociation being lower when adsorbed to vacancies [31], some O 2 molecules are likely to remain adsorbed on defective crystals without subsequent dissociation. The charge transfer between O 2 and MoS 2 is more efficient when O 2 molecules are adsorbed to vacancies, as opposed to physisorbed on the basal plane of MoS 2 [35], leading to p-doping [18,36]. Furthermore, the adsorption of other impurities onto vacancy sites could promote p-type doping of irradiated MoS 2 . For instance, in all samples, C 1s photoelectron lines are observed at~286 eV and~288.5 eV corresponding to the C-O and C=O bonds of adventitious carbon. It follows that if such moieties are present on the surface of the irradiated crystals then a charge transfer mechanism similar to that of O 2 may contribute to the observed p-type doping of MoS 2 .
Therefore, both the intrinsic properties of vacancies, i.e. acceptor states, and their interaction with adsorbed species, i.e. charge transfer, are likely to account for the p-type behaviour observed in irradiated 1L and bulk MoS 2 . Previously, it has been shown that obtaining p-type MoS 2 via electron beam irradiation is not possible due to the formation of metallic Mo nanograins [37]. However, in this work, by exploiting the momentum transfer during elastic collisions involving heavy ions, figures 1(c), (d) and Monte Carlo calculations (figures S8-S10, supplementary information) suggest that Mo atoms must also be sputtered from the lattice to suppress nanograin formation yielding p-type MoS 2 . Figure 1(b) shows that 1L MoS 2 irradiated to a fluence of 1.79 × 10 14 ions cm −2 exhibits a Mo VI :Mo IV atomic ratio of 0.94 which agrees well with the value reported by Madauß et al after irradiation of 1L MoS 2 with swift heavy ions at grazing incidence [30]. However, these authors also state that the Mo IV :S ratio remains at 0.45 before and after irradiation. Whereas, in this work, an increase from 0.48 to 0.54 for the sample irradiated to 1.79 × 10 14 ions cm −2 is observed. This result suggests that not all Mo IV metal centres bearing S vacancies are oxidised to Mo VI to form MoO 3 and further supports the charge transfer doping mechanisms outlined.
Figures 1(c) and (d) show that the Mo IV 3d 5/2 and S 2p 3/2 binding energies in bulk crystals remain invariant between~5 × 10 13 and~2 × 10 14 ions cm −2 . This finding is in agreement with previous studies as a relatively small number of defects are able to effectively 'pin' the Fermi level [29,38]. As previously outlined, the rate of S atom sputtering and defect production is higher in 1L MoS 2 under heavy ion irradiation ( figure 1(a)). Therefore, one would expect the Fermi level of 1L MoS 2 to be pinned at a lower fluence relative to the bulk crystal, namely below~5 × 10 13 ions cm −2 . However, as 1L MoS 2 is deposited onto a SiO 2 substrate, charge transfer at this interface forbids direct correlation of binding energy shifts with Fermi level pinning in 1L MoS 2 /SiO 2 samples as defects are introduced into both 1L MoS 2 and SiO 2 . Finally, valence band XPS spectra show a decrease in the separation between valence band maximum and Fermi level for bulk MoS 2 irradiated to a fluence of 5 × 10 13 ions cm −2 (figure S11, supplementary information) which further corroborates the existence of radiation-induced p-type doping in MoS 2 [29,39].
The typical Raman spectrum of pristine 1L MoS 2 shows two peaks at 385.5 cm −1 and 403.5 cm −1 , attributed to the in-plane (E ′ ) and out-of-plane (A ′ 1 ) normal vibrational modes, respectively. In the case of irradiated MoS 2 , a new peak appears: the LA (M) mode is a defect-activated peak akin to the D peak in graphene [40] and an analogy has been drawn between the LA (M) /A ′ 1 intensity (I) ratio in defective 1L MoS 2 and the I D /I G ratio in graphene [41].
In pristine 2D crystals, first-order Raman active phonons can be regarded as plane waves with wavevector q ≈ 0 i.e. they originate from the centre of the Brillouin zone. In the case of nanoparticles and nanoribbons, the phonons cannot propagate beyond the boundary of the nanostructure and can be considered as 'confined' . Due to the Heisenberg uncertainty principle, the fundamental selection rule q = 0 is relaxed in nanostructures [42] and phonon uncertainty is described by ∆q = 1/L C where L C represents the phonon correlation length [43]. In pristine 2D crystals, L C is effectively infinite, whereas, in the case of nanoparticles and nanoribbons, L C corresponds to the diameter or width of the nanostructure, respectively. Similarly to crystal boundaries, vacancies and defects disrupt the uniform propagation of phonons leading to an increase in the phonon attenuation probability and a concurrent reduction of L C . As ion fluences (σ) are conventionally expressed in ions cm −2 , one can estimate the average defect spacing (L D ) in irradiated crystals as 1 / √ σ and phonon uncertainty can be approximated as ∆q = 1/L D . Therefore, L C decreases and ∆q increases in heavy ion-irradiated crystals. Consequently, for first-order phonons, q increases as a function of fluence. Hence, the phonon frequency will change in accordance with the gradient of the corresponding phonon branch as q increases from the centre of the Brillouin zone. The Raman spectra of numerous nanostructures and irradiated crystals have been interpreted using this framework; known as the phonon confinement or spatial correlation model [41][42][43][44][45]. Figure 2 shows the frequency (ν) of the A ′ 1 (figure 2(a)) and E ′ (figure 2(b)) modes of 1L MoS 2 produced by MME and CVD. The ν A ′ 1 and ν E ′ modes monotonically blueshift and redshift, respectively, with increasing radiation exposure. This observation is in agreement with the phonon dispersion of MoS 2 as the gradients of the A ′ 1 and E ′ phonon branches are positive and negative, respectively, as q increases from the centre of the Brillouin zone; in accordance with the spatial correlation model [41]. Figure 2(c), shows the linewidth (Γ) of the A ′ 1 mode as a function of radiation exposure. When defects are introduced into the 1L MoS 2 lattice, they disrupt the propagation of the phonon and reduce its lifetime. Hence, as the Γ of a Raman signal is inversely proportional to the phonon lifetime, a linear increase in the Γ A ′ 1 mode (figure 2(c)) is observed due to a reduction in the phonon lifetime caused by the formation of defects. Figure 2(d) shows the I LA(M) /I A ′ 1 ratio of 1L MoS 2 as a function of radiation exposure. This parameter increases linearly as a function of increasing defect density, which agrees with the linear dependence of ν A ′ 1 and Γ A ′ 1 mode observed in figures 2(a) and (c). Note that unlike the ν A ′ 1 mode, which blueshifts linearly with increasing radiation exposure, the ν E ′ mode redshifts until a fluence of~1 × 10 14 ions cm −2 before becoming invariant at higher fluences (figure 2(b)). The Γ E ′ mode increases linearly above this fluence (figure S12(a), supplementary information) indicating that the lifetime of these in-plane phonons is still being progressively reduced by the increasing defect density. The physical mechanism underpinning this frequency invariance is the subject of further research. However, potential contributing mechanisms may include radiation-induced compressive strain [18,46] and significant changes to the phonon dispersion in highly defective 1L crystals. Indeed, a clear broadening of the dispersion curves is observed in MD simulations of defective MoS 2 with mono-S and mono-Mo vacancies at concentrations as low as 0.5% [47]. Whereas, in 1L MoS 2 after irradiation to a fluence of ≥1 × 10 14 ions cm −2 , S:Mo atomic ratios as low as ∼1 are observed; indicating significantly higher defect concentrations ( figure 1(a)). Moreover, at fluences ≥1.57 × 10 14 ions cm −2 the A ′ 1 frequency is blueshifted relative to the linear trend observed at lower fluences; predicted by the spatial correlation model. This could be attributed to higher vacancy concentration and consequent p-type doping; on account of the strong electron-phonon coupling of this vibration [48,49].
Returning to the XPS analysis of irradiated 1L MoS 2 crystals, despite the invariance of the S:Mo and Mo VI :Mo IV ratios between~5 × 10 13 and~2 × 10 14 ions cm −2 (figures 1(a) and (b)), Raman spectroscopy shows that the 1L MoS 2 crystals become progressively more defective within this fluence region ( figure 2(d)). Correlation of these two techniques suggests that the rate at which Mo and S atoms are sputtered become comparable in highly defective 1L MoS 2 crystals. Density functional theory calculations have shown that the T d of S atoms at the edge of a WS 2 nanoribbon are~2.8 eV lower than those in the central region; on account of their different bonding environment [28]. Therefore, changes to the T d of Mo and S atoms could partly account for this comparable sputtering rate due to the significant alteration of the Mo and S bonding environments in defective crystals irradiated to fluences ≥ 5 × 10 13 ions cm −2 . In addition, significant changes to the stoichiometry of the MoS 2 crystals are also likely to contribute to the comparable rate of S and Mo atom sputtering. For example, after irradiation to a fluence of 5 × 10 13 ions cm −2 , the S:Mo atomic ratio of~1 dictates that the probability of elastic collisions between incident ions/recoils atoms and the Mo and S lattice atoms become comparable.
Moreover, despite the considerable concentration of Mo VI centres in the irradiated 1L MoS 2 crystals ( figure 1(b)), no new Raman signals characteristic of MoO 3 are observed. This is in agreement with the expected distribution of disordered oxidised Mo VI regions, as opposed to ordered crystalline MoO 3 domains, as the ion beam is rastered homogeneously over the samples. Within the framework of the spatial correlation model, substitutional O atoms and vacancy point defects are indistinguishable as they both serve to increase the probability of phonon attenuation and decrease L c . Therefore, the morphology of the irradiated samples is likely to consist of MoS 2 single crystals interspersed with highly disordered MoO 3 domains and various types of vacancy point defects such as single Mo/S vacancies and S divacancies. However, the presence of vacancies with larger atomic configurations such as V 2Mo + 6S and V 3Mo + 5S cannot be ruled out [16]. Figure 3 shows the evolution of the vibrational properties of 1L, bi-(2L), tri-(3L), quadri-layer (4L) and bulk MoS 2 crystals as a function of radiation exposure. It should be noted that the Mulliken symbols for the out-of-plane and in-plane modes change to A 1g and E 1 2g in 2L MoS 2 on account of the different point groups of MoS 2 crystals with an odd or even numbers of layers [50] and will be used henceforth. Figure 3(a) shows the ν A1g mode for a given thickness of MoS 2 relative to its value prior to irradiation (∆ν A1g ). The corresponding data for the E 1 2g in-plane vibrational mode is available in figure S13, supplementary information. A clear difference is seen depending on the crystal thickness: the ν A1g mode of 1L and 2L MoS 2 progressively blueshifts with increasing fluence i.e. defect concentration, as previously discussed. However, the ν A1g mode of 3L and 4L does not change notably with increasing defect concentration whilst, for the bulk crystal, a redshift is observed. As the A 1g mode involves the out-of-plane vibration of S atoms, adjacent MoS 2 layers within bulk crystals exert an effective restoring force on one another. Consequently, prior to irradiation, the ν A1g mode is highest in bulk crystals due to the increased force constant (figure S14(a), supplementary information) [51]. When S vacancies are introduced into the lattice via irradiation, the effective restoring forces acting on the chalcogen atoms in the adjacent layers are attenuated, which leads to the observed redshift. Therefore, the confinement of phonons by defects is not the dominant contribution to the out-of-plane vibrational properties of defective bulk MoS 2 ; instead, these properties are strongly modulated by interlayer interactions. Note that, in pristine 1L MoS 2 , no chalcogen atoms can participate in any interlayer interactions and the ν A ′ 1 mode exhibits a value of 403.5 cm −1 ( figure 2(a)). This energy represents the limit to which the ν A1g mode can redshift in defective bulk MoS 2 due to reduced interlayer interactions. We observe that even in highly defective bulk MoS 2 crystals, irradiated to fluences up to 1 × 10 16 ions cm −2 , the ν A1g mode does not decay below value (figure S26(a), supplementary information). Additional Raman and XPS data of highly defective bulk MoS 2 , matching the three stage damage accumulation model [52][53][54], are available in figures S26 to S30, supplementary information.
The competition between (1) the reduction of L c caused by phonon confinement and (2) attenuated interlayer interactions is evidenced by the ∆ν A1g mode for 3L and 4L MoS 2 with increasing defect concentration ( figure 3(a)). In defective 3L and 4L MoS 2 , the confinement of the out-of-plane phonons by the vacancies in the nanostructure acts to increase their energy. However, even within pristine MoS 2 crystals just 2.1 nm (3L) and 2.8 nm (4L) thick, this corresponds to 67% and 75% of S atoms participating in interlayer interactions, respectively. Therefore, the attenuation of the effective restoring forces by vacancies mitigates their confining effects on phonons; hence, the invariance of the ν A1g mode in 3L and 4L MoS 2 with increasing defect density exemplifies the reciprocity between these two mechanisms at the mesoscopic scale.
Our results appear to differ from findings of the recent work by Cheng et al, who reported a redshift of the A 1g mode for all MoS 2 crystal thicknesses greater than 1L under 60 eV and 200 eV Ar + irradiation [55]. The cross section for S vacancy production in freestanding monolayer MoS 2 caused by 60 eV Ar + ions (~7 Å 2 ) is twice as large as 225 keV Xe + ions (~3.5 Å 2 ) [22]. Consequently, the number of S atoms sputtered by 60 eV Ar + ions is also twice that of 225 keV Xe + ions. Therefore, one would expect higher levels of structural disorder in 60 eV Ar + irradiated freestanding MoS 2 crystals compared to 255 keV Xe + irradiated material at a similar fluence. However, the energy of the recoil atoms originating from the SiO 2 substrates, used in both studies, will be quite different and could influence sputtering yields. Moreover, our irradiation experiment is carried out at normal incidence, whilst Cheng et al utilise an off-axis convergent ion beam. Based on MD simulations [22], the angle of incidence is expected to strongly influence the ν A1g mode of irradiated 2L MoS 2 if S atoms participating in interlayer interactions are preferentially sputtered. These findings demonstrate that through careful consideration of the irradiation conditions, the vibrational properties of defective few-layer MoS 2 crystals can be precisely engineered. Figure 3(b) shows the effect of radiation exposure on the Γ A1g mode for a given thickness of MoS 2 relative to its value prior to irradiation (∆Γ A1g ). As expected, the Γ A1g mode increases as defects are introduced into the crystal. However, a significant decrease in the Γ A1g mode is observed in 2L and 3L MoS 2 after irradiation to the lowest fluence of 2 × 10 13 ions cm −2 before progressively increasing at higher fluences. In non-irradiated MoS 2 , the Γ A1g mode of 1L MoS 2 increases from 5.3 cm −1 to 9 cm −1 for 2L crystals; this value then steadily decreases with increasing thickness until a value of 3.2 cm −1 for bulk material ( figure S14(b), supplementary information). Therefore, the Γ A1g mode is characteristic of the thickness of 2L and 3L MoS 2 [56]. The significant decrease in the Γ A1g mode after irradiation to a fluence of 2 × 10 13 ions cm −2 could be attributed to radiationinduced layer decoupling of the crystals, such that they exhibit '1L-like' Γ A1g values. An increase in the interlayer spacing of just 1 Å is sufficient to induce an indirect-direct bandgap transition in 2L MoS 2 [57]. Moreover, decoupling of few-layer MoS 2 has been observed experimentally under oxygen plasma [57] and 100 keV proton irradiation [58]. Similarly in this work, elastic collisions involving heavy ions could also induce layer decoupling.
Further evidence of potential radiation-induced layer decoupling in 2L and 3L MoS 2 is provided in figure 3(c), where the I E 1 2g /I A1g ratio in crystals of various thicknesses is plotted as a function of radiation exposure. In non-irradiated material, this ratio exhibits values of 2.1 and 1.4 for 2L and 3L MoS 2 , respectively, when using the peak heights as opposed to integrated intensities [56]. These values are characteristic of 2L and 3L crystals as they are mutually exclusive and distinct from all other MoS 2 thicknesses. After irradiation to the lowest fluence of 2 × 10 13 ions cm −2 the I E 1 2g /I A1g ratio of both 2L and 3L MoS 2 decreases significantly to values close to the defective 1L sample, which may be attributed to layer decoupling. As mentioned previously, the energy deposition profile (stopping power) of the incident ion is effectively constant throughout the 1L -4L MoS 2 crystals as the lattice atoms are constrained to within the first 0.7-2.8 nm of the incident ion track, respectively. Moreover, as the incident Xe + ions deposit only a fraction of their total kinetic energy within the 2D MoS 2 crystals, the typical flux and kinetic energies and of the Si and O recoil atoms that they are exposed to from the underlying SiO 2 substrate will be equivalent. On this basis, one might expect comparable concentrations of defects in 1L -4L MoS 2 crystals.  Figure 3(d) shows that the rate of defect production in MoS 2 is inversely proportional to the thickness of the crystal i.e. thinner 2D MoS 2 crystals are more easily damaged by heavy ion irradiation, in agreement with XPS data ( figure  1(a)). Despite the uniformity of the energy deposition profile, the neighbouring MoS 2 layers in thicker 2D MoS 2 crystals exert a repulsive potential on recoil atoms. This effectively 'stops' the recoil atom from being sputtered from the surface of the crystal and creating a permanent defect. In 1L MoS 2 , a larger fraction of recoil atoms possess sufficient kinetic energy to be sputtered from the surface of the crystal or displaced into the underlying SiO 2 substrate. Whereas, in thicker 2D crystals, the velocity of the recoil atom is reduced as it interacts with adjacent MoS 2 layers; facilitating vacancy-interstitial recombination.
Considering this stabilising effect of neighbouring MoS 2 layers, interestingly, figure 3(d) shows that the rate of defect production in bulk MoS 2 is not significantly lower than 4L crystals. This may be rationalised by considering the penetration depth of the 514.5 nm laser into bulk MoS 2 crystals. Optical absorption spectroscopy of MoS 2 shows an approximately linear increase in absorbance as the thickness of the crystal is gradually increased from 1L to 6L [59]. Extrapolation of these absorbance values and subsequent conversion to transmittance suggests that after a penetration depth of 8.4 nm the transmittance of the 514.5 nm laser falls below 10%; when accounting for the backscattering geometry of the spectrometer. Figure 4 shows the vacancy production and damage accumulation results from Monte Carlo calculations of bulk MoS 2 under 225 keV Xe + irradiation. Details of the calculations are available in the methods section. Additional results are available in figures S8-S10, supplementary information, showing that damage accumulation in bulk MoS 2 is qualitatively comparable regardless of the choice of reported T d values used [22,28]. Figure 4(a) shows the number of vacancies produced ion −1 Å −1 as a function of depth in bulk MoS 2 crystals. The number of S vacancies (V S ), Mo vacancies (V Mo ) and total vacancies (V total ) account for both direct and indirect defect production i.e. displacements of lattice atoms caused elastic collisions involving Xe + ions (direct) and Mo/S recoil atoms (indirect). The shaded region between 0 nm and 8.4 nm in figure 4(a) represents the penetration depth of the 514.5 nm laser used to acquire Raman spectra, as described earlier. It is evident that the high concentration of vacancies caused by displacement cascades contribute to the Raman spectrum of irradiated bulk MoS 2 ; even if the most defective region at a depth of 38 nm does not contribute significantly. As previously mentioned, in bulk targets, the highest concentration of defects will be localised around the Bragg peak. Hence, as the nuclear stopping power of Xe + ions traveling in MoS 2 is largest at 225 keV (figure S7, supplementary information) a sharp increase in vacancy concentration is observed between 0 nm-8.4 nm. Therefore, administering a 225 keV Xe + beam promotes both direct and indirect defect production within the surface region of bulk MoS 2 . However, it is the contribution from indirect defect production in bulk MoS 2 which accounts for its comparable rate of defect formation with respect to 4L crystals ( figure 3(d)).
As outlined previously, Si and O recoil atoms from the underlying SiO 2 substrate will also contribute to indirect defect production in 1L -4L MoS 2 crystals. However, the contribution of Mo and S recoil atoms to indirect defect production in bulk MoS 2 is more significant. The reason for this is two-fold: firstly, the nuclear stopping power of 225 keV Xe + ions in SiO 2 is 22% smaller than in MoS 2 (260.8 eV Å −1 verses 332.4 eV Å −1 ) and the T d of Si (70.5 eV) and O (28.9 eV) atoms in SiO 2 are considerably higher than Mo (20 eV-31.7 eV) and S (5 eV-6.9 eV) atoms in MoS 2 [22,28,60]. Hence, fewer recoil atoms are produced in SiO 2 under 225 keV Xe + irradiation. Secondly, when Si and O recoil atoms engage in elastic collisions with Mo and S lattice atoms, the minimum kinetic energy required for defect production (E min ) is higher compared to elastic collisions involving Mo and S recoil atoms. This is evidenced by the kinematic factor describing a head-on binary collision E min = T d (m i + m t ) 2 /4m i m t , were T d is the threshold displacement energy and m i (m t ) represent the ion (target atom) masses [22]. E min occurs when the masses of the ion and target atom are equal. Hence, a larger fraction of Mo and S recoil atoms have sufficient kinetic energy to displace lattice atoms in MoS 2 compared with Si and O recoil atoms. However, Monte Carlo calculations were performed using the Stopping Range/TRansport of Ions in Matter (SRIM/TRIM) software [61] and results must be interpreted with knowledge of three limitations. Firstly, the calculations do not account for thermal diffusion of target atoms. Hence, in situ defect annealing/vacancy-interstitial recombination will lead to a lower concentration of stable defects. Secondly, the calculations assume an amorphous target i.e. ion-channelling is not accounted for. The combined effect of these two limitations means that the distributions shown in figure 4(a) will be broader and the highest vacancy concentrations are likely to lower and shifted towards greater depths. Finally, damage accumulation is not accounted for in the calculations i.e. each ion track is simulated a pristine target. It is this assumption that gives rise to the strong linearity of the trends in figure 4(b). Figure 4(b) shows the results of damage accumulation calculations for the 0 nm-8.4 nm surface region of bulk MoS 2 , probed by Raman spectroscopy, in terms of displacements per atom (dpa) as a function of fluence for each of the aforementioned defect quantities. In line with their lower T d , the majority of displacements involve S atoms which is in agreement with XPS data ( figure 1(a)). However, the linear V s dpa trend cannot be directly translated into experimentally observed vacancy concentrations due to the omission of damage accumulation from the calculation i.e. changes in MoS 2 stoichiometry due to sputtering are not accounted for. The calculations afford considerable V Mo dpa values, hence, Mo vacancy production in heavy ion irradiated MoS 2 is highly likely and supports the acceptor state p-doping mechanism outlined earlier (figures 1(c) and (d)). Due to the limitations outlined, damage accumulation in the MoS 2 surface region caused by 225 keV Xe + irradiation will be smaller than that shown in figure 4(b) i.e. not every displaced atom will create a stable defect due to vacancy-interstitial recombination.
Ideally, controlled defect-engineering should preserve the crystallinity of 2D MoS 2 if the modified nanomaterial is intended for use in practical applications. Electron microscopy characterisation was performed in order to assess the crystallinity of pristine and heavy ion irradiated few-layer MoS 2 . Figure 5 shows the results of transmission electron microscopy (TEM) and electron diffraction measurements (additional bright field TEM images and electron diffraction patterns of irradiated MoS 2 crystals are available in figure S31, supplementary information). The electron diffraction pattern of non-irradiated MoS 2 ( figure 5(b)) exhibits the expected hexagonal pattern characteristic of the 2H polymorph of MoS 2 . The diffraction spots are well defined and have an appreciable intensity even for high order reflections which are indicative of significant long range structural order within the crystal. In few-layer MoS 2 , irradiated to a fluence of 9.86 × 10 15 ions cm −2 , broadening of the diffraction spots is observed ( figure 5(d)). Evaluating the diffraction pattern of a thinner region of the same irradiated flake (figure 5(e)), two striking differences are apparent. Firstly, there is a significant reduction in the intensity of higher order reflections indicating a greater loss of long range structural order. Secondly, the broadening of the diffraction spots is more pronounced in the direction perpendicular to the reciprocal lattice vector. The broadening of the diffraction spots is such that they begin to represent the diffuse rings of amorphous material. Similar observations have been reported for graphene during TEM experiments with in situ 30 keV He ion irradiation [62].
The amorphous carbon support film is between 15 nm and 20 nm thick, hence, slight broadening of the diffracted electron beam should be considered. Figure 5(c), shows that the thicker region of the irradiated few-layer MoS 2 crystal is supported by the carbon film (indicated by the blue arrow), whereas, the thinner region is suspended over a hole in the film (indicated by the red arrow). The diffraction spots originating from the thinner freestanding region are broader and more diffuse than those originating from the thicker supported region. Hence, the influence of the amorphous carbon support film is smaller than the effect of higher levels of structural disorder in thinner irradiated crystals. Moreover, as the irradiated crystals were transferred onto TEM grids after irradiation, the role of substrate Si and O recoils atoms on defect formation in the 2D MoS 2 crystals is the same for TEM, Raman and XPS samples.
Furthermore, at relatively low fluences of 5.92 × 10 13 ions cm −2 , electron diffraction measurements show that few-layer MoS 2 remains highly crystalline ( figure S31, supplementary information). Note that at these fluences MoS 2 is observed to be pdoped (figures 1(c) and (d)). Hence, we conclude that heavy ion irradiation is an effective tool for defectengineering of highly crystalline p-type MoS 2 .
The greater reduction in crystallinity observed for thinner MoS 2 crystals after heavy ion irradiation is corroborated by the XPS and Raman spectroscopy data: as the rate of defect production in MoS 2 is inversely proportional to the thickness of the crystal (figure 3(d)), this leads to higher sputtering yields in thinner MoS 2 crystals compared to bulk material ( figure 1(a)). Therefore, the rate at which structural disorder is introduced into the lattice during irradiation is higher for thinner crystals of MoS 2 (figures 5(d) and (e)).

Conclusion
This work represents a systematic analysis of the structural changes in MoS 2 under 225 keV Xe + irradiation. Our results show that the radiation hardness of MoS 2 is strongly influenced by the thickness of the crystal. Specifically, the rate of defect production is inversely proportional to the thickness of the crystal as evidenced by XPS, Raman and electron diffraction data. The out-of-plane vibrational properties of irradiated MoS 2 are progressively governed by interlayer interactions as the thickness of the crystal increases: the A 1g frequency blueshifts in monolayer and bilayer MoS 2 due to a reduction in the phonon correlation length caused by defect formation, whilst a redshift is observed in bulk MoS 2 due to attenuation of the effective restoring forces acting on S atoms caused by the formation of S vacancies in adjacent MoS 2 layers. Consequently, the A 1g frequency remains statistically invariant as a function of radiation exposure due to the reciprocity between these two competing mechanisms i.e. phonon confinement and attenuated interlayer interactions. We also observed that irradiated MoS 2 becomes p-doped whilst remaining highly crystalline. This observation further substantiates the use of heavy ion irradiation as a versatile route to modulate the properties of MoS 2 via defect-engineering. Our findings offer valuable insights into the radiation damage mechanisms of MoS 2 that can be extended to other semiconducting transition metal dichalcogenides and layered crystals.

Methods
Micromechanical Exfoliation: Mono-, bi-, tri-, quadri-layer and bulk MoS 2 samples were prepared via mechanical cleavage of bulk molybdenite (Manchester Nanomaterials) using dicing tape. The exfoliated flakes were deposited onto SiO 2 /Si substrates (IDB Technologies, oxide thickness ca. 300 nm) that had been heated to 130 • C in ambient conditions to suppress interfacial water inclusion. Optical microscopy was used to identify the crystals.
Chemical Vapour Deposition: Monolayer MoS 2 was synthesized by chemical vapour deposition (CVD) method using MoO 2 and sulfur as the precursors. The growth was completed at atmospheric pressure in an argon atmosphere. 10 mg of MoO 2 powder (99% Sigma-Aldrich) was placed in a ceramic boat. A SiO 2 /Si substrate is cleaved to the appropriate dimensions such that when it is placed in the boat (SiO 2 face down) it is at 1 mm distance from the MoO 2 powder. The boat with the MoO 2 and the silicon piece are placed in the centre of the furnace (Lindberg/Blue M with 1 in. quartz tube). 20 mg of sulfur is placed in a separate boat, upstream of the MoO 2 , at the edge of the furnace. The tube is flushed three times with argon carrier gas at room temperature before starting the growth. Whilst flowing 200 sccm of argon, the temperature is first increased to 400 • C without ramping limitations, and then the temperature is increased to 712 • C at 5 • C min −1 . Finally, the temperature is held at 712 • C for 15 min before cooling down to ambient conditions. Heavy Ion Irradiation: A 2.5 MV Pelletron ion accelerator (NEC model 7.5SH) was used to administer 225 keV Xe + ions to the MoS 2 specimens using average beam currents of either~35 or~245 nA. MoS 2 crystals deposited onto SiO 2 /Si wafers were mounted on a target station, described in detail elsewhere [63], fitted with a secondary electron suppressor operated at~400 Volts to enable accurate dosimetry. A cooling loop was used for all irradiations and the isotopically mixed beam was rastered over a 2.25 cm 2 or 4 cm 2 area, defined by either beamline slits or an electrically insulated tantalum window, respectively. Beam currents and accumulated charge were measured by a picoammeter and data acquisition was controlled using in-house built software. Key dosimetric parameters for all irradiations are available in table S1.
Raman and Photoluminescence Spectroscopy: Individual Raman spectra were measured in backscattering geometry under ambient conditions using a Renishaw InVia spectrometer equipped with a 100x objective lens and 2400 lines mm −1 grating, resulting in a resolution of~1 cm −1 . A 514.5 nm laser was operated at 0.15 mW power for all measurements to avoid thermal damage or local heating effects. Resonant Raman spectra of bulk MoS 2 crystals were measured using a 632.8 nm excitation wavelength and 1800 lines mm −1 grating (data is available in figure S30, supplementary information). Non-irradiated spectra were fitted using Lorentzian line shapes whilst irradiated spectra were deconvoluted using Voigt line shapes due to the contribution of multiple defectactivated peaks in the A 1g , E 1 2g and LA (M) spectral regions [41]. Raman maps were obtained in backscattering geometry using a WITec Alpha300 spectrometer equipped with 1800 lines mm −1 grating and a 100x objective lens, resulting in a theoretical spatial resolution of~330 nm when using a 514.5 nm laser. A grating of 600 lines mm −1 was used to increase the spectral range when obtaining photoluminescence maps. Lorentzian and Gaussian line shapes were used to create Raman and PL maps, respectively.
Atomic Force Microscopy: To complement layer number identification provided by Raman spectroscopy, micrographs and height profiles were obtained using atomic force microscopy. A Bruker Multimode 8 AFM was operated in PeakForce QNM mode using a scan rate of 0.5 Hz and a silicon nitride SCANASYST-AIR tip.
X-ray Photoelectron Spectroscopy: An Axis Ultra Hybrid (Kratos Analytical) equipped with an Al Kα x-ray source using 10 nA emission and operating at 15 kV bias was used to obtain XPS spectra. High resolution scans used a pass energy of 20 eV whilst survey scans used 80 eV. The typical vacuum level for measurements was between 2 × 10 -8 and 3 × 10 -8 mbar. The spectra were fitted with weighted Gaussian-Lorentzian line shapes using the CASAxps software and are calibrated to the C 1s signal of adventitious carbon at 284.8 eV.
Monte Carlo Calculations: Full damage cascade simulations of ten thousand 225 keV Xe + ions at normal incidence in bulk MoS 2 were carried out using the Stopping Range/TRansport of Ions in Matter (SRIM/TRIM) software [61]. Two sets of Mo and S threshold displacement energies, reported by Ghorbani-Asl et al [22] and Komsa et al [28], were used to calculate the average number of vacancies ion −1 Å −1 as a function of depth. In order to calculate the total number of displacements per atom (dpa) as a function of depth, the vacancies ion −1 Å −1 curve was multiplied by the fluence (ions Å −2 ) before being divided by the atomic density of MoS 2 (atoms Å −3 ). Finally, in order to calculate the dpa as a function of fluence, the vacancies ion −1 Å −1 verses depth curves were integrated between either 0 nm-8.4 nm (penetration depth of 514.5 nm laser in MoS 2 ) or 0-150 nm (the depth after which the concentration of vacancies becomes negligible) to obtain the average number of vacancies ion −1 . Following this, dpa was calculated by using the appropriate fluence (ions nm −2 ) and atomic density (atoms nm −3 ).
Transmission Electron Microscopy and Electron Diffraction: MoS 2 crystals produced by MME were transferred onto 400 mesh Cu Quantifoil grids (Agar Scientific) using an isopropanol-mediated transfer procedure described in detail elsewhere [64]. Bright field imaging and electron diffraction were carried out using an FEI Tecnai G2 20 S-TWIN Analytical TEM operating at 200 kV (dark field images were also obtained and are available in figure S32, supplementary information).
Thermal Annealing: A MoS 2 sample prepared by MME and irradiated to a fluence of 1.19 × 10 14 ions cm −2 was used to investigate the effects of thermal annealing on the vibrational properties of defective 1L and bulk MoS 2 crystals. The sample was placed into a tube furnace (Carbolite GHC 12/450) which had been flushed with nitrogen gas. The sample was then heated to a temperature of 150, 300 or 450 • C using a ramp rate of 5 • C min −1 under constant nitrogen flow and held at the desired temperature for 1 h before cooling to ambient temperature. After each temperature, the sample was removed and Raman spectra were acquired (thermal annealing data is available in figure  S33, supplementary information).