Nanoscale Probing of Interaction in Atomically Thin Layered Materials.

We combine conductive atomic force microscopy (CAFM) and molecular dynamics (MD) simulations to reveal the interaction of atomically thin layered materials (ATLMs) down to nanoscale lateral dimension. The setup also allows quantifying, for the first time, the effect of layer number and electric field on the dielectric constant of ATLMs with few-layer down to monolayer thickness. Our CAFM-assisted electrostatic technique shows that high-quality mono- and bilayer graphene is reliably produced at significant yields only by the shear type of bond breaking between layers, whereas the normal type of bond breaking exhibits a very stochastic process mainly due to the coexistence of local delamination and interlayer twist. Our dielectric constant measurements also reveal a very weak dependence on the layer number and the electric field (up to our experimental limit of 0.1 V/Å), which is in contrast with theoretical reports. Owing to unexpectedly large variations in the screening ability of pristine monolayer graphene under ambient conditions, we further demonstrate that the effective dielectric constant of monolayer graphene can be engineered to provide a broad spectrum of dielectric responses (3.5-17) through oxidation and thermal annealing, thus confirming its much higher chemical reactivity than bilayer and few layers.

Layered materials, such as graphene, MoS2, BN, and many others, are recognized as a distinct class of anisotropic solids with strong chemical bonds within each layer but weak van der Waals (vdW) interaction between the layers. This weak interlayer interaction allows them to be easily sheared parallel and/or expanded normal to the layer surface, leading to the generation of so-called atomically thin layered materials (ATLMs) with few-layer down to monolayer thickness. ATLMs display unique electrical, mechanical, electrochemical, and optical properties that are not essentially observed in their bulk layered counterparts. Therefore, over the past decade, two distinct strategies have been pursued for the synthesis of ATLMs: a top-down and a bottom-up approach. The former generally aims at overcoming the vdW forces between the adjacent layers for the exfoliation of ATLMs from their bulk crystals through mechanical (e.g., Scotch tape exfoliation,[1,2] nanoimprint-assisted shear exfoliation,[3] and electrostatic force assisted exfoliation[4−6]), chemical (e.g., liquid phase exfoliation[7,8]), and electrochemical (e.g., ion/compound intercalations[9]) processes. On the other hand, the bottom-up method depends on the chemical reaction of molecular building blocks to form covalently linked 2D networks by means of catalytic (e.g., chemical vapor deposition, CVD[10,11]), thermal (e.g., epitaxial growth[12,13]), or chemical (organic synthesis[14]) processes. Since ATLM samples with transfer-induced residues and randomly distributed nanoflakes (mainly due to the formation of defects and grain boundaries during the growth process) are inevitable in the bottom-up method, it can be expected that the highest quality samples are still produced by the top-down methods, where individual or combined external normal and lateral shear forces (applied in a direction perpendicular or parallel to the basal plane of the ATLMs, respectively) play a dominant role during the bond breaking process between layers. However, a detailed understanding of the interlayer behavior of the ATLMs under precisely controlled normal and shear loadings is still missing and thus highly desired as an essential step toward enhancing the transfer efficiency and thickness uniformity of ATLMs-based device features and controlling the number of printed flakes onto the substrate more effectively.

We exploit for the first time conductive atomic force microscopy (CAFM) with ultrahigh force–displacement resolution combined with molecular dynamics (MD) calculations to unravel the relative contributions of electrostatic attraction/repulsion, internal layer-to-layer shear, and intermolecular vdW forces to the exfoliation of the ATLMs in general and graphite specifically. Our CAFM measurements together with MD analysis of the electrostatic exfoliation mechanism suggest that anisotropic nature of the vdW interactions in the few-layer graphene (FLG) is the main barrier to the accurate control of the number of printed flakes. In particular, simultaneous monitoring of the interlayer spacing and the relative twist in FLG shows that as the interlayer spacing between the layers increases during the normal exfoliation process, their attractive vdW interaction becomes progressively weaker and weaker, allowing them more freedom to locally delaminate and thus twist relative to one another, followed by the layer sliding and separation. This facile twisting of FLG mainly triggered by the random formation of the local delamination in the normal exfoliation technique results in a random angle of rotation between the adjacent commensurate graphene flakes, making accurate control of the number of printed flakes extremely difficult, if not impossible. On the other hand, the shear exfoliation technique exhibits a very promising and controllable behavior during the printing process by eliminating the interlayer spacing variations and consequently reducing the interlayer twist angles. We also demonstrate for the first time that the relative dielectric constant of FLG is nearly independent of the layer number and the external electric field (up to our experimental limit of 0.1 V/Å), which is in obvious contrast with theoretical models in the literature. Interestingly, our dielectric constant measurements on monolayer graphene show a strong dependence on the surface reaction, which makes it an excellent electric field/charge screening material upon oxidation, and recover its charge storage capability upon thermal annealing. On the other hand, bilayer graphene and FLG exhibit very high oxidation resistance and thus are a better choice for long-term stable electronic devices with higher moisture and oxygen resistance.

Results and Discussion

Figure a presents a schematic illustration of the CAFM-assisted electrostatic manipulation setup, in which an electrically conducting Pt/Ir-coated AFM tip is used in contact mode to perform all measurements. For the scope of this paper, we focus on the exfoliation and characterization of graphene as a model system for other ATLMs. After nanostructure fabrication of 75 nm deep cylindrical mesas with a diameter of 60 nm from highly oriented pyrolytic graphite (HOPG), we utilize an in situ flattened AFM tip to uniformly adhere the selected HOPG mesa to the tip apex with a conductive adhesive polymer PEDOT:PSS(d-sorbitol). The tip with an attached mesa is then brought into contact with the SiO2/Si substrate, followed by applying a bias voltage of up to 10 V between the mesa and the highly doped Si substrate, separated by the 10 nm thick SiO2 film. Pristine graphene monolayers subjected to attractive electrostatic forces are transferred from the mesa onto the SiO2 film as the tip is gently moved away from the SiO2/Si substrate in a direction normal to the basal plane (hereafter referred to simply as normal exfoliation method) or parallel to the basal plane of graphite (referred to as shear exfoliation method). An atomically well-defined contact formed between the mesa and the substrate combined with piconewton force and sub-nanometer displacement resolution in our CAFM setup facilitates the precise evaluation of both the applied force and the vertical/lateral displacement of the mesa with respect to the substrate during the exfoliation process (see Methods and Supporting Information section S1 for more details).

Figure 1

(a) Schematic of the CAFM experimental setup used to perform shear and normal electrostatic exfoliation of FLG from nanosized HOPG mesa onto the SiO2/Si substrate. Insets show the back and bottom view of the tip with an attached HOPG nanopillar. Scale bars indicate 50 and 100 nm, respectively. (b) Shear exfoliation and (c) normal exfoliation histograms of the number of printed flakes collected from 110 and 50 samples, respectively, under different bias voltages. (d) SEM image of mono-, bi-, and trilayer graphene flakes printed by the shear exfoliation method in the form of the letter “M” at V = 10 V. (e) SEM image of monolayer and 15-layer graphene flakes printed by the normal exfoliation method at V = 9.5 V.

In our proposed setup, the exfoliation of FLG features is a combined action of the electrostatic force, applied normal load (being used to improve the conformity of the mesa to the underlying substrate morphology), van der Waals force (due to the substrate–graphene interfacial adhesion and graphene–graphene interlayer cohesion), sliding/retraction speed of the tip, surface properties of SiO2, and ambient conditions. To narrow down the range of possible experimental parameters, we carry out all measurements on the same SiO2 film at zero normal load with a relative tip–substrate speed of 10 nm/s under a clean and controlled environment (20% relative humidity at 21 °C). Our preliminary shear and normal printing measurements in the absence of bias voltage reveal no graphene exfoliation under the aforementioned experimental conditions, allowing us to elucidate the key role of the electrostatic and interlayer vdW forces in the subsequent electrostatic exfoliation process.

The histograms in Figures b and 1c show the number of printed layers as a function of the bias voltage using the shear and normal exfoliation techniques, respectively. Raman spectroscopy coupled with AFM height profile measurements is used to determine the layer number with monolayer accuracy. Ten measurements are taken for each applied bias voltage. It is evident from Figures b and 1c that the shear exfoliation method produces only 1–3 layers (predominantly mono- and bilayer graphene) at different bias voltages (Figure d), whereas the normal exfoliation method yields graphene flakes of various thicknesses (ranging from 1 to ∼20 layers) in a very stochastic manner (Figure e). We will later show in our analysis of the MD trajectories that the weaker interlayer cohesion during the normal exfoliation process facilitates the localized delamination, thereby triggering the relative twist between the adjacent commensurate graphene flakes and thus making accurate control of the number of printed flakes almost inaccessible. In contrast, the shear exfoliation method exhibits much more robust sliding behavior with the slight change in the interlayer twist angles due to enhanced corrugation of the interlayer potential energy. Figure d also shows the SEM image of mono-, bi-, and trilayer graphene flakes printed by the shear exfoliation method in the form of the letter “M” at V = 10 V, further indicating its versatility for the production of graphene flakes with high crystalline quality and uniform thickness (see the inset of Figure d). We note that, regardless of the applied bias voltage, an unexpectedly thick mesa might be produced by the shear exfoliation method provided that any twist grain boundaries exist along the thickness (c-axis) direction of the HOPG nanopillars (Figure S5). AFM measurements[15] and a combination of FIB/SEM and high-resolution TEM[16] also confirm a polycrystalline structure along the c-axis direction of HOPG with a grain thickness of 11–60 nm and 5–30 nm, respectively. Hence, during the attachment of the mesa to the glue-coated tip apex, we moved the AFM tip laterally rather than vertically to achieve a single crystalline HOPG nanopillar, which is necessary to avoid any possible shear exfoliation of thick mesas.

To investigate the atomistic details underlying our experimental results, we first need to correctly understand the role of the number of layers in the dielectric screening properties of FLG flakes. Despite the importance of such a fundamental property for any electronic material, there have been very limited studies with significant diversity in the reported values of the dielectric constant of graphitic systems, ranging from 2 to 16[17−23] (Table S1). Surprisingly, however, there is no direct experimental evidence for the dependence of the dielectric constant of FLG on the layer number and the electric field. To fill this apparent gap, we here report the relative dielectric constant of FLG on the SiO2/Si substrate under different electric fields using dc electrostatic force microscopy (DC EFM). Figure a shows the AFM topography image of 1–8 graphene layers mechanically exfoliated from HOPG onto a 10 nm thick SiO2/Si substrate. A sharp needle is used to gently scratch through the thin SiO2 film and expose the underlying Si for the dielectric constant measurement of the SiO2 film as a validation of our subsequent experimental results. The height profile along the green line in the topography is shown in Figure b. In our setup, the contact potential difference (VCPD), the capacitance gradient (∂C/∂z), and the vdW force (FvdW) between the tip and sample surface are first measured by acquiring the total force (sum of electrostatic force, Fel, and FvdW) on the Pt/Ir-coated tip (SCM-PIT, Nanoworld, with the spring constant of 1.9 N/m) while sweeping the bias voltage between −3 and 3 V on the sample surface at different tip–sample distances. The total force can be given bywhere the first term represents the electrostatic force. This parabolic equation with three fitting parameters (∂C/∂z, VCPD, and FvdW) is used to fit F–V curves, such as the ones shown in Figure c on a 4LG/SiO2/Si sample. It is evident from the offset of the parabolic F–V curves along the y-axis that the contribution of FvdW to the total force is negligible when a bias voltage is applied, in particular, at z > 10 nm. Hence, throughout the experiments described in the following, we scan over the sample from a distance farther than 10 nm to only measure the electrostatic force. The fitting parameter VCPD also reveals a dependence on the tip–sample distance in such a way that VCPD of 4LG varies from 294 mV (at z = 6.1 nm) to 342 mV (at z = 19.7 nm). Using this method, we measured in Figure b VCPD between the tip and sample surface along the same green line in Figure a at z = 10 nm, clearly indicating the layer-dependent surface potentials in FLG up to four layers.

Figure 2

(a) AFM topography image of 1–8LG onto a 10 nm thick SiO2/Si substrate with the corresponding layer numbers labeled. (b) Height profile (blue line) and contact potential difference VCPD profile (red line) corresponding to the green line in panel a. (c) Total force–voltage curves taken on the 4LG/SiO2/Si substrate at each tip–surface distance. Circles are experimental data, and the lines are parabolic fits using eq at a constant lift height. Three fitting parameters ∂C/∂z (aF/nm), VCPD (V), and FvdW(nN) are given for each curve. (d) Measured electrostatic force versus tip–Si distance taken on the bare Si surface (gray circles), on the 10 nm thick SiO2 film (blue triangles), and on the 28LG (red squares) at V = 10 V. The lines are theoretical fittings to eq . Top inset shows that as the tip moves across the sample surface in constant height, the tip experiences a larger electrostatic force on 28LG than that on Si and SiO2. Bottom inset shows the cross-section of 3D finite element calculation of the electrostatic potential distribution between the tip and the 28LG/SiO2 sample (see Figure S7 for the corresponding electric field distribution). (e) Relative dielectric constant as a function of the layer number under relatively low and high bias voltages. The application of the bias voltage ≤3 V makes the dielectric response extremely weak in our setup. The dashed line is a guide to the eyes and represents the dielectric constant of the bulk HOPG. (f) Dependence of the relative dielectric constant of 1–3LG and bulk HOPG on oxygen reaction at V = 10 V.

In order to precisely quantify the relative dielectric constant of FLG, the electrostatic force acting on the tip needs to be calculated by integrating the Maxwell stress tensor over the surface of the probe. Since an accurate analytical model that can exactly reproduce the tip–sample electrostatic interaction is not available, we carry out three-dimensional (3D) finite element electrostatic simulations using COMSOL Multiphysics (AC/DC Electrostatics module) to calculate the Maxwell stress tensor from the electrostatic potential distribution (bottom inset of Figure d) obtained by solving the Poisson equation in a cylindrical space (Supporting Information section S2). We first calibrate the apex geometry of the probe by taking electrostatic force–distance (Fel–z) curves on a conductive surface (e.g., highly doped silicon or HOPG) close to the graphene flakes. However, we note that only local electrostatic force at the tip apex depends strongly on the tip–sample distance within the range 10–150 nm and thus the global electrostatic contribution from the cantilever shank and the cone is negligible. As such, in Figure d the electrostatic force on the tip apex is obtained by subtracting the electrostatic force of the cantilever shank/cone at z > 200 nm from the total force. All Fel–z curves were fitted with our finite-element calculations over a 10–150 nm tip–sample distance at V = 10 V, using the effective apex radius R as the only fitting parameter, while the nominal half cone angle was fixed at θ = 15°. From ten Fel–z measurements (similar to the one shown in Figure d on the bare silicon surface), we found R to be 28 ± 0.5 nm, which is consistent with the nominal value ∼20 nm provided by the manufacturer.

After VCPD and R were determined as a prerequisite for the accurate quantification of the dielectric constant of the FLG, we next measure the Fel–z curves on the graphene flake of different thicknesses, followed by matching the finite-element results to the experimental data using the only fitting parameter ϵr. We illustrate in Figure d the Fel–z curves on the bare silicon (for the sake of tip calibration), on the 10 nm thick SiO2 (for comparison purposes), and on the 28LG (a thicker flake was chosen for more clarity in the figure). From several measurements on different areas of the sample, we obtain ϵr = 20.1 ± 1.9 for the 28LG and ϵr = 3.86 ± 0.67 for the ultrathin SiO2 film (in good agreement with the corresponding bulk material, 3.8,[24] and ultrathin films, ∼4.0 ± 0.9[25,26]). For further comparison, we revisited the dielectric constant of SWCNTs on 2 nm thick SiO2/Si substrate, reported by Lu et al. using a combination of scanning force microscopy and finite element electrostatic simulations.[23] As shown in Figure S8, 3D modeling of an SWCNT of diameter 3 nm as a hollow cylinder rather than a solid cylinder leads to the dielectric constant of ∼22.5, which is more than twice as large as that of a solid SWCNT of the same diameter. This modified value for the dielectric constant of SWCNTs is more consistent with that of 28LG.

We now perform a series of similar dielectric measurements on the 1–8 LG of Figure a under ambient conditions, and the extracted dielectric constants are shown in Figure e. Although the dielectric screening ability of 1LG is relatively weaker (∼20%) than that of bulk HOPG, the overall dielectric response of FLG samples to the external electric field is almost independent of the number of layers. Interestingly, the presence of a relatively strong electric field of E = 0.1 V/Å (or equivalently 10 V/10 nm) does not show any systematic change in the dielectric response of FLG. These observations are in sharp contrast with density functional theory (DFT) calculations of effective dielectric constant of freestanding 2–10LG[27] where ϵr varies from ∼3 (for 2LG) to ∼8 (for 10LG) at E = 0.1 V/Å and becomes electric field dependent for E > 0.01 V/Å.

A relatively large variation in the measured dielectric constant of monolayer graphene under ambient conditions (Figure e) motivates us to study the possible effect of surface reaction on the dielectric response of the FLG. To do so, we oxidized the FLG using a modified Hummer’s method in which the FLG/SiO2/Si substrate was dipped into the diluted oxidant solution (60% H2SO4:0.01 M KMnO4 = 1:1) for up to 5 min, followed by deionized water rinse and an N2 dry. In Figure f, our measurements on the FLG under different exposure times reveal a strong dependence of the dielectric constant of monolayer graphene on the surface reaction which makes it an excellent charge screening material upon 300 s oxidation, whereas bilayer graphene and FLG exhibit very high oxidation resistance. We also observed that vacuum thermal annealing of monolayer graphene at 400 °C for 5 h can fully recover its charge storage capability, making it a unique material with a wide range of dielectric response upon oxidation/thermal annealing.

We next perform classical MD simulations using the LAMMPS simulator[28] to gain atomistic insight into the electrostatic shear/normal exfoliation mechanisms. Eight circular graphene layers with AB stacking and radius ∼2.5 nm are placed at a distance of 3.0 Å above an amorphous SiO2 substrate while the flattened tip is modeled by a tapered silicon (001) layer, as illustrated in Figure a. The 8-LG stack can be printed by displacing the tip upward (to the right) for the normal exfoliation case (shear exfoliation case) with a constant speed. To hold the system in space, 2 Å of the SiO2 substrate from the bottom was treated as rigid throughout the simulation (see Supporting Information section S3 for full details of the atomistic simulation setup). We use our recently proposed spatial discrete model[29] to find the charge distribution within each graphene flake and through the 8-LG thickness and then assign, for the first time, the electric charge of each carbon atom by substituting their position coordinates into the relevant charge density profile (Supporting Information section S4). Figure b illustrates the layer-by-layer charge density profiles in the 8-LG system when a total excess charge density of Q = 1013 cm–2 is induced (see Supporting Information and Figure S12 for the corresponding Fermi level profiles). Our spatial discrete model in Figure S14b also suggests that almost 87%, 91%, and 95% of the total excess charge density reside within the two innermost layers of the 8-LG system upon application of Q = 1012 cm–2 (equivalent to a bias voltage of ∼0.46 V), 1013 cm–2 (∼4.6 V), and 5 × 1013 cm–2 (∼23 V), respectively, implying that the gate-induced electric field can be felt very weakly by layers N > 2. This is consistent very well with our CAFM measurements in Figure e that the relative dielectric constant (which is a measure of charge storage capability and electric field screening in a material) is almost independent of the electric field and the layer number, in particular, for N > 2.

Figure 3

(a) Atomic structure of the 8-LG/SiO2 system. The background color and the arrows in the figure correspond to the local electric fields (color can be read from the scale bar and the length of arrows between layers is proportional to the field intensity). Left inset: density of states in the four innermost graphene flakes versus the electronic band energy, in which the transparent area represents the average induced charge density. Right inset: top view of AB-stacked circular flakes cut out of the rectangular sheet with a mixture of armchair and zigzag edges. (b) Charge density profiles of an 8-LG system for Q = 1013 cm–2, where each dashed line represents the average charge density ⟨q⟩ = Q in layer i. (c) 3D discrete charge density profile of the innermost flake (i = 1) in the 8-LG system for Q = 1013 cm–2 where q1 is the charge density on atom j belonging to the innermost flake.

We next assign the charge of each atom by substituting their radial coordinates into the charge density profile of each layer. In Figure c, we provide 3D discrete charge density profile of the innermost flake in the 8-LG system for Q = 1013 cm–2, indicating the charge variations at the zigzag and armchair edges, as previously confirmed by scanning gate microscope measurements and the charge–dipole model[30,31] (Figure S17). As the last piece in the puzzle of the electrostatic MD simulations, we calculate the attractive electrostatic force of each atom using the well-established concept of the parallel plate capacitor model, which is already verified by experiments for FLG systems[32] and shown to substantially reduce the simulation cost, allowing the detailed study of the problem (see Supporting Information section S5 for full details of the attractive and repulsive force implementations in our MD simulations). After relaxation of the uncharged system at 300 K for 50 ps, we assign the charge of each atom and equilibrate the charged system at 300 K using a Nosé–Hoover thermostat for 10 ps. Then, the attractive electrostatic forces are applied to each atom and the system is again equilibrated for another 10 ps. For the normal (shear) exfoliation process, the tapered silicon (001) layer is pulled in the z direction (x or y direction) with a constant speed of 1 × 10–2 Å/ps (1.5 × 10–2 Å/ps). Newton’s equations of motion are integrated using the velocity Verlet algorithm with a time step of 1 fs. This time step yielded the total energy variation of <0.01% during the whole period of simulations.

In order to provide a quantitative demonstration of the normal and shear electrostatic printing of the FLG, Figure a shows the number of printed flakes as a function of the total induced charge density in the 8-LG. As an illustration, snapshots from the MD simulations of the normal exfoliation for Q = 8.5 × 1012 cm–2 at t = 1 ns and the shear exfoliation for Q = 9.5 × 1012 cm–2 at t = 3 ns are shown in Figures b and 4c, respectively. From Figure a, the minimum induced charge density on 8-LG required for the normal and shear exfoliation of graphene flake is 8.5–9 (1012/cm–2), which is in good agreement with our experimental results for the normal exfoliation (Q ≈ 12.9 × 1012 cm–2 for hs = 10 nm and V = 6 V) and the shear exfoliation (Q ≈ 10.8 × 1012 cm–2 for hs = 10 nm and V = 5 V) and also with the other experimental results for the normal exfoliation of 18 nm wide FLG nanoribbons and 1.4 μm diameter pillars (Q ≈ 3.7 × 1012 cm–2 for hs = 50 nm and V = 8.5 V),[4] the shear exfoliation of 5 μm wide square mesas and 25 μm wide ribbons (Q ≈ 12.4 × 1012 cm–2 for hs = 52 nm and V = 30 V),[5] and also the normal exfoliation of sub-20 nm wide nanoribbons (Q ≈ 8.6 × 1012 cm–2 for hs = 5 nm and V = 2 V)[33] where the total charge density is approximated as Q = ε0εsV/(ehs) according to the parallel plate capacitor model. It is also observed from Figure a that the overall number of printed layers in the shear exfoliation model increases by the increase of the induced charge density, reasonably consistent with our experimental results in Figure b. However, a constant number of printed layers for Q ranging, for instance, from 9.5 to 11 (1012 cm–2) are hypothesized to primarily be the result of the electrostatic screening effect.

Figure 4

(a) Number of printed layers as a function of the induced charge density for both normal and shear exfoliation techniques. Snapshots from MD simulation of (b) the normal exfoliation for Q = 8.5 × 1012 cm–2 at t = 1 ns and (c) the shear exfoliation for Q = 9.5 × 1012 cm–2 at t = 3 ns. (d) A portion of the MD trajectory for the normal exfoliation of the 8-LG system when Q = 10.5 × 1012 cm–2. Variation of the interlayer rotation/distance between layers, labeled 5 and 6, and between 6 and 7 as a function of simulation time. The separation of the layer 7 from 6 is initiated at t ≈ 0.45 ns (highlighted by magenta dashed line). (e) Corresponding snapshot of the MD simulation for such exfoliation taken at t = 0.6 ns. Local delamination is marked in transparent red circles.

Unlike the case of shear exfoliation, it is observed that the number of printed flakes in the normal exfoliation technique does not necessarily increase with the increase of the bias voltage, leading to a random number of printed flakes, as already observed in Figure c. This counterintuitive observation can be understood in terms of anisotropic nature of the vdW interactions in FLG where the interlayer shear strength τs within the basal plane competes with the tensile strength σs (i.e., interfacial cohesion strength) normal to the basal plane during the normal exfoliation/printing course, whereas in the shear exfoliation technique, the interlayer shear strength is primarily responsible for initiating flake sliding and separation. In order to better understand how these physical parameters and their possible interplay can hinder or facilitate the FLG exfoliation, we next establish a quantitative characterization of the interlayer interactions of graphite.

Recent experimental observations on the relative sliding motion of graphite demonstrated that the interlayer shear strength of the AB-stacked (commensurate) graphite flakes (τsc ≈ 140 MPa) is drastically reduced by more than 2 orders of magnitude for their non-AB-stacked (incommensurate) counterparts (τsic ≈ 0.25–2.5 MPa) due to the superlubricity phenomenon in graphite.[34] From experimental measurements[35] and atomistic results,[36] a very slight interlayer rotation (∼2 degrees) between two adjacent commensurate graphene flakes can cause the interlayer shear strength (i.e., interlayer friction) to suddenly decrease by over 50% (Figure S18). This clearly indicates that the interlayer shear strength is very sensitive to the in-plane rotation. In addition, the tensile strength of polycrystalline (incommensurate) graphite normal to the basal plane was measured to be in the range σsic ≈ 10.3–20.7 MPa,[37] which is 1 order of magnitude greater than τsic. To the best of our knowledge, there is no direct experimental measurement of the tensile strength for crystalline (commensurate) graphite σsc. Although a slight difference in the measured interfacial adhesion energy (i.e., basal plane cleavage energy) of the incommensurate (0.37 J m–2) and commensurate graphite (0.39 J m–2)[38] implies that the values of their corresponding out-of-basal-plane elastic modulus (C33) could be relatively close to one another, their tensile strength could exhibit a remarkably different behavior (see Table S2 for more comprehensive data obtained from a wide range of experimental methods and a detailed discussion about the interlayer mechanical properties of FLG/graphite).

Keeping this quantitative description of the vdW interaction of graphite in mind, an evaluation of the MD trajectories and electrostatic interactions indicates that, during the normal exfoliation process, the interlayer shear strength and out-of-basal-plane tensile strength are highly coupled through vdW interactions between the adjacent graphene flakes. Monitoring of the interlayer spacing (Δd = d – d) and the interlayer rotation (Δθ = θ – θ) of the graphene flakes in our simulations (as an illustration, see Figure d for the normal exfoliation of the 8-LG system when Q = 10.5 × 1012 cm–2 and the corresponding snapshot of the MD simulation for such a system in Figure e) reveals that, as the graphene flakes are continuously being expanded during the normal exfoliation process, their attractive vdW interaction becomes progressively weaker and weaker, leading to the facile twisting and sliding of the graphene flakes. This, coupled with our MD observations that the interlayer rotation θ varies within the range −2.5° < θ < 2.5° before the exfoliation is initiated, indicates that adjacent graphene flakes with a larger interlayer rotation are more susceptible to sliding relative to one another under even relatively low shear stress levels. Figure e clearly shows that the interlayer shear stress, mainly induced by the local delamination of the layers during the normal exfoliation/printing course, leads to a complete separation between layers, labeled 7 and 6, rather than, for instance, between 6 and 5 due to the larger interlayer rotation between 7 and 6, as shown in Figure d. Interestingly, our MD results suggest that the normal exfoliation process is always initiated at the edges rather than the middle of the graphene flakes, which can be attributed to the greater electrostatic attractive and repulsive forces caused by the charge accumulation on the edges.

Our analysis of the simulation trajectories further reveals that the shear exfoliation method can effectively suppress the interlayer rotation whose value does not exceed ±0.5° before complete exfoliation is achieved. Given that the interlayer potential corrugation (i.e., the interlayer potential energy variation) is a measure of how easily adjacent layers can slide and rotate relative to one another, our MD calculations for a bilayer system show that the potential corrugation increases when the interlayer spacing is reduced by imposing the attractive electrostatic forces (Figure S19). This finding suggests that, compared to the normal exfoliation, a larger potential corrugation and thus a smaller interlayer rotation in the shear exfoliation technique are caused by the absence of interlayer spacing variations induced by the upward pulling forces.

Conclusions

In summary, we reported the first combined theoretical and nanoscale experimental study on the shear and normal exfoliation of ATLM systems, providing fundamental insights into the accurate control of the number of printed layers. Both experimental observations and MD simulations confirmed that the accurate control of the number of printed flakes is not feasible using the normal exfoliation method. We attributed this result to an intrinsic competition between the interlayer shear strength (which is highly influenced by the interlayer twist angle) and the out-of-plane tensile strength (which strongly depends on the interlayer spacing and local delamination) during the normal exfoliation course. Instead, the ability of the shear exfoliation method to eliminate the interlayer spacing variations and simultaneously suppress the interlayer twist angles (due to the larger interlayer potential corrugation) provides much better control over the desired number of the printed flakes, making it superior to the normal exfoliation method. Our electrostatic force measurements on FLG/SiO2/Si samples suggest a constant relative permittivity nearly independent of the layer number and the external electric field (up to our experimental limit of 0.1 V/Å), which is in sharp contrast with theoretical models. We also demonstrated that the dielectric constant of monolayer graphene can be tuned from 17 to 3.5 upon oxidation and recovered its charge storage capacity by thermal treatment. Notably, bilayer graphene and FLG can retain their chemical inertness under oxidation and thus are well-suited for fabrication of long-term stable electronic devices with higher moisture and oxidation resistance. Our findings about the complex behavior of the vdW interactions between the graphene layers and the way the interlayer shear strength and normal strength change is a general and fundamental result and thus can be used for any other types of exfoliations. While we have specifically focused on the FLG systems, our analyses should be extensible to the electrostatic exfoliation (and mechanical exfoliation as a whole) of other ATLMs, leading to the effective production of 2D materials for the use in high-performance ATLM-based electronic and mechanical devices.

Methods

Sample Preparation

A ∼100 nm thick bilayer of poly(methyl methacrylate) (PMMA) 495K (60 nm)/950K (40 nm) is spin coated onto the freshly cleaved surface of 1 mm thick HOPG substrate (SPI-1 grade with a mosaic spread value of 0.4°), and each layer is baked for 10 min at 120 °C to evaporate the solvent and then patterned by electron beam lithography. After development of the exposed PMMA area in 1:3 MIBK/NMP, a 10 nm thick aluminum film is deposited by thermal evaporation, followed by a lift-off step. To thin down the unprotected HOPG area, oxygen plasma etching is carried out in a reactive ion etching system using pure O2 as the reactive gas. Cylindrical mesas with a radius of 30 nm and etch depth of 75 nm emerge from the HOPG substrate during the plasma etch. After plasma etching, the sample is soaked in 0.1 mol/L KOH water solution for ∼3 min to remove the Al layer, followed by an annealing process at 600 °C under constant Ar/H2 flow for 1 h to remove any resist/metallic residues from the HOPG substrate (Figure S1). Next, the Pt/Ir5-coated tip with the normal spring constant of 2.96 N/m (as measured by the thermal noise method) was scanned in contact mode on the SiO2/Si substrate for 30 min at a load of 2 μN and 70% relative humidity, followed by cleaning the flattened tip via a polishing over ultrasmooth monolayer graphene for 15 min at a load of 200 nN to achieve a residue-free contact surface (Figure S2). Our preliminary SEM observation of the tip at the apex area suggests a very flat triangular shape (Figure S3). Using an approach–retract technique, the flattened tip is coated with a very thin layer of conductive polymer glue by putting the tip apex in gentle contact (at zero normal force) with the prebaked 25 nm thick PEDOT:PSS(d-sorbitol) film on an electrically grounded SiO2/Si substrate. Applying a negative bias voltage of 5 V to the probe results in the formation of raised features in the film, followed by the mass flow of the locally softened polymer toward the tip apex due to localized Joule heating and strong electric field gradient (Figure S4). The location of each mesa is then determined by switching the operational mode of the AFM to noncontact mode, which allows us to avoid any contact between the glue-coated flattened tip and the mesa surface during the image scanning. Although the tip apex is flat, the noncontact mode can still provide us with desired resolution imaging for the subsequent attachment of the mesa to the tip. Switching the mode of operation back to the contact mode, the glue-coated tip apex is then moved to the center of the selected mesa at an applied normal force of 200 nN, and subsequently the mesa/apex contact area is annealed at 95 °C for 30 min using a thin film heater beneath the HOPG substrate. We then move the tip laterally along a single basal plane of graphite, leading to easy shear of the upper section of the mesa (attached to the tip apex) relative to the lower one (fixed to the HOPG substrate), thanks to the extremely low friction of graphite at an incommensurate contact interface (Figure a and Figure S5).

Layer Number Identification

We performed Raman measurements under ambient conditions to identify the number of layers with monolayer accuracy. To avoid laser-induced heating, the laser power at the sample was set to be below 1 mW. Several Raman spectra of each N-LG sample were collected to ensure the repeatability of the results. In our Raman spectra, we did not observe rotating modes R (∼1483–1496 cm–1) and R′ (∼1622–1626 cm–1) nor any change in the integrated intensity of the G and 2D bands in the printed bilayer and multilayer graphene, which is indicative of AB stacking (Figure S6).

Finite Element Simulation

We perform three-dimensional (3D) finite element electrostatic simulations using COMSOL Multiphysics (AC/DC Electrostatics module) to solve the following Poisson equation in a cylindrical space, ∇. (εrε0E) = 0, where εr is the relative permittivity of SiO2, air, or FLG (depending on the subdomain to which the equation is applied); ε0 is the permittivity of vacuum (= 8.854 × 10–12 C2 N–1 m–2); E (= −∇V) is the electric field vector (E, E, E); and V is the electric potential. We also set the following boundary conditions: the electrical potential (e.g., V = 10 V) is defined on the surface of the probe while the bottom surface of the SiO2 substrate is electrically grounded (V = 0). The Neumann condition (dV/dn = 0) is used on the lateral and upper sides of the simulation box (Figure S7). Then, we calculate the electrostatic attractive force exerted on the tip by the integration of the Maxwell stress tensor over the surface of the probe. In a given subdomain, the Maxwell stress tensor (σ) can be expressed by , where I is the identity tensor. We consider a cylindrical simulation box of radius 250 nm and height >300 nm. The probe is modeled as a solid truncated cone of height 250 nm and the half cone angle 15° with a semispherical apex of radius R positioned in the truncated region (tangent to the cone surface, forming continuity in the geometry) and at a distance z from the substrate. The FLG film is modeled as a solid cylinder of radius 250 nm, height h, and relative dielectric constant, εr, over the 10 nm thick SiO2 substrate of εr = 3.8 (which is obtained in the absence of FLG). The air surrounding the probe is modeled as empty space of εr = 1. We adopt an extrafine tetrahedral mesh for the entire simulation box, except for the tip apex and FLG and SiO2 thin films where an extremely fine tetrahedral mesh is applied for better numerical accuracy.

Molecular Dynamics Simulations

In order to provide more accurate atomistic models of the electrostatic exfoliation process, we chose graphene from a wide variety of ATLMs in our MD simulations (and thus in our experiments) because there are very well-established empirical potentials to accurately model few-layer graphene (FLG). We adopt reactive empirical bond order (REBO) potential function[39] to model the intralayer carbon–carbon interactions within the same graphene layer while the free graphene edges are passivated by hydrogen. A registry-dependent (RD) interlayer potential that can accurately describe the overall cohesion, corrugation, equilibrium spacing, and compressibility of FLG is implemented in the LAMMPS code to model the carbon–carbon interaction between graphene flakes.[40] For the same MD simulation but different interlayer potentials (LJ or RD potential), both the number and the orientation of printed flakes were completely different, indicating that the potential corrugation (which cannot be described by the LJ potential) plays a crucial role in determining the intrinsic resistance to interlayer sliding and controlling the exfoliation behavior of the FLG under external electrostatic loads (Supporting Information section S3). Tersoff potential and Stillinger–Weber potential are utilized for the modeling of SiO2 substrate and silicon (001) layer, respectively. Given that the graphene–SiO2 interaction is physisorption in nature, it has been proposed that the short-range vdW interaction is the predominant mechanism at the graphene–SiO2 interface rather than O–C and Si–C covalent bonds.[41−44] As a result, we use a standard 12–6 LJ potential for describing Si–C and O–C interactions according to the Universal Force Field (UFF) model and the Lorentz–Berthelot mixing rules. Although the extreme flexibility of graphene (which makes its interaction with SiO2 more liquid-like than solid-like) and the surface properties of SiO2 play a role in the exfoliation of graphene, the Si–C and O–C interaction parameters, used in this paper, alone do not lead to the graphene exfoliation, allowing us to elucidate the key role of electrostatic and interlayer vdW forces in the exfoliation process. The minimum graphene–SiO2 interfacial adhesion strength required to print monolayer graphene onto the substrate in the absence of electrostatic forces can be obtained for εSi–C = 13.36 meV and εO–C = 5.163 meV which are 1.5 times greater than the interaction energy values we used in this article (see Figure S9). Nevertheless, we will demonstrate later that the electrostatic force can significantly facilitate the print of the graphene flakes onto the substrate with the weak surface adhesion. The glue between the tip and graphene flakes is simply modeled by applying the LJ potential between the silicon layer and the topmost graphene flake using a larger Si–C interaction energy (i.e., εSi–C = 17.8 meV).

Spatial Discrete Model

We now examine the charge distribution of a system with N layers of graphene using our recently proposed spatial discrete model.[29] As schematically shown in Figure a, applying a bias voltage V between the highly doped Si substrate and N-layer graphene (N-LG) induces a total excess charge density of Q in N-LG, whose layer i can carry a charge density of Q such that the following constraint holds: Q = ∑Q. Based on the method of images, the induced excess charge density in a finite-size N-LG stack with a circular shape of radius R can be distributed over the ith layer as where is the charge distribution profile, normalized to its average value ⟨f⟩ for generality purposes; is a dimensionless parameter; r denotes the radial coordinate of atoms; is a polynomial function of which only depends on the ratio of the graphene size to the dielectric thickness and is determined by using the method of images, followed by solving the Love equation; and α (>0) is to determine the amount of charge density at the edge of the ith layer (i.e., ). The charge distribution in the N-LG/SiO2/Si system tends to minimize the total energy Ut which is given as the sum of energy stored in SiO2 as the dielectric medium (Ud), electrostatic energy between the graphene layers (Ue), and the band-filling energy in each layer (Ub). The first two terms can be given by Ud = Q2hs/2ϵ0ϵs and Ue = (dg/2ϵ0ϵg)∑(Q – ∑Q)2, where hs and ϵs are the thickness and dielectric constant of the SiO2 substrate, respectively, dg is the interlayer distance, and ϵg is the dielectric constant in N-LG. Assuming that the electronic band structures remain unchanged under an external electric field, , where D is the density of states (DOS) and is the average value of the Fermi energy profile across the layer i in terms of the constant Fermi energy . In N-LG systems, D is obtained from the summation of DOS for each energy band with double spin and double valley degeneracies as , where Nb (= N/2) is the number of bands in , j = 2l – 1 for even multilayers, ℏ is the reduced Planck constant, vf (= 3γ0a/2ℏ) is the Fermi velocity, γ0 is the nearest neighbor hopping parameter, γ1 is the nearest neighbor interlayer coupling constant, and a = 1.42 is the C–C bond length. Finally, the average charge density of each layer which can be given by is obtained by minimizing the total energy of the system with respect to and α as the variational parameters under the constraint that Q = ∑Q. For our numerical calculations, we take γ0 = 3.14 eV and γ1 = 0.4 eV as typical values of bulk graphite and use the average measured value of ϵg = 20 for the N-LG system.

We next assign the charge of each atom by substituting their radial coordinates into the charge density profile of each layer. To do so, the point charge on atom j belonging to the flake i can be determined by multiplying the corresponding charge density q to the triangular area (Ac = 3√3a2/4) surrounding the associated atom, and can be given by where and with the index notations i and j varying from 1 to N (N being the total number of graphene layers) and 1 to M (M being the total carbon atoms in each layer), respectively; (= r/R) is a dimensionless parameter; r denotes the radial coordinate of atom j in the ith layer which carries the corresponding charge density of q; R is the radial coordinate of the atom at the edge of the layer i; and Q = (1/M)∑∑q is thus the total induced charge density in the FLG. As already shown in Figure b, due to the assumption of infinitely sharp edges, the analytical model fails to predict the charge distribution at the very edge of the flake where different types of graphene edge states (e.g., zigzag and armchair) exist. However, our proposed spatial discrete model can successfully account for the charge distribution at the zigzag and armchair edges of the graphene flake based on position coordinates of each atom.

We finally calculate the attractive electrostatic force of each atom using the well-established concept of the parallel plate capacitor model. As shown in Figures b and 3c, the charge accumulation is confined to a small region close to the graphene edge and thus the electric field can be assumed to be relatively uniform over the majority of the graphene flake area. As a result, the force of each atom which is proportional to the square of its corresponding charge density can reliably be written as F = Acq2/2ϵ0ϵSiO which can successfully capture both the fringe field and screening effects on the attractive electrostatic force which acts on each individual atom in the direction perpendicular to the substrate. Furthermore, the repulsive electrostatic forces due to the like charges on all carbon atoms are computed using the Coulomb pair potential.