Preparation and Modeling of Graphene Bubbles to Obtain Strain-Induced Pseudomagnetic Fields

It has been both theoretically predicted and experimentally demonstrated that strain can effectively modulate the electronic states of graphene sheets through the creation of a pseudomagnetic field (PMF). Pressurizing graphene sheets into bubble-like structures has been considered a viable approach for the strain engineering of PMFs. However, the bubbling technique currently faces limitations such as long manufacturing time, low durability, and challenges in precise control over the size and shape of the pressurized bubble. Here, we propose a rapid bubbling method based on an oxygen plasma chemical reaction to achieve rapid induction of out-of-plane deflections and in-plane strains in graphene sheets. We introduce a numerical scheme capable of accurately resolving the strain field and resulting PMFs within the pressurized graphene bubbles, even in cases where the bubble shape deviates from perfect spherical symmetry. The results provide not only insights into the strain engineering of PMFs in graphene but also a platform that may facilitate the exploration of the strain-mediated electronic behaviors of a variety of other 2D materials.


Introduction
Two-dimensional (2D) materials have garnered significant attention owing to their unique structure-a 2D hexagonal lattice composed of carbon atoms-which imparts remarkable electronic properties to them [1][2][3][4][5].In addition, the atomical thinness of 2D materials endows them with exceptional flexibility and mechanical deformability, further enabling the manipulation of their electronic structure through strain engineering methods [5][6][7][8][9][10].For instance, it has been found that out-of-plane deformation can induce local changes in the electronic properties of graphene [9,10], leading to a shift in the electron's Bloch wave vector, akin to the behavior observed in the presence of an external magnetic field [4,11,12].This phenomenon has been termed strain-induced pseudomagnetic field (PMF) [4,11,13].The intensity of such PMFs (which can be even greater than 300 Tesla) can surpass by several orders of magnitude that of the external magnetic field generated with a superconducting magnet, opening up exciting avenues for exploring novel quantum phenomena in graphene systems [4,14,15].More broadly, strain engineering is a valuable tool for modifying the conductivity, magnetism, and various other electronic characteristics of 2D materials [16][17][18].
Here, we first briefly discuss the advantages and disadvantages of typical bubbling methods in terms of durability and the controllability of the size and shape of the resulting bubbles.We then present a fast, simple bubbling technique relying on oxygen plasma chemical reactions.This method effectively demonstrates the rapid generation of out-ofplane deflections of graphene sheets.To accurately compute the strain field and resulting PMFs within the pressurized graphene bubbles, we develop a numerical scheme that is applicable to bubbles of any shape.The effect of the slippage of the graphene sheet over its supporting substrate on the PMFs is discussed.The experimental approach to inducing graphene bubbles and the computational method for analyzing the resulting strain fields can be applied to a variety of other 2D materials for the study of their fascinating strainmediated electronic and optoelectronic behaviors.

A Brief Overview of Bubbling Methods
Nanobubbles can be created by inserting gas molecules between 2D materials and a substrate.Specifically, protons with specific energy can penetrate the 2D material and form gas molecules between layers or between layers and substrates, resulting in dome-shaped bubbles.For instance, Tedeschi et al. [42] and He et al. [43] reported that, via irradiating bulk crystals of transition metal dichalcogenides (TMDs) and hexagonal boron nitride (hBN) with low-energy hydrogen ions or hydrogen plasma, hydrogen ions could penetrate the uppermost layer to form gas-condensed bubbles filled with H 2 .Yasuda et al. utilized the permeability of graphene defect sites to hydrogen protons and the electrochemical hydrogen evolution reaction to construct hydrogen bubbles instead of using proton irradiation [44].The height and radius of hydrogen proton-induced bubbles can be controlled within a range from a few nanometers to several hundred nanometers or even micrometers.In addition, the medium can also comprise other particles such as Ne + , Xe + , Ar + , and He + , but the resulting bubbles usually exhibit a height range that spans from subnanometer to several nanometers [45][46][47][48].Controlling the height, shape, and position of the bubbles using these methods above has proven challenging [49].To overcome this challenge, Polimeni and coworkers have recently achieved selective proton irradiation via creating a layer of H-opaque mask followed by proton irradiation, which led to a top array of single-layer controlled-size/position circles containing pressurized hydrogen [42,50,51].
In fact, in earlier stages of relevant strain engineering studies, PMFs from graphene bubbles were formed spontaneously via directly transferring the 2D material sheet onto a target substrate [4,39].During the transfer process, gas, liquid, or foreign contaminants (such as residues of organic matter used for transfer) could become trapped on the surface of either the 2D material or the substrate [41,[52][53][54][55][56][57][58][59][60][61][62].After contact, these substances could become squeezed and aggregated at the interface due to the attractive van der Waals (vdW) interaction at the 2D material/substrate interface.It has also been reported that this type of bubble can be formed in the presence of defects within the material or at the material interface, through which gas and liquid molecules may migrate spontaneously to the layer/layer interface [63][64][65][66][67][68][69][70].Recent attempts at the fabrication of 2D material bubbles have involved the intentional introduction of nanoparticles at the interface, giving rise to artificial means to control the process [71][72][73].
The density and size of such spontaneous bubbles tend to be random.In order to control the formation of bubbles, extensive efforts have employed external energy input to facilitate the intercalation of extra water or gas molecules into the interface of the 2D materials.For instance, an electron beam can induce the dissociation of water sandwiched between graphene layers to produce hydrogen gas, thereby triggering gasfilled graphene bubbles [74,75].This also illustrates the potential of graphene as an effective chemical reactor and hydrogen storage container.Proton irradiation or chemical etching can also cause a chemical reaction in silicon dioxide substrate beneath the graphene, leading to bubble formation [76].Applying a voltage allows the water within the interface to decompose into gas, resulting in bubbles with controllable size and shape [77,78].Laser irradiation is another method that can activate specific chemicals.For example, chlorine trifluoride (ClF3) or liquid nitrogen (LN) embedded between graphite layers can induce bubbling under laser stimulation [79,80].
The construction of bubbles using 2D material can also be intentionally achieved through introducing liquid at the vdW interface [53,64,66,[81][82][83]. Pantelis et al. demonstrated the construction of water-filled blisters at the MoS 2 /graphene interface through adjusting the air humidity, reaching heights of nearly 30 nm [81].Although the introduction of liquid at the layer/substrate interface during wet transfer is unintentional, it provides a powerful platform to study the chemical and physical behavior of liquids in nanoscale confined spaces.Additionally, blisters can be manipulated through changing liquid composition and volume.Vasu et al. utilized this method to encapsulate salts or molecular solutions between graphene/graphene or graphene/graphite interfaces [84].
Methods utilizing pressure chambers to controllably prepare bubbles have increasingly been utilized for studying the physical, chemical, and mechanical properties of strained 2D materials.One such method was first reported by Bunch et al. [85].Briefly, the 2D material was mechanically exfoliated or transferred onto a substrate with pre-patterned microcavities and then, the assembly was placed into a high-pressure chamber [85].Within the chamber, gas molecules diffused along the interface between the material and the substrate, as well as through the pores of the substrate, until the pressure equalized (usually 4-6 days but highly dependent on the gas used) [35,38]; a pressure difference was thus created after taking out the sample, which induced gas-pressurized 2D material bubbles [86][87][88][89][90][91].In addition, electrostatic force can be generated via applying a voltage, eliminating the need for a high-voltage chamber to drive the 2D material to bulge [92].
Table 1 presents a summary of the typical height, lateral size, controllability of shape and position, and durability of bubbles induced by the aforementioned methods.It should be noted that the table does not include all methods, due to unavailable data.Material thickness is disregarded, and durability is measured through the absence of a significant height decrease.The data rely solely on reported values, potentially underestimating the upper limits of corresponding methods.It is clear from Table 1 that there remain shortcomings such as lengthy bubble manufacturing cycles, complex device setups, and poor durability.These limitations hinder their utility in strain engineering applications as well as nanodevice applications including valleytronic devices, pressure sensors, and so on [2,[93][94][95][96].

Experimental Method
Here, we report an alternative method that effectively combines the advantages of spontaneous bubbles and pressurized bubbles.Figure 1 shows a schematic diagram of the oxygen plasma-assisted preparation of graphene bubbles proposed in this work.Based on micro-nano processing technology, photolithography and reactive ion etching (RIE, ETCHLAB 200, SENTECH Inc., Berlin, Germany) were employed to fabricate an array of circular holes on the surface of a silicon wafer (with a silicon oxide layer thickness of 300 nm).Before use, the substrate with pre-patterned microcavities was subjected to preprocessing via oxygen plasma with a power of 100 W for 5 min.Note that oxygen plasma treatment (Plasma cleaner CPC-A, CIF International Group Co., Ltd, Compton, CA, USA) has previously been used to enhance the quality of ohmic contact via effectively cleaning the surface of 2D materials [62].Next, the graphene was micro-mechanically exfoliated onto the treated substrate via Scotch tape (3M Co., Cottage Grove, MN, USA).Subsequently, the assembly was heated at 90 • C for 6 min to reduce the stickiness of the tape and increase the yield of the few-layer graphene.Note that this setting was flexible-we found heating at 90-110 • C for 4-7 min worked very well.Upon immediate removal of the tape after heating, graphite flakes of varying thicknesses remained on the substrate, and graphene bubbles of different heights and different layers were formed spontaneously.Finally, we used scanning electron microscopy (SEM, Quanta 600, FEI Inc., Hillsboro, OR, USA) to characterize the topology of the prepared samples.In the upper right image in Figure 1, bubbles can be seen forming on the holes covered by graphene.This observation was consistent across the few-layer graphene samples examined in other locations, confirming the effectiveness of our method and highlighting its advantages for mass production.Here, we demonstrate only three-layer graphene sheets, while emphasizing that this method should also be effective for other 2D materials with various numbers of layers.
The underlying mechanism is as follows.Pre-patterned substrates with microcavities invariably accumulate organic contaminations during micro/nanofabrication and subsequent storage [97,98].Ultraviolet light produced in the plasma is highly effective in breaking down most organic bonds present in surface contaminants.This process aids in the breakdown of organic contaminations.At the same time, these organic compounds undergo vigorous oxidative chemical reactions when exposed to energetic oxygen species created in the plasma, as shown in Figure 1.Specifically, oxygen molecules undergo collision, dissociation, ionization, and recombination effects under the influence of high-speed electrons and the photoelectric effect, thereby generating oxygen plasma.Oxygen plasma comprises various particles, including O + 2 , O − , O − 2 , O − 3 , and excited-state energetic oxygen species such as O * 2 , among others.Notably, excited oxygen species can react with the organic layer on the substrate surface, resulting in the following chemical reactions: organic , leading to the spontaneous pressurization of the atomic sheet that seals these molecules into the cavity.The substrate with pre-patterned microcavities is pre-processed with oxygen plasma, after which the graphene is micro-mechanically exfoliated onto the treated substrate.Peeling the tape from the heated assembly can lead to the spontaneous bubbling of the transferred graphene sheet immediately at these prepatterned cavities.The upper right panel shows an SEM image capturing graphene bubbles at a particular deflection angle.
The underlying mechanism is as follows.Pre-patterned substrates with microcavities invariably accumulate organic contaminations during micro/nanofabrication and subsequent storage [97,98].Ultraviolet light produced in the plasma is highly effective in breaking down most organic bonds present in surface contaminants.This process aids in the breakdown of organic contaminations.At the same time, these organic compounds undergo vigorous oxidative chemical reactions when exposed to energetic oxygen species created in the plasma, as shown in Figure 1.Specifically, oxygen molecules undergo collision, dissociation, ionization, and recombination effects under the influence of highspeed electrons and the photoelectric effect, thereby generating oxygen plasma.Oxygen plasma comprises various particles, including O , O , O , O , and excited-state energetic oxygen species such as O * , among others.Notably, excited oxygen species can react with the organic layer on the substrate surface, resulting in the following chemical reactions: organic O * → CO H O [99], leading to the spontaneous pressurization of the atomic sheet that seals these molecules into the cavity.
We note that oxygen plasma treatment is commonly employed for surface cleaning and substrate modification.However, due to the unique structure of the substrate, organic matter is adsorbed on both the substrate surface and the inner wall of cavities, albeit without uniformity.When oxygen plasma bombards the surface, the chemical reaction within the cavities persists longer than on the surface itself.This phenomenon arises from the relatively confined spatial characteristics of cavities, which concentrate chemical reactions.Furthermore, organic contaminations on the surface might experience physical removal due to the impact of high-energy particles, which is similar to the effect of sandblasting.In contrast, the depth within the microcavity limits the physical removal of contamination, thereby ensuring the chemical reactions occurring inside the cavity.Notably, although the primary purpose of heating is to remove adhesive and diminish tape-material adhesion, heating is also likely to make the chemical reaction in the cavity more complete, leading to increased gas filling into graphene bubbles.

Surface Topography
Atomic force microscopy (AFM, Brucker Multimode 8) was employed to characterize the morphology of bubbles.The scanning was conducted in tapping mode with a scanning frequency of 0.7 Hz.An AFM cantilever featuring a silicon tip, with a resonance frequency of 150 kHz and a stiffness of 2.8 N/m, was utilized for the measurements.Figure 2a illustrates the height of bubbles fabricated via the oxygen plasma-assisted method, with a bubble height of about 130 nm.A comparative analysis with bubbles prepared using a bulge The substrate with pre-patterned microcavities is pre-processed with oxygen plasma, after which the graphene is micro-mechanically exfoliated onto the treated substrate.Peeling the tape from the heated assembly can lead to the spontaneous bubbling of the transferred graphene sheet immediately at these prepatterned cavities.The upper right panel shows an SEM image capturing graphene bubbles at a particular deflection angle.
We note that oxygen plasma treatment is commonly employed for surface cleaning and substrate modification.However, due to the unique structure of the substrate, organic matter is adsorbed on both the substrate surface and the inner wall of cavities, albeit without uniformity.When oxygen plasma bombards the surface, the chemical reaction within the cavities persists longer than on the surface itself.This phenomenon arises from the relatively confined spatial characteristics of cavities, which concentrate chemical reactions.Furthermore, organic contaminations on the surface might experience physical removal due to the impact of high-energy particles, which is similar to the effect of sandblasting.In contrast, the depth within the microcavity limits the physical removal of contamination, thereby ensuring the chemical reactions occurring inside the cavity.Notably, although the primary purpose of heating is to remove adhesive and diminish tape-material adhesion, heating is also likely to make the chemical reaction in the cavity more complete, leading to increased gas filling into graphene bubbles.

Surface Topography
Atomic force microscopy (AFM, Brucker Multimode 8) was employed to characterize the morphology of bubbles.The scanning was conducted in tapping mode with a scanning frequency of 0.7 Hz.An AFM cantilever featuring a silicon tip, with a resonance frequency of 150 kHz and a stiffness of 2.8 N/m, was utilized for the measurements.Figure 2a illustrates the height of bubbles fabricated via the oxygen plasma-assisted method, with a bubble height of about 130 nm.A comparative analysis with bubbles prepared using a bulge device revealed a conspicuous distinction.Particularly, bubbles prepared using a bulge device (similar to that in [85]) exhibited numerous minuscule bubbles at the graphene/substrate interface, as shown in Figure 2a,b.In contrast, these minuscule bubbles were absent in bubbles prepared using the oxygen plasma-assisted method.This difference confirmed the essential difference in the principles of the two bubbling methods.When using a bulge device, gas molecules are required to diffuse into the circular cavity through the channel between the material and the substrate.This diffusion causes gas molecules to accumulate at the interface, especially in the local contaminated areas, forming many nanobubbles and reducing interfacial adhesion.Furthermore, these smaller nanobubbles act as pathways for further gas diffusion, compromising the stability during the decay stage of the bubble.Conversely, the method proposed in this study effectively addresses this issue through directly generating gas from within the cavity, as shown in Figure 2d.Good interfacial integration should be able to slow down the escape of gas molecules.In addition, since bubble formation occurs through a rapid reaction, the preparation time of 2D materials bubbles is significantly reduced from over 4-6 days [13,35] to a few minutes.Therefore, this preparation method is high-throughput, with the number and size well controlled through the pre-patterning of the substrate, which can even be rectangular or other irregular shapes [37].
ing many nanobubbles and reducing interfacial adhesion.Furthermore, these smaller nanobubbles act as pathways for further gas diffusion, compromising the stability during the decay stage of the bubble.Conversely, the method proposed in this study effectively addresses this issue through directly generating gas from within the cavity, as shown in Figure 2d.Good interfacial integration should be able to slow down the escape of gas molecules.In addition, since bubble formation occurs through a rapid reaction, the preparation time of 2D materials bubbles is significantly reduced from over 4-6 days [13,35] to a few minutes.Therefore, this preparation method is high-throughput, with the number and size well controlled through the pre-patterning of the substrate, which can even be rectangular or other irregular shapes [37].(c,d) bubble prepared with the oxygen plasma-assisted method reported in this work.Note that both methods can lead to 1-3 layer graphene bubbles of heights ranging from 0 nm to 200 nm and a diameter of 3.5 µm.

Deflating Behavior
It has been reported that even the smallest gas molecules (such as helium) and ions (such as lithium) cannot penetrate defect-free graphene [90].At the graphene/SiO2 interface, however, bubbles prepared through gas pressurization via the bulge device can release most of the filled gas within 10 h for hydrogen-filled bubbles and within 7 days for nitrogen-filled bubbles [100,101].In cases of strong interface adhesion, gas molecule diffusion is slow and mainly occurs through channels within the silicon dioxide.We then  (c,d) bubble prepared with the oxygen plasma-assisted method reported in this work.Note that both methods can lead to 1-3 layer graphene bubbles of heights ranging from 0 nm to 200 nm and a diameter of 3.5 µm.

Deflating Behavior
It has been reported that even the smallest gas molecules (such as helium) and ions (such as lithium) cannot penetrate defect-free graphene [90].At the graphene/SiO 2 interface, however, bubbles prepared through gas pressurization via the bulge device can release most of the filled gas within 10 h for hydrogen-filled bubbles and within 7 days for nitrogen-filled bubbles [100,101].In cases of strong interface adhesion, gas molecule diffusion is slow and mainly occurs through channels within the silicon dioxide.We then investigated the durability of bubbles prepared via the bugle device and the method reported in this work.We focused on four sets of bubble samples.
Over the subsequent 31 days, we conducted AFM scans on the same samples using consistent scanning parameters and a cantilever.Our results show that the height of the three-layer thickness graphene bubbles prepared using the oxygen plasma-assisted method decreased at different rates, but to a smaller extent, as shown in Figure 3.As a comparison, bubbles of the same thickness prepared using a bulge device experienced rapid height reduction, nearly flattening after 5 days.Based on the observation in Figure 2, it is reasonable to hypothesize that the difference attributed to the fast gas-release channel caused by adhesion imperfections at the graphene/SiO 2 interface was largely suppressed in bubbles fabricated via our reported method.
method decreased at different rates, but to a smaller extent, as shown in Figure 3.As a comparison, bubbles of the same thickness prepared using a bulge device experienced rapid height reduction, nearly flattening after 5 days.Based on the observation in Figure 2, it is reasonable to hypothesize that the difference attributed to the fast gas-release channel caused by adhesion imperfections at the graphene/SiO2 interface was largely suppressed in bubbles fabricated via our reported method.

Figure 3.
Height measurements of 3-layer thickness bubbles measured over 31 days.All bubbles showed signs of gradual deflation, indicating that the gas of the bubbles was able to escape through the graphene/SiO2 interface but at a much slower rate in bubbles prepared via the oxygen plasmaassisted method than the direct bulging method.For the same bubble, the height was measured along three different line scanning directions, and the error was found within 1 nm.

Calculation of Strain and Pseudomagnetic Fields
We then moved on to assessment of the strain and potential pseudomagnetic fields in these prepared bubbles.Strain measurement within nanobubbles of moderate deflections can be achieved through utilizing approximate theories such as Föppl-von Kármán (FvK) equations.To solve such nonlinear partial differential equations, we followed the numerical scheme reported in [102], where a numerical approach based on the Chebyshev spectral method was used to solve the strain tensor of arbitrary symmetric 2D material bubbles.This method relied solely on morphological data obtained from AFM and the Poisson's ratio of the 2D material.We made modifications to the boundary conditions in this numerical scheme to match the practical situation in typical 2D material bubbles (to be discussed shortly) and then calculate the strain components in the 2D material sheet: where  and  are the 2D Cartesian coordinates,  and  are the in-plane displacements (to be solved), and ℎ is the out-of-plane deflection of the sheet (that can be directly measured from the bubble).Given the uncertainty regarding slip conditions at the graphene/SiO2 interface [39], we focused on two limiting cases (i.e., no-shear and no-slip conditions) to solve these strain components as well as the corresponding pseudomagnetic fields.In particular, an effective gauge field can be induced from the in-plane strains: via oxygen plasma-assisted via a bulge device Figure 3. Height measurements of 3-layer thickness bubbles measured over 31 days.All bubbles showed signs of gradual deflation, indicating that the gas of the bubbles was able to escape through the graphene/SiO 2 interface but at a much slower rate in bubbles prepared via the oxygen plasmaassisted method than the direct bulging method.For the same bubble, the height was measured along three different line scanning directions, and the error was found within 1 nm.

Calculation of Strain and Pseudomagnetic Fields
We then moved on to assessment of the strain and potential pseudomagnetic fields in these prepared bubbles.Strain measurement within nanobubbles of moderate deflections can be achieved through utilizing approximate theories such as Föppl-von Kármán (FvK) equations.To solve such nonlinear partial differential equations, we followed the numerical scheme reported in [102], where a numerical approach based on the Chebyshev spectral method was used to solve the strain tensor of arbitrary symmetric 2D material bubbles.This method relied solely on morphological data obtained from AFM and the Poisson's ratio of the 2D material.We made modifications to the boundary conditions in this numerical scheme to match the practical situation in typical 2D material bubbles (to be discussed shortly) and then calculate the strain components in the 2D material sheet: where x and y are the 2D Cartesian coordinates, u x and u y are the in-plane displacements (to be solved), and h is the out-of-plane deflection of the sheet (that can be directly measured from the bubble).
Given the uncertainty regarding slip conditions at the graphene/SiO 2 interface [39], we focused on two limiting cases (i.e., no-shear and no-slip conditions) to solve these strain components as well as the corresponding pseudomagnetic fields.In particular, an effective gauge field can be induced from the in-plane strains: In the low-energy Dirac Hamiltonian, the coupling constant A PMF ≈ 7 µm • T can be further linked to the hopping energy, Fermi velocity, and electron charge [103].This gauge field induces shifts in the Dirac cones of graphene at points K and K' in opposite directions, akin to the effect of a perpendicularly applied magnetic field.These PMFs can be associated with the strain-induced gauge field as follows:

Numerical Scheme
Due to the significantly smaller thickness of the graphene sheet compared with the dimensions of the film and the moderate deflection of graphene sheets used in the production of typical bubbles, the material law remains linear and the classical large-deflection thin plate theory (i.e., FvK equations) can be applied to establish control equations for in-plane displacement, strain, and stress [104].The equilibrium equation and constitutive equation can be written as: where Y is the in-plane Young's modulus (i.e., Young's modulus times thickness), ν is Poisson's ratio, µ is the in-plane shear modulus.We used the Airy stress function χ to express stresses as follows: so that equilibrium Equation ( 4) was satisfied automatically.Since the in-plane strain needs to satisfy the coordination relationship, the strain is expressed in the form of a stress function and brought into the coordination equation to obtain the compatibility equation of the FvK equations in the following form: The above equation is equivalent to the following two coupled equations: Following our experiments, we considered a square region characterized by a central nanobubble with a length l as the computational region, as illustrated in Figure 4a.We assumed that the film was flat within a small neighborhood at the boundaries of the square region, i.e., at these boundaries, the partial derivatives of the height along the x and y directions were zero.Since the boundary was far enough away from the center of the bubble, we assumed that the films in the neighborhood were under the same simple stress state.At the same time, we assumed that the tension applied at the film boundary was constant.Therefore, in the boundary neighborhood, the stress tensor of the film can be expressed as: where T xx , T yy , and T xy represent the constant normal and shear stresses applied in the x and y directions, respectively.We note that the far-field stresses were set to be zero in [102].At the boundary, the relationship between stress components and the stress function χ is established as: This leads to the expression for φ, the Laplacian of χ: where  ,  , and  represent the constant normal and shear stresses applied in the  and  directions, respectively.We note that the far-field stresses were set to be zero in [102].At the boundary, the relationship between stress components and the stress function  is established as: This leads to the expression for , the Laplacian of : Considering these nontrivial boundary conditions, we integrate the stress components to derive  ,  and  ,  : ,   ,    , Using  ,  and  ,  , we find  ,  , where       represents the Laplace equation solution: where  ,  , and  are integration constants.For the stress/strain field distribution, we focused on the second-order partial derivatives of  , hence ignoring the second-order terms.Thus, at the boundary, we simplify  ,  to: The height information of a nanofilm can be extracted from AFM scanning images.However, these images often contain noise, resulting in an irregular height function with numerous artifacts.To address this, Gaussian filtering was applied to the interpolated data obtained from the original AFM scans.As shown in Figure 4b,c, this filtering process effectively removed extraneous noise, resulting in a smoother image and more accurate calculation results.
Following [102], we applied the Chebyshev spectral method to solve the FvK equations.This method relies on Chebyshev polynomials as basis functions, enabling accurate approximation of functions with fewer grid points compared to other numerical methods.Considering these nontrivial boundary conditions, we integrate the stress components to derive A(x, y) and B(x, y): Using A(x, y) and B(x, y), we find χ(x, y), where χ 0 = c 1 x + c 2 y + c 3 represents the Laplace equation solution: where c 1 , c 2 , and c 3 are integration constants.For the stress/strain field distribution, we focused on the second-order partial derivatives of χ, hence ignoring the second-order terms.Thus, at the boundary, we simplify χ(x, y) to: The height information of a nanofilm can be extracted from AFM scanning images.However, these images often contain noise, resulting in an irregular height function with numerous artifacts.To address this, Gaussian filtering was applied to the interpolated data obtained from the original AFM scans.As shown in Figure 4b,c, this filtering process effectively removed extraneous noise, resulting in a smoother image and more accurate calculation results.
Following [102], we applied the Chebyshev spectral method to solve the FvK equations.This method relies on Chebyshev polynomials as basis functions, enabling accurate approximation of functions with fewer grid points compared to other numerical methods.The Chebyshev spectral method represents functions as finite series expansions of smooth basis functions.In this approach, differentiation operations are transformed into algebraic operations, simplifying the solution process.Taking the [−1, 1] × [−1, 1] interval as an example, with Chebyshev points defined as: x i = cos(iπ/N), i = 0, 1, 2, . . ., N andy j = cos(jπ/N), j = 0, 1, 2, . . ., N for convenience, we may rewrite the Laplacian operator as: where D N is an (N + 1) × (N + 1) matrix, and ⊗ denotes the Kronecker product.As a result, the FvK equations can be rewritten as: where These equations are then solved through successively taking the inverse of the modified matrices.It should be noted here that this numerical scheme is limited to 2D materials of FvK-type behavior and smooth distributions on strains (i.e., it is not capable of describing detailed dislocation stress and pinning stress at the material-material interface).

Numerical Results
For a case where there is no shear between the film and the substrate, and the graphenesubstrate interface is completely slippery, the boundary condition at the far field can be simplified to the film boundary stress being free with no tension.The resulting strain field and PMFs are shown in Figure 5. Figure 5a,b illustrates that the stress distribution induced via bubbling in graphene decreased with increasing distance from the center and exhibited larger strain gradients near the edge of the bubble.Such non-uniform lattice distortion can alter the low-energy electronic band structure of graphene.The magnitude of the electron transition between carbon atoms mimics the impact of an actual magnetic field perpendicular to the graphene plane.As depicted in Figure 5c, the PFMs generated by the axially symmetric graphene bubble exhibited three-fold symmetry, being most pronounced near the bubble's periphery and attenuating rapidly towards its center.This result agrees with previous calculations using simplified theories or different methods, in both distribution and magnitude [105,106].We also found that the maximum magnitude of PFM created this way was proportional to the square of the height-to-radius ratio of the graphene bubble and inversely proportional to the radius of the bubble.In this sense, the potential exists to achieve large PFMs in future studies through preparing bubbles of small radii and large heights using the oxygen plasma-assisted method we report.
for convenience, we may rewrite the Laplacian operator as: where  is an  1  1 matrix, and ⊗ denotes the Kronecker product.As a result, the FvK equations can be rewritten as: where  ,   can be measured from the height image of the bubble.To incorporate boundary conditions, the original matrices are rescaled based on experimental dimensions, resulting in equations satisfying the boundary conditions: and    (20) These equations are then solved through successively taking the inverse of the modified matrices.It should be noted here that this numerical scheme is limited to 2D materials of FvK-type behavior and smooth distributions on strains (i.e., it is not capable of describing detailed dislocation stress and pinning stress at the material-material interface).

Numerical Results
For a case where there is no shear between the film and the substrate, and the graphene-substrate interface is completely slippery, the boundary condition at the far field can be simplified to the film boundary stress being free with no tension.The resulting strain field and PMFs are shown in Figure 5. Figure 5a,b illustrates that the stress distribution induced via bubbling in graphene decreased with increasing distance from the center and exhibited larger strain gradients near the edge of the bubble.Such non-uniform lattice distortion can alter the low-energy electronic band structure of graphene.The magnitude of the electron transition between carbon atoms mimics the impact of an actual magnetic field perpendicular to the graphene plane.As depicted in Figure 5c, the PFMs generated by the axially symmetric graphene bubble exhibited three-fold symmetry, being most pronounced near the bubble's periphery and attenuating rapidly towards its center.This result agrees with previous calculations using simplified theories or different methods, in both distribution and magnitude [105,106].We also found that the maximum magnitude of PFM created this way was proportional to the square of the height-to-radius ratio of the graphene bubble and inversely proportional to the radius of the bubble.In this sense, the potential exists to achieve large PFMs in future studies through preparing bubbles of small radii and large heights using the oxygen plasma-assisted method we report.(a,b) and PMFs (c) induced through graphene bubbles prepared via the oxygen plasma-assisted method under no-shear conditions.The geometry of the bubble in this demonstration is described in Figure 4.
We have discussed a limiting case of zero shear resistance at the graphene-substrate interface.In cases where the substrate underneath, such as silicon oxide, can provide a pinning effect, no-slip conditions would appear more appropriate [39,107].However, a complex boundary layer arises from geometric irregularities of the bubble footprint.Hence, we designated the bubble boundary at a height of 1 nm as its cut-off bottom boundary.As the displacement of the bubble boundary is zero at no-slip conditions, we approximated the boundary condition via applying a biaxial in-plane pulling force at the far field membrane boundary, ensuring zero displacements of the bubble boundary.This approach rendered the governing equations and boundary conditions identical inside the bubble boundary equivalent to the no-slip situation.However, as the boundary was fixed, we disregarded the region outside the bubble boundary.It is important to note that although the geometric equation is nonlinear, the nonlinearity solely appears in the out-of-plane height term, which can be accurately determined through measurement.
Specifically, we began by fitting the contour of the axisymmetric nanobubble using the polynomial function: A constant tension T is applied at infinity to enforce a zero displacement at the boundary r = a.This requires: according to the analytical calculation of the axisymmetric version of Equation (19).Using w(r) to fitted the x and y directions of the bubble respectively, we obtained the corresponding parameters a 1x , a 2x , a 3x and a 1y , a 2y , a 3y .These parameters were then substituted into the Equation ( 22) to derive T x and T y .We approximated that through applying T x and T y , respectively, at the membrane boundary to constrain the displacement of the nanobubbles at the boundary to zero.Subsequently, the strain field within the boundaries was obtained, as shown in Figure 6a,b.
complex boundary layer arises from geometric irregularities of the bubble footprint.Hence, we designated the bubble boundary at a height of 1 nm as its cut-off bottom boundary.As the displacement of the bubble boundary is zero at no-slip conditions, we approximated the boundary condition via applying a biaxial in-plane pulling force at the far field membrane boundary, ensuring zero displacements of the bubble boundary.This approach rendered the governing equations and boundary conditions identical inside the bubble boundary equivalent to the no-slip situation.However, as the boundary was fixed, we disregarded the region outside the bubble boundary.It is important to note that although the geometric equation is nonlinear, the nonlinearity solely appears in the out-ofplane height term, which can be accurately determined through measurement.Specifically, we began by fitting the contour of the axisymmetric nanobubble using the polynomial function: A constant tension  is applied at infinity to enforce a zero displacement at the boundary  .This requires: according to the analytical calculation of the axisymmetric version of Equation (19).Using   to fitted the x and y directions of the bubble respectively, we obtained the corresponding parameters  ,  ,  and  ,  ,  .These parameters were then substituted into the Equation ( 22) to derive  and  .We approximated that through applying  and  , respectively, at the membrane boundary to constrain the displacement of the nanobubbles at the boundary to zero.Subsequently, the strain field within the boundaries was obtained, as shown in Figure 6a,b.Consequently, the solution involving constant tension at the membrane boundaries was equivalent to superposing the solution with no-shear conditions and the solution for a membrane boundary subjected to constant tensions.The strain field of the latter remained constant with a strain gradient of zero, thus not contributing to the PMFs.Furthermore, since we did not consider the strain distribution outside the bubble boundary, it was set to zero.The resulting PMF is shown in Figure 6c.It can be seen that the same Figure 6.Numerical results of strain fields (a,b) and PMFs (c) induced through graphene bubbles prepared via the oxygen plasma-assisted method under no-slip conditions.The geometry of the bubble in this demonstration is described in Figure 4.
Consequently, the solution involving constant tension at the membrane boundaries was equivalent to superposing the solution with no-shear conditions and the solution for a membrane boundary subjected to constant tensions.The strain field of the latter remained constant with a strain gradient of zero, thus not contributing to the PMFs.Furthermore, since we did not consider the strain distribution outside the bubble boundary, it was set to zero.The resulting PMF is shown in Figure 6c.It can be seen that the same bubble under no-slip conditions (Figure 6a,b) can provide large strain magnitudes compared to that under no-shear conditions (Figure 5a,b).Again, we found good agreement in both distribution and magnitude between the results in Figure 6 and previous calculations using simplified theories and clamped boundary conditions [105,106], while the numerical scheme discussed here could be applied to the solution of PMF of bubbles of arbitrary shapes.

Conclusions
In this study, we propose a rapid bubbling method based on chemical reactions between oxygen plasma and organic matter, enabling the rapid induction of in-plane strain and out-of-plane deformation in graphene and other 2D materials.Compared with other bubbling preparation methods, this approach has the advantages of short preparation time, simple operation, controllable bubble shape and position, and good durability.Furthermore, we introduce a numerical method to calculate the strain of the bubble and its induced threefold symmetric pseudomagnetic field under both no-shear and no-slip conditions.This

Figure 1 .
Figure1.Schematic diagram of bubbles prepared via oxygen plasma-assisted method.The substrate with pre-patterned microcavities is pre-processed with oxygen plasma, after which the graphene is micro-mechanically exfoliated onto the treated substrate.Peeling the tape from the heated assembly can lead to the spontaneous bubbling of the transferred graphene sheet immediately at these prepatterned cavities.The upper right panel shows an SEM image capturing graphene bubbles at a particular deflection angle.

Figure 1 .
Figure1.Schematic diagram of bubbles prepared via oxygen plasma-assisted method.The substrate with pre-patterned microcavities is pre-processed with oxygen plasma, after which the graphene is micro-mechanically exfoliated onto the treated substrate.Peeling the tape from the heated assembly can lead to the spontaneous bubbling of the transferred graphene sheet immediately at these prepatterned cavities.The upper right panel shows an SEM image capturing graphene bubbles at a particular deflection angle.

Figure 2 .
Figure 2. AFM height images of bubbles obtained via two preparation methods: (a,b) bubble prepared with a bulge device;(c,d) bubble prepared with the oxygen plasma-assisted method reported in this work.Note that both methods can lead to 1-3 layer graphene bubbles of heights ranging from 0 nm to 200 nm and a diameter of 3.5 µm.

Figure 2 .
Figure 2. AFM height images of bubbles obtained via two preparation methods: (a,b) bubble prepared with a bulge device;(c,d) bubble prepared with the oxygen plasma-assisted method reported in this work.Note that both methods can lead to 1-3 layer graphene bubbles of heights ranging from 0 nm to 200 nm and a diameter of 3.5 µm.

Figure 4 .
Figure 4. (a) Schematic diagram of the bubble model: the topography of bubble (b) before smoothing and (c) after smoothing.Scale bars: 1 µm.

Figure 4 .
Figure 4. (a) Schematic diagram of the bubble model: the topography of bubble (b) before smoothing and (c) after smoothing.Scale bars: 1 µm.

∂x∂y 2 can
be measured from the height image of the bubble.To incorporate boundary conditions, the original matrices are rescaled based on experimental dimensions, resulting in equations satisfying the boundary conditions: Lφ = f and Lχ = φ(20)

Figure 5 .
Figure 5. Numerical results of strain fields (a,b) and PMFs (c) induced through graphene bubbles prepared via the oxygen plasma-assisted method under no-shear conditions.The geometry of the bubble in this demonstration is described in Figure 4.

Figure 6 .
Figure 6.Numerical results of strain fields (a,b) and PMFs (c) induced through graphene bubbles prepared via the oxygen plasma-assisted method under no-slip conditions.The geometry of the bubble in this demonstration is described in Figure4.

Table 1 .
Summary of typical characteristics of bubbles prepared via different methods.