When layers of van der Waals materials are deposited via exfoliation or viscoelastic stamping, nanobubbles are sometimes created from aggregated trapped fluids. Though they can be considered a nuisance, nanobubbles have attracted scientific interest in their own right owing to their ability to generate large in-plane strain gradients that lead to rich optoelectronic phenomena, especially in the semiconducting transition metal dichalcogenides. Determination of the strain within the nanobubbles, which is crucial to understanding these effects, can be approximated using elasticity theory. However, the Föppl–von Kármán equations that describe strain in a distorted thin plate are highly nonlinear and often necessitate assuming circular symmetry to achieve an analytical solution. Here, we present an easily implemented numerical method to solve for strain tensors of nanobubbles with arbitrary symmetry in 2D crystals. The method only requires topographic information from atomic force microscopy and the Poisson ratio of the 2D material. We verify that this method reproduces the strain for circularly symmetric nanobubbles that have known analytical solutions. Finally, we use the method to reproduce the Grüneisen parameter of the E′ mode for 1L-WS2 nanobubbles on template-stripped Au by comparing the derived strain with measured Raman shifts from tip-enhanced Raman spectroscopy, demonstrating the utility of our method for estimating localized strain in 2D crystals.

In recent years, 2D van der Waals crystals have emerged as a promising class of materials for optoelectronic applications owing to their stunning physical properties and unique tunability that is enabled by stacking different 2D crystals into vertical heterostructures. When atomically thin materials are deposited on a substrate by exfoliation, delamination of the material from the host substrate is a common occurrence, leading to the formation of structural imperfections such as nanoscale wrinkles, folds, and bubbles.1 Nanobubbles, which form when trapped materials between the atomic layer and the substrate coalesce into a local pocket of pressure, are among the most ubiquitous of these imperfections. Though often undesired for device applications, their ability to modify the electronic structure of the atomic layer has spurred significant interest in the physics of the nanobubbles themselves.1–4 In the case of the molybdenum and tungsten families of the transition metal dichalcogenides (TMDs), nanobubbles significantly alter the optical bandgap and energy of excitonic photoluminescence by up to several hundred millielectronvolts.5 The inhomogeneity of the strain has also been shown to induce complex phenomena such as exciton funneling2,6,7 and single-photon emission at cryogenic temperatures,3 motivating efforts toward developing strained 2D TMD-based devices for solar light harvesting and optical quantum information processing.

Determining the magnitude and distribution of strain within individual nanobubbles is critical to understanding their fundamental properties. This characterization is often done by combining atomic force microscopy (AFM) to measure topographical parameters—such as the maximum height and radius of a nanobubble—with models that are based on elasticity theory and treat the monolayer TMD as a thin plate that is subjected to a transverse load.1,5,8 This approach, however, presents a challenge as the governing equations that describe the behavior of a thin plate, the Föppl–von Kármán (FvK) equations, are highly nontrivial to solve owing to their coupled nonlinearity.8,9 Many studies have made use of Hencky’s solution for thin circular plates with clamped edges under a uniform transverse load.10 This solution, for example, yielded estimated values for the Grüneisen parameters of graphene that agreed well with those derived from ab initio density functional theory.11 Hencky’s solution has also been employed to estimate the van der Waals adhesion energies of both graphene12 and MoS2,13 as well as the shift of the optical bandgap of MoS2 under biaxial strain.5 

However, Hencky’s solution is ultimately limited in applicability. In nanobubbles that are formed by trapped material, the edges are not clamped but adapt according to the balance of the attractive forces of the monolayer to the substrate and the pressure of the trapped material. Nanobubbles, for instance, can even be moved relative to the top layer and substrate by means of contact AFM.14 For nanobubbles of this type, Khestanova et al. solved the FvK equations for circularly symmetric nanobubbles in MoS2, graphene, and hBN. Remarkably, the authors found a universal distribution of strain and self-consistent scaling phenomena,1,8 which indicates that an ideal nanobubble of a van der Waals material on a given substrate will have a constant shape and strain distribution, regardless of size. The magnitude of the strain was found to depend only on the square of the nanobubble aspect ratio, hmaxr, which is determined by the attractive force of the monolayer to the substrate and stiffness of the monolayer material. Thus, strain within the nanobubble can be tuned by changing the substrate, which provided means to estimate the approximate strain in circular MoS2 nanobubbles.2 

While these applications of elasticity theory to circular nanobubbles have yielded great insight into adhesion energies and the effect of strain on excitonic photoluminescence of TMDs, the restriction to nanobubbles of circular symmetry limits the use of these models on real samples. Nanobubbles in TMDs are not generally circularly symmetric and can take on arbitrary and amorphous shapes,14 which may have significantly different strain distributions compared to the relatively simple case of circular nanobubbles. Recently, advanced atomistic theory has shown that the effect of strain from nanobubbles in monolayer TMDs can be dramatically different than the results of classical elasticity theory for nanobubbles with large aspect ratios.15 According to classical plate theory, the maximum strain of a circularly symmetric nanobubble is at its center,16,17 which would then correspond to the maximum detuning of band-edge electronic states.5 By contrast, more sophisticated calculations of MoS2 nanobubbles use valence force field theory combined with tight-binding models that include effects of strain as well as changes in dielectric screening when the material separates from the substrate. These models predict that the greatest strains and shifts in energy levels are concentrated at the nanobubble edges.15 Further, the atomistic theory predicts that for nanobubbles of sufficiently high aspect ratios, the atomic lattice will form wrinkles around the nanobubble boundary. A key advantage of the more microscopic approach is that assumptions of circular symmetry are unnecessary, and arbitrary nanobubble morphologies can be considered. Using high-resolution scanning near-field optical microscopy (SNOM), we recently showed that the nanoscale spatial distributions of excitonic energies and photoluminescence intensities within WSe2 nanobubbles indeed correspond to the nanoscale distributions of strain that are predicted by the atomistic theory,4 experimentally confirming the importance of understanding nanoscale strain distributions in nonideal, highly nonsymmetric nanobubble shapes.

The atomistic approach offers a significantly more accurate description of the nanoscopic strain and its effect on the electronic band structure, but the requisite calculations are computationally expensive, limiting its use for routinely evaluating complex strain profiles in TMDs and other atomically thin van der Waals materials. Further, the universal scaling observed in Ref. 1 shows that the FvK equations do describe the average mechanical properties of nanobubbles to a reasonable degree. And estimates of the overall energetic shift of excitonic photoluminescence under biaxial strain by the mechanical blister tests in Ref. 5 agree with predictions from ab initio density functional theory.18 These early successes lead us to ask: can the FvK equations be utilized to describe localized strain effects of excitons in non-symmetric TMD nanobubbles, enabling more rapid analysis of complex geometries?

Here, to better understand the applicability of the FvK equations to TMD nanobubbles, we develop a numerical method to estimate the strain in nanobubbles of arbitrary symmetry. This method relies only on the experimentally measured AFM topography and the Poisson ratio of the material to calculate the full strain tensor for nanobubbles of arbitrary shape. We verify that this method reproduces the strain for circularly symmetric bubbles of different height profiles, where the analytical solutions are available. We then apply this method to nanobubbles formed in monolayer WS2 deposited on template-stripped gold. Enabled by the gold substrate, we use nanoscale photoluminescence and Raman measurements to show that the predicated strain tracks with the PL intensity and energy. Lastly, we show that our method can estimate nanoscale strain distributions of large areas of many isolated and partially overlapping nanobubbles simply with AFM imaging, filling a key gap in scale between atomistic modeling and pure elasticity theory for understanding and predicting nanoscopic strain distributions of nanoscale to microscale strain-engineered 2D systems. Using this approach, we rederive the Grüneisen parameter for the E′ mode of monolayer WS2, which agrees well with previous experiments.19,20 Our method is thus fully consistent with established plate theory models, and by using the measured topography as a source function, the method captures structural details such as kinks and wrinkles, which are not captured in the known analytical solutions. This procedure thus provides a facile way to map strain distributions in monolayers and vertical heterostructures of TMDs and other van der Waals layered materials.

A cross section of a hypothetical nanobubble is shown in Fig. 1(a). Trapped material forms a pocket of pressure, which is the source of a deformation and thus strain of the nanobubble. The 2D material of the nanobubbles is deflected in the vertical (z) direction by an amount that is clearly much greater than the thickness of the plate, as is commonly observed in nanobubbles in 2D materials.1 The FvK equations are a coupled set of equations that describe the stresses and deflection of a thin plate that is subject to such a displacement, where the out-of-plane deflection is large compared to the thickness of the plate.1,8,9,17,21

FIG. 1.

Comparison of the numerically calculated strain of circularly symmetric bubbles of different profiles to the exact strain as derived in Rostami et al.8 (a) Diagram of a nanobubble of monolayer WS2. (b) Diagram of the xz-cross section of the deformation induced in a thin plate by a transverse load. The variables ux and uz are the body displacements of the plate from the unstrained (i.e., normal) geometry. (c) Normalized height profiles of the three nanobubble cases: Gaussian (blue), parabolic-like (red), and pure bending (yellow-orange). (d) The corresponding calculated strain values from each case of nanobubble geometry in (c) as a function of radial distance. The solid line is the exact radial strain as derived in Rostami et al.,8 and the symbols are the numerical values calculated by solving for the Airy stress function in Eqs. (7) and (9). For the complete formula of the height profile and the strain, see Appendix D of Rostami et al.8 

FIG. 1.

Comparison of the numerically calculated strain of circularly symmetric bubbles of different profiles to the exact strain as derived in Rostami et al.8 (a) Diagram of a nanobubble of monolayer WS2. (b) Diagram of the xz-cross section of the deformation induced in a thin plate by a transverse load. The variables ux and uz are the body displacements of the plate from the unstrained (i.e., normal) geometry. (c) Normalized height profiles of the three nanobubble cases: Gaussian (blue), parabolic-like (red), and pure bending (yellow-orange). (d) The corresponding calculated strain values from each case of nanobubble geometry in (c) as a function of radial distance. The solid line is the exact radial strain as derived in Rostami et al.,8 and the symbols are the numerical values calculated by solving for the Airy stress function in Eqs. (7) and (9). For the complete formula of the height profile and the strain, see Appendix D of Rostami et al.8 

Close modal

The question that we aim to answer here is: given a known out-of-plane deflection and Poisson ratio of a thin plate, what is the corresponding in-plane strain induced in the plate? In the following, we summarize the relevant aspects of the theory to numerically solving this problem. More details are provided in the supplementary material, and a full treatment can be found in the many classic texts on elasticity theory such as Refs. 9 and 21. The key relationship is the dependence of the stress on the displacement. This relationship is described with two sets of equations: one relating the stresses to the strains and the other relating the strains to the displacements. The stress–strain relationships are described by Hooke’s law,

σxx=E1υ2εxx+υεyyσxy=E1+υεxy
(1)
σyy=E1υ2εyy+υεxx.

Here, σij and εij are the components of the stress and strain tensors, respectively, E is Young’s modulus, and υ is the Poisson ratio of the material. The strain–displacement equations are as follows:

εxx=uxx+12hx2,εxy=12uxy+uyx+hxhy,
(2)

and

εyy=uyy+12hy2,

where ux and uy are in-plane body displacements, and h = uz is the resulting height (the displacement in z) as summarized in Fig. 1(b). The condition for equilibrium, which is derived in Ref. 9 and summarized in the supplementary material, is that the divergence of the stress tensor is zero, yielding two coupled equations in terms of the body displacements:

σxxx+σxyy=0
(3a)

and

σxyx+σyyy=0.
(3b)

For the general problem, a third equation is needed to relate in-plane stresses to the applied transverse load (in this case, the pressure of the trapped material inside the nanobubble) to completely solve the system. But for the problem of calculating the strain of nanobubbles with known deflections (i.e., known h), these equations form a complete set for the unknown stresses (or displacements). In the case of nanobubbles with circular symmetry, the problem reduces further to a single second-order ordinary differential equation that has been solved analytically for nanobubbles of various deflection profiles.8,17

For a general nanobubble, it is useful to recast Eqs. (1)–(3) into one equation by use of the Airy stress function, χ, defined by9,22

σxx=2χy2,σxy=2χyx,
(4)

and

σyy=2χx2.

Because of the equality of mixed partial derivatives, Eq. (3) is satisfied automatically. Thereby, any solutions we find using the Airy stress function will also be solutions of Eq. (3). Using Eqs. (1) and (2) and differentiating twice, we arrive at a fourth-order equation for the Airy stress function

22χ=E2hx22hy22hyx2.
(5)

Thus, the problem can be put in the form of an inhomogeneous biharmonic equation23 for the Airy stress function with a source function, the local Gaussian curvature, that is composed of products of the second-order partial derivatives of the measured topographic height multiplied by Young’s modulus. Computationally solving (5) amounts to defining the biharmonic operator ∇22 for the solution domain, applying boundary conditions, and inverting the operator. For this work, we employed a spectral collocation method,23–25 solving for the Airy stress function as a linear combination of Chebyshev polynomials. Full details of the calculation and discussion of the boundary conditions are located in the supplementary material.

Before applying the technique to the topography of real nanobubbles, we verified that the method correctly reproduces the strain for nanobubbles with circular symmetry, where analytical solutions have been derived.17Figure 1(c) shows three radial cross sections of the symmetric nanobubbles that were derived originally in Ref. 17. Figure 1(d) shows the corresponding radial strain for the exact solution (solid lines) vs that derived from the numerical solution of Eq. (5) (with boundary conditions defined by the topography). As can be seen, the numerical solutions reproduce the analytical solutions for each case of nanobubble topography.

We next aim to verify that our computational method correctly solves the FvK equations. To verify our solutions, it is necessary to measure the strain independently of the estimate obtained from the AFM topography. For the atomically thin TMDs, both the photoluminescence and Raman scattering signals shift to lower energy with increasing tensile strain.2,5,15,20,26,27 In principle, this strain dependence would then allow us to hyperspectrally map the local energies of the photoluminescence and Raman signals and use their relative shifts as a proxy for the local tensile strain. The extracted strain could then be contrasted with our calculation method that is outlined above. However, because of the nanometric size of many nanobubbles, with radii of order ∼100 nm or less, this experimental approach, namely its resolution, is severely limited by the diffraction limit of conventional optical microscopes, which imposes a resolution limit of ∼335 nm for the wavelengths used in this study.

To overcome this limitation, apertureless scanning near-field optical microscopy (a-SNOM) is a class of scanning-probe-based super-resolution imaging that can achieve resolutions far better than that of conventional optical microscopy.28 These techniques are powerful for investigating inelastic optical scattering of low-dimensional semiconductors and have been successfully harnessed to probe a variety of nanoscale phenomena in monolayer TMDs.29–33 For this work, we used tip-enhanced optical microscopy, which achieves nanoscopic resolution by exciting localized surface plasmons at the apex of a metallized AFM probe. The resulting local enhancement of the incident optical field can exceed factors of 103 on length scales down to a few nanometers or below in certain cases,34,35 enabling hyperspectral tip-enhanced photoluminescence and Raman scattering (TEPL36 and TERS,37 respectively) imaging with resolutions <10 nm. Thus, SNOM offers an ideal approach to investigate optical changes due to nanoscopic strain within single nanobubbles.

We have chosen to study nanobubbles in monolayer WS2 as the benchmark for our approach. Monolayer WS2 possesses excellent optical properties and its photoluminescence and Raman dependencies on strain are well-characterized both experimentally20 and theoretically,18 making it an ideal material for experimental verification of the strain estimation methodology that is presented here. Our sample was prepared by exfoliation of a bulk crystal of WS2 onto gold. A template-stripped gold substrate was used, which produces metallic films with very low surface roughness with rms values of <1 nm,38 allowing for excellent adhesion of the exfoliated WS2 onto the substrate. More details of the sample preparation are provided in Ref. 38. Furthermore, the use of the gold substrate also allows both TERS and TEPL to operate in “gap-mode,” which produces even larger field enhancements compared to dielectric or insulating substrates.28 

An excitation wavelength of 637.25 nm was used for the TERS and TEPL characterization measurements presented here. The incident excitation power of the laser was ∼70 μW. The laser wavelength is near resonant with the A exciton of unstrained monolayer WS2, which lies at ∼620 nm.32 While this choice of the excitation does limit the observable photoluminescence, the resonance condition greatly enhances the Raman scattering response39 and, as will be shown, the effect of strain is substantial enough that, over the nanobubble, significant TEPL is observed. Thus, this choice of excitation allows simultaneous measurement of the TEPL and TERS signals, eliminating the risk of erroneous correlations due to modification of the sample between measurements. All AFM and TERS/TEPL measurements were conducted on a customized commercial AFM and optical microscope (Horiba OmegaScope-R) that was coupled to a micro-Raman spectrometer (Horiba Xplora).

Figure 2 shows the tip-enhanced optical microscopy comprised of a sum of integrated TERS and TEPL signals [Fig. 2(a)] alongside a semi-contact AFM micrograph [Fig. 2(b)] of a field of nanobubbles of monolayer WS2 on gold. The field contains dozens of the nanobubbles, the majority of which have irregular, highly noncircular shapes. Additionally, several nanobubbles are separated by distances where the strain fields could potentially overlap, further complicating solutions of the FvK equations for a single nanobubble. As can clearly be seen by comparing Figs. 2(a) and 2(b), each resolvable nanobubble in the AFM micrograph has a substantial optical response. Notably, the emission intensity is non-uniform over the nanobubble surface with certain nanobubble regions exhibiting greater intensity than neighboring regions. A possible explanation for this localized emission within individual nanobubbles is funneling of excitons to points of greatest strain, which was theoretically predicted7 and recently measured2 in monolayer MoS2. We comment more on this possibility below.

FIG. 2.

Nano-optical and AFM characterization of an area of monolayer WS2 exfoliated on template-stripped gold. (a) Map of the integrated TEPL and TERS intensity from the monolayer WS2. (b) AFM micrograph of the same region revealing a dense population of nanobubbles of WS2 that show a one-to-one correspondence to the bright spots of emission in (a). The white arrow indicates the nanobubble featured in Fig. 3. Scale bars are 500 nm.

FIG. 2.

Nano-optical and AFM characterization of an area of monolayer WS2 exfoliated on template-stripped gold. (a) Map of the integrated TEPL and TERS intensity from the monolayer WS2. (b) AFM micrograph of the same region revealing a dense population of nanobubbles of WS2 that show a one-to-one correspondence to the bright spots of emission in (a). The white arrow indicates the nanobubble featured in Fig. 3. Scale bars are 500 nm.

Close modal

On inspecting the distribution of nanobubbles in Fig. 2(a), one observes that there are many small nanobubbles dispersed throughout the map. A closer look reveals that the larger nanobubbles have relatively few of these smaller nanobubbles immediately around them. We hypothesize that this behavior is the result of the attractive force between neighboring nanobubbles, which was predicted theoretically by Khestaonva et al.1: the greater strain field of a larger nanobubble can attract and merge with surrounding smaller bubbles, effectively clearing the area around it. While the data in Fig. 2 are not sufficient to prove this effect unambiguously, the consistency of this feature does imply an underlying physical mechanism that leads to a systematic pattern beyond random placements that occur during fabrication.

We next solve Eq. (5) using the AFM topography of a single large nanobubble in Fig. 2(b) and compare the result to the corresponding TEPL and TERS data from Fig. 2(a). For this test, we choose the largest nanobubble located in the top center of the field (indicated by the white arrow), which shows a relative lack of close neighbors. In Fig. 3(a), we show the resulting map of the trace of the strain tensor, using a Poisson ratio of 0.22 for 1L-WS2.40 We see that maximum strain of ∼0.85% is located at the center of the nanobubble, which is also at the point of maximum deflection when compared to the AFM topography [inset of Fig. 3(a)]. This correlation between deflection and strain is expected from the analysis of circularly symmetric plates that is presented in Fig. 1. However, in stark contrast to the predictions of a circular plate, our strain analysis predicts a secondary maximum located near the edge of the nanobubble.

FIG. 3.

Strain analysis of individual nanobubbles of monolayer WS2. (a) Trace of the strain tensor calculated from the second Föppl–von Kármán equation for the nanobubbles given in Fig. 1. The inset contains the corresponding AFM micrograph. (b) High-resolution map of the integrated TEPL and TERS signals from the nanobubble in (a). (c) Sample emission spectra from localized regions of the nanobubble [denoted by the same colored crosses in (b)]. The localized emission is comprised of the TERS signal, TEPL from strain shifted exciton states, and TEPL from localized exciton (LX) states. Scale bars are 50 nm.

FIG. 3.

Strain analysis of individual nanobubbles of monolayer WS2. (a) Trace of the strain tensor calculated from the second Föppl–von Kármán equation for the nanobubbles given in Fig. 1. The inset contains the corresponding AFM micrograph. (b) High-resolution map of the integrated TEPL and TERS signals from the nanobubble in (a). (c) Sample emission spectra from localized regions of the nanobubble [denoted by the same colored crosses in (b)]. The localized emission is comprised of the TERS signal, TEPL from strain shifted exciton states, and TEPL from localized exciton (LX) states. Scale bars are 50 nm.

Close modal

To confirm this prediction, independent signatures of such a secondary region of maximized strain should be identifiable as energetic shifts of the photoluminescence and Raman response that are colocalized with the strained region. The corresponding nano-optical data for this nanobubble are shown in Figs. 3(b) and 3(c). The integrated TEPL and TERS intensity is shown in Fig. 3(b), and select point spectra are shown in Fig. 3(c). Indeed, the reddest emission correlates with the area, where our model predicts greatest strain.

The maximum emission intensity occurs near the secondary maximum of strain. However, exciton funneling effects would result in regions of maximum strain exhibiting the brightest PL emission intensity, as the highest strained regions correspond to the deepest potential wells and thus the energetic minima for diffusing excitons.2,7 This apparent discrepancy can be accounted for when the excitation geometry of the gap-mode TEPL/TERS measurement is considered. The enhancement of the electric field due to the plasmonic probe decays exponentially as the size of the gap between the probe and underlying gold film is increased.41 Thus, for points of a nanobubble farther away from the gold substrate, the increased emission intensity due to exciton funneling is offset and eventually overcome by the large fall-off of excitation intensity and reduction in emission enhancement due to the increasing gap size. This effect is seen by diminished TERS intensity relative to the TEPL intensity at points toward the nanobubble peak in Fig. 3 as demonstrated by the yellow spectra that is recorded at the apex of the nanobubble. These reductions in the ratio between the TERS and TEPL intensities are signatures of a corresponding reduction in the excitation intensity due to this enlargement of the gap.

Finally, we show that our strain-analysis method can be applied over the large field of nanobubbles shown in Fig. 2. Such a large area is computationally prohibitive to analyze with atomistic simulations, but clearly cannot be analyzed with elastic plate theory alone. Figures 4(a)–4(c) show three square fields of nanobubbles from subregions of Fig. 2 with their corresponding maps of TEPL and TERS intensity for ease of reference (inset). Figure 4(a) includes the nanobubble analyzed in detail in Fig. 3(a). As can be seen by comparison of Figs. 4(a) and 3(a), the two solutions have nearly the same strain morphology, where the minor discrepancies are attributed to the fact that the solution of Fig. 3(a) implicitly assumed that the nanobubble is completely isolated [the larger field solution of Fig. 4(a) relaxes this assumption].

FIG. 4.

Strain analysis of larger areas with multiple nanobubbles of monolayer WS2. (a)–(c) The trace of the strain tensor calculated from the second Föppl–von Kármán equations subregions of the area shown in Fig. 2(a). The insets show the integrated intensity of the TEPL and TERS of the regions analyzed. (d) The energy of the measured E′ Raman mode of the monolayer WS2 as a function of tensile strain averaged over multiple nanobubbles. The error bars are estimated using the standard error of the mean. A linear fit is denoted in red and yields a slope of −1.75 cm−1/%, corresponding to a Grüneisen parameter 0.49 ± 0.23, agreeing well with other experimental values. All scale bars are 200 nm.

FIG. 4.

Strain analysis of larger areas with multiple nanobubbles of monolayer WS2. (a)–(c) The trace of the strain tensor calculated from the second Föppl–von Kármán equations subregions of the area shown in Fig. 2(a). The insets show the integrated intensity of the TEPL and TERS of the regions analyzed. (d) The energy of the measured E′ Raman mode of the monolayer WS2 as a function of tensile strain averaged over multiple nanobubbles. The error bars are estimated using the standard error of the mean. A linear fit is denoted in red and yields a slope of −1.75 cm−1/%, corresponding to a Grüneisen parameter 0.49 ± 0.23, agreeing well with other experimental values. All scale bars are 200 nm.

Close modal

To show the quantitative accuracy of the strain solution, we leverage the high number of nanobubbles to calculate the Grüneisen parameter for the E′ phonon mode of monolayer WS2. The Grüneisen parameter, γi, relates the change in frequency of a crystal vibration to a hydrostatic stretching (without any shear) for a single vibrational mode, i, and is defined by19 

γi=VωiωiV,
(6)

where V is the volume of the unit cell, and ωi is the angular frequency of the vibrational mode. Since the E′ mode in monolayer WS2 is confined to the plane of the material, we can express the Grüneisen parameter as the change in the frequency due to pure in-plane biaxial strain:20 

γi=1ωiωiϵ.
(7)

Here, ϵ refers to the fractional change in area due to stretching.

For even a perfectly circular nanobubble, the strain is in general not purely hydrostatic, but varies from a pure hydrostatic stretching at the nanobubble center to a pure shear at the nanobubble edge (see supplementary material as well as Ref. 17). Raman scattering signals from regions of the nanobubble off the center will therefore reflect a mix of a pure shear and a pure stretch. From the integrated TERS intensity in supplementary material Fig. 3, it is clear that the strongest Raman is from such intermediate regions, and therefore, the mode shifts are not described by the Grüneisen parameter alone.20 However, unlike the hydrostatic component, which provides a net tensile strain, shearing results in both positive and negative strains. When averaging over a large number of nanobubbles, we expect the strain to be approximately hydrostatic as the different signed shear components will average to 0. We show that this is indeed the case in supplementary material Fig. 2, which reports the average values of the shear components of the strain in the nanobubbles as a function of the hydrostatic strain. The mean values of the shear strain components are <0.1% when the data are sorted into bins that cover 0.16% tensile strain. Therefore, the strain when averaged over many nanobubbles is approximately pure stretching without shear, and the corresponding slope that describes how the average energy of the E′ mode depends on tensile strain under similar binning procedures should reproduce the Grüneisen parameter.

Figure 4(d) shows the mean Raman shift of the E′ mode vs the hydrostatic strain for six strain bins of width ∼0.16% (see supplementary material for further discussion on the chosen bin size). The error bars are the standard error of the mean. A linear fit to the data is plotted in red, which yields a corresponding slope of −1.75 cm−1/%. Inserting this value into Eq. (7), we obtain a Grüneisen parameter of 0.49 ± 0.23, which agrees well with the experimental values of 0.5 that was measured via a diamond anvil cell19 and 0.54 that was measured using uniaxial plate-bending techniques.20 

In conclusion, we have demonstrated a facile numerical recipe using the AFM topography to estimate the strain of nanobubbles in 2D materials that accurately estimates the strain at nanoscopic lengths, where inhomogeneities dramatically effect optoelectronic properties,29 but are too large for advanced atomistic calculations. This method is generally applicable to nanobubbles of arbitrary shape and can be easily applied to complex areas that are densely populated with both isolated and interacting nanobubbles. Using state-of-the-art nano-optical characterization tools, we show that the spatial shifts of the photoluminescence energy agree with the estimated strain of a single noncircular nanobubble. By combining the strain maps of dozens of irregular nanobubbles, we reproduce the Grüneisen parameter of the E′ in-plane Raman active mode for the monolayer WS2, showing the quantitative accuracy of our method for truly “nano” nanobubbles with radii on the order 10–100 nm. Although we have applied this method to nanobubbles, we note that Eq. (5) does not make an assumption of a “bubble” or of a trapped material, and the recipe is generally applicable to other deflections of 2D materials such as wrinkles. Because of the wide availability of powerful numerical solvers for partial differential equations and the ubiquity of AFM microscopes for topographic characterization, our approach provides a practical and accessible recipe for estimation of strain in low-dimensional materials.

See the supplementary material for methods and additional supporting data and figures.

The data that support the findings of this study are available from the corresponding author upon reasonable request. Our code for calculating the strain tensor components of an arbitrarily shaped nanobubble in a 2D semiconductor based on AFM topography data is available free of charge at https://github.com/T1618/Nanobubble-Strain-Tensor-Calc.

N.J.B. and P.J.S. gratefully acknowledge support from the National Science Foundation under Award No. NSF-1838403. J.W.K. gratefully acknowledges support from the National Science Foundation under Award Nos. NSF-1437450 and NSF-1363093. T.P.D. and P.J.S. acknowledge partial support from the Programmable Quantum Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under Award No. DE-SC0019443. D.J. acknowledges primary support for this work by the U.S. Army Research Office under Contract No. W911NF-19-1-0109 and the National Science Foundation (Grant No. DMR-1905853). Work at the Molecular Foundry was supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

1.
E.
Khestanova
,
F.
Guinea
,
L.
Fumagalli
,
A. K.
Geim
, and
I. V.
Grigorieva
, “
Universal shape and pressure inside bubbles appearing in van der Waals heterostructures
,”
Nat. Commun.
7
,
12587
(
2016
).
2.
A. V.
Tyurnina
etet al, “
Strained bubbles in van der Waals heterostructures as local emitters of photoluminescence with adjustable wavelength
,”
ACS Photonics
6
,
516
524
(
2019
).
3.
G. D.
Shepard
etet al, “
Nanobubble induced formation of quantum emitters in monolayer semiconductors
,”
2D Mater.
4
,
021019
(
2017
).
4.
T. P.
Darlington
 et al, “
Imaging strain-localized exciton states in nanoscale bubbles in monolayer WSe2 at room temperature
,”
Nat. Nanotech.
(to be published) (
2020
).
5.
D.
Lloyd
etet al, “
Band gap engineering with ultralarge biaxial strains in suspended monolayer MoS2
,”
Nano Lett.
16
,
5836
5841
(
2016
).
6.
E.
Blundo
 et al, “
Evidence of the direct-to-indirect band gap transition in strained two-dimensional WS2, MoS2, and WSe2
,”
Phys. Rev. Res.
2
,
012024
(
2020
).
7.
J.
Feng
,
X.
Qian
,
C.-W.
Huang
, and
J.
Li
, “
Strain-engineered artificial atom as a broad-spectrum solar energy funnel
,”
Nat. Photonics
6
,
866
872
(
2012
).
8.
Z.
Dai
 et al, “
Interface-governed deformation of nanobubbles and nanotents formed by two-dimensional materials
,”
Phys. Rev. Lett.
121
,
266101
(
2018
).
9.
L. D.
Landau
and
E. M.
Lifshiftz
,
Theory of Elasticity
(
Pergamon Press
,
1959
).
10.
H.
Hencky
, “
On the stress state in circular plates with vanishing bending stiffness
,”
Z. Math. Phys.
63
,
311
317
(
1915
).
11.
J.
Zabel
etet al, “
Raman spectroscopy of graphene and bilayer under biaxial strain: Bubbles and balloons
,”
Nano Lett.
12
,
617
621
(
2012
).
12.
S. P.
Koenig
,
N. G.
Boddeti
,
M. L.
Dunn
, and
J. S.
Bunch
, “
Ultrastrong adhesion of graphene membranes
,”
Nat. Nanotechnol.
6
,
543
546
(
2011
).
13.
D.
Lloyd
etet al, “
Adhesion, stiffness, and instability in atomically thin MoS2 bubbles
,”
Nano Lett.
17
,
5329
5334
(
2017
).
14.
M. R.
Rosenberger
etet al, “
Nano-“Squeegee” for the creation of clean 2D material interfaces
,”
ACS Appl. Mater. Interfaces
10
,
10379
10387
(
2018
).
15.
C.
Carmesin
etet al, “
Quantum-dot-like states in molybdenum disulfide nanostructures due to the interplay of local surface wrinkling, strain, and dielectric confinement
,”
Nano Lett.
19
,
3182
3186
(
2019
).
16.
Y.
Ma
etet al, “
Extended Hencky solution for the blister test of nanomembrane
,”
Extreme Mech. Lett.
22
,
69
78
(
2018
).
17.
H.
Rostami
,
F.
Guinea
,
M.
Polini
, and
R.
Roldán
, “
Piezoelectricity and valley chern number in inhomogeneous hexagonal 2D crystals
,”
npj 2D Mater. Appl.
2
,
15
(
2018
).
18.
R.
Frisenda
etet al, “
Biaxial strain tuning of the optical properties of single-layer transition metal dichalcogenides
,”
npj 2D Mater. Appl.
1
,
10
(
2017
).
19.
W.
Fan
 et al, “
Vibrational spectrum renormalization by enforced coupling across the van der Waals gap between MoS2 and WS2 monolayers
,”
Phys. Rev. B
92
,
241408
(
2015
).
20.
A. M.
Dadgar
etet al, “
Strain engineering and Raman spectroscopy of monolayer transition metal dichalcogenides
,”
Chem. Mater.
30
,
5148
5155
(
2018
).
21.
S.
Timoshenko
and
S.
Woinowsky-Krieger
,
Theory of Plates and Shells
(
McGraw-Hill
,
1959
).
22.
H. S.
Seung
and
D. R.
Nelson
, “
Defects in flexible membranes with crystalline order
,”
Phys. Rev. A
38
,
1005
1018
(
1988
).
23.
N.
Mai-Duy
and
R. I.
Tanner
, “
A spectral collocation method based on integrated Chebyshev polynomials for two-dimensional biharmonic boundary-value problems
,”
J. Comput. Appl. Math.
201
,
30
47
(
2007
).
24.
L. N.
Trefethen
,
Spectral Methods in MATLAB
(
Society for Industrial and Applied Mathematics
,
2000
).
25.
N.
Jewell
, MATLAB Central File Exchange,
2020
.
26.
A.
Steinhoff
,
M.
Rösner
,
F.
Jahnke
,
T. O.
Wehling
, and
C.
Gies
, “
Influence of excited carriers on the optical and electronic properties of MoS2
,”
Nano Lett.
14
,
3743
3748
(
2014
).
27.
S. B.
Desai
etet al, “
Strain-induced indirect to direct bandgap transition in multilayer WSe2
,”
Nano Lett.
14
,
4592
4597
(
2014
).
28.
D. F.
Ogletree
etet al, “
Near-field imaging: Revealing optical properties of reduced-dimensionality materials at relevant length scales
,”
Adv. Mater.
27
,
5692
(
2015
).
29.
W.
Bao
etet al, “
Visualizing nanoscale excitonic relaxation properties of disordered edges and grain boundaries in monolayer molybdenum disulfide
,”
Nat. Commun.
6
,
7993
(
2015
).
30.
K.-D.
Park
etet al, “
Hybrid tip-enhanced nanospectroscopy and nanoimaging of monolayer WSe2 with local strain control
,”
Nano Lett.
16
,
2621
2627
(
2016
).
31.
A.
Krayev
etet al, “
Dry transfer of van der Waals crystals to noble metal surfaces to enable characterization of buried interfaces
,”
ACS Appl. Mater. Interfaces
11
,
38218
38225
(
2019
).
32.
C.
Kastl
etet al, “
The important role of water in growth of monolayer transition metal dichalcogenides
,”
2D Mater.
4
,
021024
(
2017
).
33.
K. K. H.
Smithe
etet al, “
Nanoscale heterogeneities in monolayer MoSe2 revealed by correlated scanning probe microscopy and tip-enhanced Raman spectroscopy
,”
ACS Appl. Nano Mater.
1
,
572
579
(
2018
).
34.
M.
Richard-Lacroix
,
Y.
Zhang
,
Z.
Dong
, and
V.
Deckert
, “
Mastering high resolution tip-enhanced Raman spectroscopy: Towards a shift of perception
,”
Chem. Soc. Rev.
46
,
3922
3944
(
2017
).
35.
A.
Bhattarai
etet al, “
Tip-enhanced Raman nanospectroscopy of smooth spherical gold nanoparticles
,”
J. Phys. Chem. Lett.
11
,
1795
1801
(
2020
).
36.
E. J.
Cheng
,
L.
Novotny
, and
X. S.
Xie
, “
Near-field fluorescence microscopy based on two-photon excitation with metal tips
,”
Phys. Rev. Lett.
82
,
4014
4017
(
1999
).
37.
R. M.
Stöckle
,
Y. D.
Suh
,
V.
Deckert
, and
R.
Zenobi
, “
Nanoscale chemical analysis by tip-enhanced Raman spectroscopy
,”
Chem. Phys. Lett.
318
,
131
136
(
2000
).
38.
D.
Jariwala
,
A.
Krayev
,
J.
Wong
,
A. E.
Robinson
,
M. C.
Sherrott
,
S.
Wang
,
G.-Y.
Liu
,
M.
Terrones
, and
H. A.
Atwater
, “
Nanoscale doping heterogeneity in few-layer WSe2 exfoliated onto noble metals revealed by correlated SPM and TERS imaging
,”
2D Mater.
5
,
035003
(
2018
).
39.
E.
del Corro
etet al, “
Atypical exciton–phonon interactions in WS2 and WSe2 monolayers revealed by resonance Raman spectroscopy
,”
Nano Lett.
16
,
2363
2368
(
2016
).
40.
K.
Liu
etet al, “
Elastic properties of chemical-vapor-deposited monolayer MoS2, WS2, and their bilayer heterostructures
,”
Nano Lett.
14
,
5097
5103
(
2014
).
41.
L.
Chen
and
B.
Hecht
,
Principles of Nano-Optics
(
Cambridge University Press
,
2012
).

Supplementary Material