Facile and quantitative estimation of strain in nanobubbles with arbitrary symmetry in 2D semiconductors verified using hyperspectral nano-optical imaging

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.¹ 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.⁵ The inhomogeneity of the strain has also been shown to induce complex phenomena such as exciton funneling^2,6,7 and single-photon emission at cryogenic temperatures,³ 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.¹⁰ 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.¹¹ Hencky’s solution has also been employed to estimate the van der Waals adhesion energies of both graphene¹² and MoS₂,¹³ as well as the shift of the optical bandgap of MoS₂ under biaxial strain.⁵ 

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.¹⁴ For nanobubbles of this type, Khestanova et al. solved the FvK equations for circularly symmetric nanobubbles in MoS₂, 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 MoS₂ nanobubbles.² 

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,¹⁴ 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.¹⁵ 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.⁵ By contrast, more sophisticated calculations of MoS₂ 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.¹⁵ 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 WSe₂ nanobubbles indeed correspond to the nanoscale distributions of strain that are predicted by the atomistic theory,⁴ 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.¹⁸ 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 WS₂ 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 WS₂, 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.¹ 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

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,

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:

and

where u_x and u_y are in-plane body displacements, and h = u_z 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:

and

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 by^9,22

and

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

Thus, the problem can be put in the form of an inhomogeneous biharmonic equation²³ 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 ∇²∇² 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.¹⁷ Figure 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.²⁸ 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 10³ on length scales down to a few nanometers or below in certain cases,^34,35 enabling hyperspectral tip-enhanced photoluminescence and Raman scattering (TEPL³⁶ and TERS,³⁷ 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 WS₂ as the benchmark for our approach. Monolayer WS₂ possesses excellent optical properties and its photoluminescence and Raman dependencies on strain are well-characterized both experimentally²⁰ and theoretically,¹⁸ 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 WS₂ onto gold. A template-stripped gold substrate was used, which produces metallic films with very low surface roughness with rms values of <1 nm,³⁸ allowing for excellent adhesion of the exfoliated WS₂ 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.²⁸ 

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 WS₂, which lies at ∼620 nm.³² While this choice of the excitation does limit the observable photoluminescence, the resonance condition greatly enhances the Raman scattering response³⁹ 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 WS₂ 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 predicted⁷ and recently measured² in monolayer MoS₂. We comment more on this possibility below.

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.¹: 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-WS₂.⁴⁰ 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.

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.⁴¹ 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].

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 WS₂. 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 by¹⁹ 

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 WS₂ 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:²⁰ 

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.²⁰ 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⁻¹/%. 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 cell¹⁹ and 0.54 that was measured using uniaxial plate-bending techniques.²⁰ 

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,²⁹ 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 WS₂, 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.