Tunable strain and bandgap in subcritical-sized MoS2 nanobubbles

Nanobubbles naturally formed at the interface between 2D materials and their substrate are known to act as exciton recombination centers because of the reduced bandgap due to local strain, which, in turn, scales with the aspect ratio of the bubbles. The common understanding suggests that the aspect ratio is a universal constant independent of the bubble size. Here, by combining scanning tunneling microscopy and molecular dynamics, we show that the universal aspect ratio breaks down in MoS 2 nanobubbles below a critical radius ( ≈ 10 nm), where the aspect ratio increases with increasing size. Accordingly, additional atomic-level analyses indicate that the strain increases from 3% to 6% in the sub-critical size range. Using scanning tunneling spectroscopy, we demonstrate that the bandgap decreases as a function of the size. Thus, tunable quantum emitters can be obtained in 2D semiconductors by controlling the radius of the nanobubbles.


Introduction
The advent of semiconducting honeycomb-structured transition-metal dichalcogenide (TMDC) membranes has stimulated significant research activity in the two-dimensional (2D) materials community [1].Owing to their atomic-scale thickness and their beneficial electronic and optical properties, such as their direct band gap [2] and strong spin-orbit coupling [3], TMDCs are a promising alternative to graphene.Additionally, observations of single-photon emitters in various TMDCs suggest potential quantum applications [4].This makes them appealing for fundamental studies as well as for electronics, spintronics, and optoelectronics applications (see recent reviews [1,4] and citations within).
In the TMDC family, molybdenum disulfide is the most studied material because it can be synthesized in single layers, it is relatively stable in air, and its band gap (1.86 eV) is larger than that of silicon (1.12 eV) [2], which is advantageous for low-power transistors [5].Additionally, because of its 2D nature, its electronic properties can be easily tailored by strain engineering.[6,7] Moderate biaxial tensile strains (∼2%) can decrease the bandgap size and trigger a direct-to-indirect bandgap transition [8][9][10].At larger levels of strain (∼10%), first-principle calculations predict the closure of the bandgap, thus leading to a semiconductor-to-metal transition [9,11].
For designing new devices with bespoke properties, the creation of van der Waals (vdW) heterosystems has recently attracted considerable research interest because of their potential application in different fields [12][13][14].The formation of a vertical interface can modify the electronic and magnetic properties of a membrane and induce exotic physical phenomena such as superconductivity in twisted bilayer graphene (tBLG) [15], the anomalous quantum Hall effect in tBLG aligned to hexagonal boron nitride (hBN) [16], or Mott-like states in MoSe 2 /hBN/MoSe 2 [17].At the same time, obtaining 2D materials with large areas is essential for the assembly of real devices.In this context, the assembly of MoS 2 exfoliated on Au substrates is a model system.The relatively strong interaction of MoS 2 with Au leads to the exfoliation of high-quality monolayer membranes that exhibit a large area [18].On (111)-oriented substrates, a moiré pattern develops at the MoS 2 -Au interface, modulating the electronic properties of the surface through local variations in the MoS 2 -Au bonds [19][20][21].In particular, the hybridization of the S 3p z orbitals with the Au d bands induces distortions in the valence band of MoS 2 [20] and bandgap renormalization [21].
The adventitious contamination present on individual layers before the assembly of TMDC/substrate systems leads to the formation of stable bubbles with micro-and nanometer sizes [22].The vdW forces that bind the adjacent layers facilitate the coalescence of trapped matter in minute bubbles, thus creating large inter-bubble areas that are atomically sharp and free of contamination via a process known as "selfcleansing."[23] Recent investigations have shown that the presence of bubbles leads to variations in the local properties of 2D materials [24,25].In addition, nanobubbles formed in different 2D heterostructures can serve as quantum-dot emitters because of the interplay between the strain-induced confinement [25] and changes in the local dielectric environment [26].Notably, the aspect ratio between the height and radius of commonly occurring bubbles exhibits a universal scaling, which results from the balance between the membrane-substrate adhesion and the elastic properties of the membrane [22,[27][28][29].Detailed knowledge of the strain levels emerging in these bubbles is fundamental to understanding the local, strain-tunable electronic properties of monolayer TMDCs and related 2D materials.
In this study, we assess the shape and strain of intrinsic nanobubbles formed in MoS 2 monolayers exfoliated on Au(111).Room-temperature scanning tunneling microscopy (STM) is employed to obtain atomically resolved images of MoS 2 nanobubbles with radii of < 25 nm.Our STM measurements show that, below a critical radius (≈10 nm), the MoS 2 bubbles exhibit decreasing aspect ratios with decreasing size.We compare our experimental results with molecular dynamics (MD) simulations, which support the aspect ratio scaling breakdown observed in our STM experiments.Our MD results also reveal the presence of ultra-high pressures (up to 10 GPa) inside the smaller nanobubbles.Larger MoS 2 bubbles, conversely, exhibit pressures of the order of 100 MPa.The radial pair distribution function (RPDF) of the matter trapped inside the MD bubbles suggests a transition from a gas-like phase in larger nanobubbles to a liquid-like phase inside smaller bubbles.Additional atomic-level (STM and MD) analyses show that the strain follows the evolution of the aspect ratio, increasing with increasing bubble size from a minimum of ≈3% for the smallest bubbles.Above a critical radius of ≈10 nm, the strain stabilizes around a universal value as high as 8%.Strain maps also show the evolution of the local strain in the bubbles.The tunability of the strain is additionally supported by scanning tunneling spectroscopy (STS) measurements of the bandgap in MoS 2 nanobubbles.Our measurements indicate a bandgap decrease rate of approximately -140 meV/% strain in bubbles with radii of < 10 nm.

STM measurements of MoS 2 nanobubbles
The nanobubbles are analyzed using STM topographic images of MoS 2 monolayers exfoliated using the gold-assisted method [18] on Au(111) films grown on mica (see Methods).The samples are prepared in air and then transferred to ultra-high vacuum (UHV).We find that annealing at 400 K is sufficient to desorb most contaminants and obtain high-resolution STM images, as shown in Fig. 1.
Figure 1a shows the topographic scan of a 70 nm × 70 nm area displaying a continuous moiré pattern that is indicative of MoS 2 -covered Au terraces [30,31].The nanobubbles are imaged as bright protrusions, which are predominantly found at the step edges of the Au substrate.This is consistent with the generally accepted notion that the bubbles originate from molecules trapped between the 2D membrane and the substrate, as the step edges are the energetically favored adsorption sites.The coverage of the protrusions on the surface ranges between 2% and 6%.In magnified STM scans, as shown in Fig. 1b, the atomic structure of MoS 2 nanobubbles is resolved.We observe a hexagonal lattice that is consistent with that of the top S atoms in 2H-MoS 2 [32,33].Although the moiré pattern is visible on the adjacent terrace (Fig. 1a), no superstructure is observed on the bubbles by STM (Fig. 1b).This confirms the decoupling of MoS 2 from the Au substrate over the nanobubble area.The interatomic distances, which measure 3.17 ± 0.03 Å on the moiré-characterized regions, range between 2.9 Å and 3.5 Å on the surface of the bubbles.The variation in these distances is due to local variations in strain, as we discuss below.
The bubbles exhibit elliptical/lenticular shapes and bell-like profiles (Fig. 1c).We measure their maximum height, h, and radius, R, along their longest axis.To ensure consistent measurements for nanobubbles formed at step edges and on terraces, we consider all bubbles as only being formed on the highest adjacent terrace.That is, for bubbles formed at step edges, h and R are measured by projecting the level of the highest terrace to the other side of the bubble (see Fig. 1c).Then, we obtain the aspect ratio, h/R, which is shown as a function of R for more than 170 bubbles in Fig. 1d.For R > 10 nm, the aspect ratio measured by STM exhibits a plateau at h/R = 0.21±0.01.The plateau is consistent with the universal scaling of bubbles demonstrated using different 2D materials, where the h/R value is determined solely by the adhesion energy and mechanical properties of the membrane [22].We confirm this h/R value for large bubbles by atomic force microscopy (AFM) at ambient pressure (h/R = 0.22 ± 0.02), thus excluding any influence of the external pressure or apparent-height effects in the STM measurements (see Supplementary section S1).
In bubbles with R ≤ 10 nm, we find h/R values that decrease with decreasing R. The lowest h/R value found in our experiments is 0.07 for the smallest bubbles (R ≈ 2 nm).Thus, the universal scaling of the aspect ratio appears to break down for small nanobubbles with sizes well below those investigated by Khestanova et al. [22].We define a critical radius, R c , below which the h/R breakdown occurs, and estimate it to be ≈ 10 nm from STM measurements.
We fit the profiles of the MoS 2 bubbles to the general shape proposed by Blundo et al. [27], where the height of the bubble, z, is defined as a function of the radial distance from the center, r * , as z(r . Fitting the profiles from different experimental and simulated large-sized bubbles (R > 100 nm) formed in monolayer membranes, Blundo et al. constructed a relation between h/R and the exponent q, whereby an exponent of q ≈ 2.3 is expected for h/R = 0.22.In our MoS 2 bubbles, we find a value of q ≈ 2.3 for R > R c .This agreement implies the absence of residual strain and defects in our membrane, as we would expect them to modify the shape of the nanobubbles.For smaller bubbles (R < R c ), we obtain lower values of q (as low as 1.83, see Supplementary section S2), which is consistent with the observed decrease in the h/R ratio in this regime.

MD simulations of MoS 2 nanobubbles
We perform all-atom MD simulations to study the formation of MoS 2 nanobubbles by employing the LAMMPS code [34].Our computational setup is designed to capture the creation and relaxation of MoS 2 bubbles, in which a monolayer MoS 2 membrane is accommodated on top of a flat (111)-oriented Au substrate whose surface contains adhered contamination-modeled in our simulations as nitrogen atoms.During the relaxation, the contaminants trapped between the MoS 2 and the substrate (see Supplementary Fig. S3) form a MoS 2 bubble at the center of the periodic MoS 2 /Au(111) system (see the Supplementary video).For further details, see Methods and Supplementary section S3.
We also simulate the formation of MoS 2 bubbles at Au(111) steps of monoatomic and biatomic height.Supplementary section S4 provides a comparative analysis of relaxed MD bubbles with R ≈ 17 nm (where R > R c ) formed at Au step edges and on flat Au(111).According to our simulations, the resulting h/R is 0.20 and 0.22 when the bubble is formed at a biatomic step and on a flat surface, respectively.Thus, our MD results indicate that the maximum error in our experimental h/R measurements is ≈ 10% when comparing bubbles formed at steps and on a flat surface.
Figure 2a shows the profiles of relaxed MD bubbles for different R. Further fitting of the bubble profiles yields values of the exponent q that range from approximately 1.6 to 2.2 as h/R increases from approximately 0.12 to 0.21 (see Supplementary Table S1).In addition, our MD simulations clarify the role of the vdW-type adhesion between the membrane and the substrate in determining the resulting aspect ratio of relaxed bubbles.As expected, we find that increasing the adhesion between the MoS 2 and Au substrate, which was modeled using Lennard-Jones (LJ) potentials (see Supplementary sections S5 and S6), gradually increases the aspect ratio of MoS 2 nanobubbles.Additional results and discussions are provided in Supplementary section S7.
Following the h/R evolution shown in Fig. 2b, our atomistic simulations also support the scaling breakdown revealed by the STM measurements for R < R c .Our MD simulations indicate a slightly larger R c (≈ 13 nm) than that observed in our STM analysis (R c ≈ 10 nm) (see Figs. 1d and 2b).
We map the evolution of the hydrostatic pressure of the trapped atoms, P , as a function of the MoS 2 nanobubble size.Supplementary section S8 describes the approach we adopt to extract P from our MD simulations.The resulting evolution of P as a function of R shown in Fig. 2c indicates ultra-high pressures of up to 10 GPa inside MoS 2 bubbles with R < 10 nm.This is consistent with previous reports on graphene nanobubbles of comparable size [35,36].P levels are reduced toward pressures on the order of 100 MPa when R > 10 nm (see the inset to Fig. 2c).A similar P (R) behavior was confirmed in recent experimental and analytical studies of nanobubbles formed in diverse vdW heterosystems, where the pressure inside the nanobubbles is inversely proportional to h [22,37].
Figure 2d displays the RPDFs, g(r), of the trapped atoms as a function of R. In general terms, the RPDF represents the probability of finding a particle at a distance r from a particle placed at position r=0.The g(r) shape indicates the phase of the trapped material inside the MD nanobubbles.In this framework, the RPDFs plotted in Fig. 2d suggest the presence of an "ordered" liquid-like phase inside nanobubbles of small size, as evidenced by the second and third peaks [38] of the g(r) − r curves for R < R c .The transition to the liquid phase in small nanobubbles is associated with the accumulation of extreme pressures (P > 2 GPa) inside the bubbles.Recent density functional theory (DFT) results for graphene nanobubbles support the evolution of the RPDFs with increasing P suggested by our MD results [37].The RPDFs for R > R c , where nanobubbles exhibit much lower P levels, tend to follow those usually found in disordered gas-like fluids, as indicated by the softening of the secondary peaks.These findings agree with the continuum analysis reported by Iakovlev et al. [39], which indicates that gas-filled graphene nanobubbles develop a higher aspect ratio than liquid-containing ones.

Strain in MoS 2 nanobubbles
We calculate the local atomic-level strains, ε local , across our relaxed MD nanobubbles by computing the individual atomistic displacements during the bubble formation and relaxation (Supplementary section S9).Figures 3a and 3b display the resulting ε local maps from MoS 2 membranes containing small (R < R c ) and large (R > R c ) nanobubbles, respectively.Contrary to the common notion that the maximum strain in 2D bubbles is located around the top [22,24], our maps display a doughnut-like strain pattern in MoS 2 bubbles, where the ε local global maxima are detected close to the bubble periphery and the local minima are detected at the center (see Fig. 3c).These strain patterns (Fig. 3a,b) coincide with those recently reported for MoS 2 [26], WS 2 [29], and WSe 2 [25,40] bubbles.Our MD analysis also reveals that the development of such a strain pattern leads to a complex deformation state consisting of uniaxial, biaxial, and shear components (Supplementary section S10).
Figure 3d shows the average strain, ε bubble = ⟨ε local ⟩, obtained inside the simulated MoS 2 bubbles as a function of the radius, where ⟨ε local ⟩ denotes the average value of the ε local data in the bubble.The average (in-bubble) strain gradually increases from ≈ 3% to ≈ 8% when R increases from 2.5 nm to R c .STM measurements of the atomic-scale strains occurring in our experimental MoS 2 nanobubbles confirm the strain scaling found in our simulated bubbles with R < R c .For R > R c , we obtain an average strain of ≈ 8%, which is found to be independent of R.
In connection with this strain scaling, our MD results reveal that the ε local distributions follow different shapes in the R < R c regime from those observed in the (universal) large-bubble regime (R > R c ) (see Figs. 3e,f).In this context, small MoS 2 nanobubbles tend to develop broad ε local distributions, whereas the local strains of large bubbles adhere to right-skewed distributions with a sharp peak at ε local ≈ ⟨ε local ⟩.Overall, these results evidence that the in-bubble strain scales with R in a similar manner to the aspect ratio (compare Figs. 3d and 1d).

Bandgap of MoS 2 nanobubbles
We acquire spectra near the center of 11 bubbles with different sizes and measure their local bandgap by STS.Because strain can be used to tune the bandgap size, the size of the bandgap can in turn be used to measure the strain.The dI/dV spectra are shown in Fig. 4a.Here, the bandgap is identified as a region of zero conductance.The presence of this zero-signal region on top of bubbles confirms the decoupling of MoS 2 from the substrate, on which the strong MoS 2 -Au hybridization normally manifests as a non-zero signal inside the gap and leads to bandgap renormalization [21] (see Supplementary section S12).Because the states far from the Γ point (such as those at the band edge of unstrained MoS 2 ) are difficult to detect by STS [20], we use the onsets of the tunneling conductance from its zero level on the two sides of the Fermi level as a measure of the bandgap edges.The measured bandgap corresponds to the quasiparticle bandgap, rather than the optical bandgap which is lowered by the exciton binding energy [41][42][43].The energies of the conduction band minimum (CBM) and valence band maximum (VBM) are plotted as a function of R in Fig. 4b.The VBM upshifts in energy from approximately -1.6 to -1.4 eV with increasing R, whereas the CBM downshifts from approximately 0.6 to 0.4 eV.As a result, the local bandgap of our nanobubbles decreases from E g = 2.2 eV for bubbles with R = 2.4 nm to ≈ 1.8 eV for those with R = 7.9 nm, as shown in Fig. 4c.The data are fitted to a linear regression, displayed as a visual guide in the figure .A linear approximation is consistent with the evolution of ε bubble with increasing R for R < 10 nm (see Fig. 3f), whereas a plateau is expected for larger R.
We acknowledge that the large error in our room temperature STS measurements is due to thermal broadening (∆E = (3.3k B T ) 2 + (2.5 V mod ) 2 = 0.125 eV, where the thermal broadening of 3.3 k B T accounts for approximately 0.09 eV at T = 300 K, and V mod = 30 mV is the modulation applied to the bias).Two measurements on the same bubble of radius R = 2.4 nm (Fig. 4) indicate the magnitude of our error.

Discussion
In large bubbles formed in 2D materials, a constant aspect ratio is determined solely by the balance between the substrate-2D membrane adhesion and mechanical properties of the membrane [22].However, in the STM experiments and MD simulations, we find that this universal scaling is broken for MoS 2 nanobubbles with R < R c , which exhibit lower aspect ratios.At a first glance, such behavior may be regarded as a nanoscale effect caused by the increasing edge-to-bulk ratio (the 2D equivalent of the surfaceto-volume ratio) as the bubble size decreases.Here, the relevant edge and bulk effects would be the membrane-substrate adhesion and the bending rigidity of the membrane, respectively.However, the adhesion tends to minimize R and the bending rigidity limits the membrane curvature and h.Hence, in this simplified scenario, greater h/R values would be expected when the bubble size decreases.This is opposite to our observations for R < R c .
Turning to membrane theory, the aspect ratio of the bubbles is driven by the minimization of the total energy, Here, E bend is the out-of-plane bending energy, E el is the in-plane elastic energy of the membrane, E b is the free energy of the trapped substance, and E vdW is the energy required to separate the membrane and substrate.The universal aspect ratio (expressed as h/R = (πγ/5c 1 Y ) 1/4 , where γ is the total adhesion, Y is Young's modulus, and c 1 is a constant) is valid above a critical size, where the bending rigidity of the membrane can be neglected [22].At smaller sizes, however, the minimization of the total energy is determined by the bending rigidity term.As a result, for small radii, the h(R) relation shifts from a linear dependency to a higher-order one [22,35].
Along these lines, a model derived by Lyublinskaya et al. [44] predicts an Rdependent aspect ratio evolution in small bubbles.Using non-linear plate theory and membrane theory, the authors showed that the energy minimization is dominated by the bending term below a critical radius that depends on temperature.This leads to a h/R ∼ R dependence followed by a h/R ∼ R 1−η/2 one (with η = 0.8 for a 2D membrane) in the small-size limit.For radii above R c , the bending energy is renormalized and the energy minimization is dominated by the elastic energy.In this range, the h/R ratio is determined by the pressure-induced tension (tension-dominated regime), which leads to the well-established universal scaling [44].Hence, in the context of membrane theory, the scaling breakdown observed in our MoS 2 nanobubbles is explained by the onset of the bending-dominated regime for R < R c .Notably, by combining STM experiments and MD simulations, Villarreal et al. [35] reported an opposite behavior at very low radii, observing increasing h/R ratios as the bubble size decreases.Their nanobubbles (0.5 nm < R < 3 nm) are formed by sputtering noble-gas atoms through a graphene layer grown on Pt(111).However, in this system, a minimum bubble height, corresponding to a monolayer of trapped gas, is obtained for R > R c (R c ≈ 0.4 nm for graphene).The aspect ratio enters a h/R ∝ 1/R regime rather than the bending-dominated one, and therefore the breakdown observed in this case has a different origin.
Our STM measurements and MD simulations also demonstrate an increase in strain with increasing sizes of nanobubbles for R < R c .Although the computational costs and the constraints on STM resolution (Supplementary section S11) limit the number of data points at larger radii, the average strain appears to reach a plateau centered at ε bubble = 8.2 ± 0.9%.This is consistent with the strain scaling as ε ∝ h 2 /R 2 in bubbles formed in 2D materials [22] (see Figs. 1d and 2b).In a similar MoS 2 /Au(111) system, Pető et al. [8] found the best fit for their Raman spectra and the bandgap measured by STS for a biaxial strain of 2% applied to flat MoS 2 .However, our MD analyses of strain (Supplementary Fig. S9) indicate that MoS 2 nanobubbles exhibit a complex distribution of deformation states that differ from the purely biaxial scenario proposed by Pető et al. [8].Thus, the combination of firstprinciples simulations and Raman measurements may lead to a misestimation of the average strain in nanobubbles.Recent studies also reported lower characteristic values of strain (ε bubble ≈ 0.6 − 2%) in large MoS 2 bubbles, obtained by photoluminescence (PL) and Raman spectroscopy [24] or mechanical models [22,45].Similar values were also found in WS 2 [29,46].This discrepancy might be due to the PL/Raman spot size that is often close to one micrometer, which is larger than the typical size of bubbles formed in 2D materials.Thus, averaging with the flat, unstrained membrane must be considered.This idea is supported by our complementary simulations of the straininduced Raman E mode shift distribution, which compare the spectra of flat MoS 2 areas that contain bubbles with bubble-only areas (see Supplementary section S13).Furthermore, using STM topographies of nanobubbles, much higher strains (≈ 10%) have been reported [35].We note that the breaking point of pristine MoS 2 layers is ≈ 11% [47], although lower values have been reported [48].
Combining the STS data with our strain analyses, we estimate a decrease rate in the bandgap as a function of bubble size of approximately -70 meV/nm for R < R c (2 nm≤ R ≤ 8nm).The strain increases linearly by ≈ 3% in the same R range (see Fig. 3f), which yields a bandgap decrease rate of approximately -140 meV/% strain for R < R c .This value agrees well with the -124 meV/% recently obtained by combining PL measurements and DFT calculations [7] and is in overall agreement with other reported bandgap decrease rates under biaxial strain (ranging between -45 and -99 meV/% for the optical bandgap measured by PL [24,48,49] and between -190 and -230 meV/% for the indirect bandgap obtained from DFT calculations [8,50]).For MoS 2 under uniaxial strain, rates as high as -136 meV/% strain were observed by combining PL/Raman measurements and DFT calculations [51].In the regime of linearly increasing strain obtained from the STM measurements and MD simulations for R < R c , our STS analysis confirms a continuous tuning of the bandgap with increasing bubble size.We expect the bandgap size to reach a constant minimum for R > R c , following the evolution of the strain with an inverse slope.Thus, the results in this work show that, in the subcritical size limit, the strain and bandgap of the nanobubbles can be tailored by controlling their radii.

Methods Experimental
MoS 2 was exfoliated on Au(111) using a preparation method that predominantly yields monolayer flakes with sizes of at least hundreds of micrometers [18].The samples were prepared at ambient pressure by flame-annealing the gold substrates (300 nm epitaxial Au(111) grown on mica, Georg Albert PVD, Heidelberg, Germany) and exfoliating a freshly cleaved MoS 2 crystal on top of it, in rapid sequence.The samples were analyzed using an optical microscope to ensure sufficient MoS 2 coverage before being placed into a UHV chamber (with a pressure < 1 × 10 −10 mbar), minimizing their exposure to ambient air to just a few minutes.To remove contamination from the surface while preventing structural modifications or the creation of defects [52], the MoS 2 /Au(111) samples were annealed at a low temperature (400 K) for one hour before each STM measurement.STM was performed at room temperature using a variable-temperature SPECS SPM Aarhus 150 equipped with a tungsten tip.To ensure consistent measurements of R and h, we only considered images with lateral sizes of 100-150 nm acquired at bias voltages of 1-1.5 V.The former requirement is necessary to avoid different non-linear background effects due to long-range scanning, and the latter is required to exclude offsets due to bias-dependent apparent heights (Supplementary section S1).The post-processing and analysis of the STM images were carried out using Gwyddion software (http://gwyddion.net)[53].

Computational
We constructed MD cells containing the MoS 2 /Au(111) heterosystem with characteristic areas ranging from 10 nm × 10 nm to 140 nm × 140 nm (see Supplementary Table S2).The as-built MoS 2 membranes have a bell-like shape that allows the introduction of trapped atoms between the MoS 2 and the Au substrate (see Supplementary Figs.S3a and b).Periodic boundary conditions are imposed on the two sides of the computational domains, whereas the top and bottom walls are non-periodic.The fixed vertical dimension of the MD cells is set to h 0 +20 Å, where h 0 is the height of the as-built MoS 2 bell.Additional details of MD cell construction are provided in Supplementary section S3.
To describe the intralayer atomic interactions in the MoS 2 membranes, we use the reactive atomic bond order potential parameterized by Liang et al. [54], which generates consistent results with those from DFT calculations.The Au substrate is composed of two rigid, (111)-oriented atomic layers with a face-centered cubic lattice parameter of a = 4.08 Å [55].The atomic interactions in the substrate are excluded from the calculations so that the positions of Au atoms remain fixed during the MD simulations.We choose N atoms to simulate the trapped contamination between the MoS 2 membrane and the Au substrate.The interatomic forces inside the bubbles are modeled using the 12-6 LJ potential (see Supplementary Eq. (S3) and Table S4).We also employ this model to describe the interactions between the trapped atoms with the MoS 2 membrane and the Au substrate, where the cross LJ parameters are obtained using the Waldman-Hagler mixing rules; see Supplementary section S5.The MoS 2 membrane and the Au(111) substrate mutually adhere through weak vdW forces that are also described by the LJ potential.The employed LJ values lead to ultra-high adhesion energies between the MoS 2 and Au(111) (see the discussion in Supplementary section S7).The LJ cutoff distance is set to 9 Å.
The energy of the as-built MoS 2 /N/Au(111) system (Supplementary Fig. S3a) is initially minimized using the conjugate gradient method.Then, the system is relaxed over a 500-ps run.The dynamics of the trapped atoms follow the microcanonical NVE ensemble whereas the atoms of the MoS 2 membrane follow the canonical NVT ensemble with the Nosé-Hoover thermostat maintaining the temperature of MoS 2 at approximately 300 K.The computational timestep is set to 0.5 fs.Atomic visualization and data analysis of the MD results are conducted using the OVITO package (www.ovito.org)[56].

Fig. 1
Fig. 1 STM measurements of MoS 2 /Au(111) nanobubbles.(a) STM topography of monolayer MoS 2 exfoliated on Au(111) displaying the characteristic moiré-patterned terraces separated by monoatomic steps.The bright protrusions correspond to nanometer-sized MoS 2 bubbles.Setpoint for the measurement: bias voltage V bias = 1.5 V, tunneling current It = 0.3 nA.(b) STM image of a MoS 2 nanobubble of R = 1.9 nm.Setpoint: V bias = 0.5 V, It = 0.6 nA.(c) Profile of a MoS 2 nanobubble (R = 4.2 nm), marked with a green line in (a).(d) The evolution of h/R with R in MoS 2 nanobubbles as obtained by STM and AFM measurements.The gray line obtained by exponential regression provides a visual guide.Each data point corresponds to a unique MoS 2 bubble.

Fig. 2
Fig. 2 Properties of MoS 2 nanobubbles obtained from MD simulations.(a) Bubble profiles with varying sizes, R. The profile for R = 21.0 nm is fitted by adopting the model proposed by Blundo et al. [27].(b) Evolution of the h/R ratio as a function of R. (c) Calculated hydrostatic pressure, P , inside relaxed MoS 2 nanobubbles of different sizes (Supplementary section S3).The gray line obtained by exponential regression provides a visual guide.Each data point in (b, c) accounts for one MD bubble.(d) RPDFs of the atoms trapped between the MoS 2 and the Au(111) substrate for different bubble sizes.The color scheme follows that given in (a).

Fig. 3
Fig. 3 Strain in MoS 2 nanobubbles.ε local maps of MoS 2 from (a) 25 nm × 25 nm and (b) 130 nm × 130 nm periodic MD cells containing a relaxed bubble of R = 6.7 nm and R = 26.4nm, respectively.For visual guidance, the bubble's perimeters are drawn in (a, b) with a dashed yellow contour.(c) Strain and height profiles-marked with a green line in (a)-across the MoS 2 nanobubble.(d) ε bubble vs. R obtained from STM and MD.The strain estimation procedure in STM MoS 2 nanobubbles is described in Supplementary section S11.The MD ε bubble calculations correspond to the averaged (in-bubble) ε local values obtained inside the MoS 2 bubbles.ε local distributions of the MoS 2 bubbles with R > Rc and with R < Rc(≈ 13 nm) are given in (e, f), respectively.

Fig. 4
Fig. 4 Bandgap measurement by STS.(a) dI/dV spectra acquired near the center of 11 nanobubbles of increasing radii.The inset is a representative STM image of a bubble (R = 9.9 nm), and the location of the tip for the measurements is marked with a cross.Two spectra with different tip apexes are shown for a bubble of R = 2.4 nm.The setpoint for the measurements was V bias = 1.0-1.2V, It = 0.5 nA.(b) The energy of the valence and conduction band onsets extracted from the dI/dV spectra, as shown in (a) for the bubble of R = 2 nm.(c) The calculated bandgap (Eg = E CBM −E VBM ) as a function of R. Notice that, for R > 10 nm, stabilization at a constant value is expected.