Atomic understanding of structural deformations upon ablation of graphene

We investigate the atomic rearrangement in graphene under femtosecond pulse illumination with reactive molecular dynamics simulations and compare with ultra-fast laser ablation experiments. To model the impact of the laser pulse irra-diation, heat is locally applied to a selected area of the graphene layer and the resulting structural deformation is simulated as a function of time, providing a detailed understanding of the bond breaking process under laser illumination and subsequent re-equilibration after the pulse is turned off. Analysis of the atomic dynamics indicates that the types of defects formed depend on the pulse energy and exposure duration. By varying the exposed area, we determine that the shape of the ablated area is not only a function of the pulse energy, but also of the beam spot size and pulse repetition. Furthermore, we apply a machine learning approach to extrapolate our simulated data to experimental length scales and reproduce the trends in ablated area as a function of temperature. Our study provides a first step towards understanding the design parameters for graphene nano-patterning.

edges with long dangling chains can lead to undesired field enhancement or varying reproducibility of parameters, and which bear the risk of chemical contamination. To gain full control of the shape and size of the patterned graphene, a thorough understanding of the underlying physical mechanisms involved during the ablation process is required.
The process of laser ablation in general involves [22][23][24] (1) light absorption by electrons at the surface layer; (2) transfer of energy from the surface to the target interior by nonlinear electronic heat conduction; (3) energy transfer from the electrons to the lattice, resulting in lattice expansion and formation of metastable electronic states; and (4) structural deformation and formation of cracks. When the laser pulse is ultrafast, these processes are temporally decoupled; therefore, we can study them separately as excitation, energy transfer, and material removal processes. [24] A number of previous studies utilizing analytical approaches have provided an understanding of the energy transfer during the ablation process at a macroscopic level; for example, continuum and hydrodynamic models, [25][26][27][28][29][30][31] finite element simulations, [32] time-dependent Schrödinger equation, [33] ) and twotemperature models. [24,[34][35][36][37] While such studies provide important insight for controlling the ablation rate, to precisely control the shape and size of hole formation, it is beneficial to study and understand the mechanisms on the atomic scale; that is, the breaking and re-forming of the strong carbon-carbon bonds. [38,39] Coupling of first-principles time-dependent density functional theory (TDDFT) with electrodynamics allows for the description of electronic excitations upon ablation, and subsequent transfer of energy to atoms. [22,[40][41][42] However, such an approach is limited in the number of atoms and time-scale studied and so cannot describe large-scale atomic motion. Classical molecular dynamics (MD) simulations can describe these larger scales and have been utilized to elucidate the melting and material removal processes upon ablation of metals and organic solids. [43][44][45][46][47][48][49] However, classical MD fails when it comes to bond-breaking or bond-forming, [48,49] and thus is not accurate for describing covalently bonded materials, for example, as in our case of graphene ablation. In this situation, reactive molecular dynamics (RMD) [50][51][52][53][54] is a promising approach because it accurately describes bondbreaking and formation. RMD has recently been applied to study laser-induced graphene formation, [55] as well as nanogroove formation in graphene upon application of a laser pulse, [56] but has not yet been applied to understanding hole formation via laser pulse.
In this manuscript, we utilize RMD to understand the structural distortions that occur during the process of laser ablation and compare our findings with experimental mea-surements. We apply an intense local heating within a ∼ 1-10 nm radius associated with the laser beam and investigate the deformation of free-standing graphene as the applied temperature (or laser intensity) and heating area (or beam size) are modified. We show the different atomistic distortions and defects that lead to hole formation, and reveal that the shape and extent of the ablated area can be controlled by modification of the laser power and beam size. Lastly, we utilize a machine learning model to extrapolate our simulated results to the experimental micronsized scale, and show good agreement between simulated and measured trends in the ablation area as a function of temperature. This study is a step towards advancing precision graphene patterning by developing an understanding of the atomic-scale processes associated with ablation.

Computational details
All molecular dynamics (MD) simulations were performed using the LAMMPS simulation package [57] with the ReaxFF potential, [50,[58][59][60][61] which has been previously used to study the bond breaking and oxidation of graphene. [54,55,62,63] The MD simulations were carried out at a constant number of atoms and volume with a time step of 0.1 fs. As described in the introduction section, MD simulations are used to study the last step of laser ablation upon receiving the laser pulse and electron excitation. Specifically, MD is used to understand structural deformation and formation of cracks. So to mimic the laser ablation, we used a local Nosé−Hoover thermostat [64,65] to heat up a selected area of a given radius (r), while the atoms outside of this area were kept without a thermostat. Here, the simulated heating profile is rectangular, with an abrupt drop-off whereas we expect a Gaussian drop-off experimentally; this may influence the motion of atoms and ablation shape near the edge of the beam. Because experiments are performed with multiple pulses applied to the sample, we used two periods of thermal ablation each of 100 fs duration with a 10 ps relaxation period in between each pulse. Such time-scales are consistent with the experiment; see Figures S1 and S2. Graphene was placed in a cuboid simulation box containing 34,932 carbon atoms, with periodic boundary conditions applied. Because graphene can be stretched or compressed depending on the experimental conditions, we systematically studied the influence of tensile and compressive forces on the ablation process and structural deformation by changing the length ( ) of the simulation box from 29.2 to 31.2 nm, with equilibrium being 30.1 nm.
The length of the box perpendicular to the graphene plane was fixed at 20.0 nm, allowing the evaporated carbon atoms to sputter away from the surface. To select the atoms in the heated region, we drew a virtual circle with the desired radius ( ) at the center of graphene layer, and coupled the atoms located inside the circle to the thermostat. A range of radius from 0.5 to 9.0 nm, and temperature ( ) from 5000 to 10000 K were considered.
Additionally, we investigated the role of oxygen in the ablation process and determined that, on the time scale of our simulations ( < 0.1 ns), oxygen does not significantly alter the process or the resulting deformed structures, see Figure S3. Hence, in the simulations results presented herein, the oxygen molecules are excluded from the simulations to reduce the computational cost.

Experimental details
For the experimental graphene ablation studies, the graphene sample was irradiated with a 100 fs Ti:sapphire laser at a repetition rate of 80 MHz and a wavelength of 819 nm. Each ablated spot size was illuminated for 15 seconds with varying power levels from 1.0 to 1.8 W, corresponding to fluence values between 19 and 35 mJ cm -2 , resulting in multi-pulse ablation, see ref. [66]. A single layer of commercially purchased chemical vapor deposition graphene was transferred to a silicon wafer with a silicon oxide layer on top and used for the studies. Circular ablated areas, clear lines as well as defect regions were induced through this micromachining process. All experiments were conducted under ambient environmental conditions. We used image processing techniques to determine the boundary of the ablated regions and calculate the ablated area; a detailed description is provided in the Supporting Information.
To estimate the temperature profile generated by the laser pulses, we employ a two-temperature model. The two-temperature model offers one of the most accurate methods to investigate the heat transfer dynamics for femtosecond ablation in metals. [67] For the material removal during the ablation, a critical point phase separation approach is selected. [68] In the simulation, a Gaussian source with the same pulse properties as the experiments is used and its single-pulse interaction with graphene is investigated. For an accurate prediction of the temperature profile, the temperature dependency of the thermophysical parameters of graphene, that is, heat capacity, electron thermal conductivity and the electron-phonon coupling are incorporated based on partial sets of parameters from several studies. [2,[69][70][71][72] Based on this approach, peak electron temperatures associated with different laser fluence values can be derived so that a basic comparison to the MD simulation temperatures can be established.

Dynamics of the ablation process
We first present the dynamics of atoms within graphene in the presence of repeated thermal pulses in Figure 1. A thermal pulse with a temperature = 6000 K and a duration of 0.1 ps is applied, resulting in a rapidly increased local temperature within the targeted area; see Figure 1A. After termination of the pulse, within less than 4 ps, the accumulated heat gradually dissipates through the rest of the graphene layer. All heated atoms, however, do not lose the additional thermal energy at the same rate. As the temperature snapshots in Figure 1A show, at time = 0.5 ps, the outer atoms are cooled within ∼ 0.5 ps, while the inner atoms remain at high temperature until ≳ 2 ps. This delay in the cooling process of internal atoms disrupts the C-C bonds and hence the hexagonal ring structure of those high energy atoms. Thus, the probability of sputtering of carbon atoms from the graphene structure as well as dislocation and defect formation is higher at the central region.
With the application of the second pulse, these distortions are enhanced such that post-equilibrium, there is a clear ablated region at the center of the graphene sheet. Figure 1B presents a representative ablated region, displaying the broken bonds, dangling carbon chains, and deformation of the hexagonal graphene rings, showing examples of disrupted C-C bonds with 1, 2, 3, and 4 connections. As shown in Figure 1C, there is a time-dependence of the formation of these deformations. Within the first ps after the pulse, the number of the pristine graphene's sp 2 carbon atoms, bonded to three atoms ( 3 ), is dramatically reduced while the population of one ( 1 ), two ( 2 ), or four ( 4 ) connected carbons is enhanced. Interestingly, some of the carbon atoms reorganize over time (∼ = 2 ps after application of pulse) such that the number of 3 atoms increases. A similar behavior is observed after interaction with the second pulse. The population of the sp 3 carbon atoms ( 4 ) temporarily increase with the thermal pulse but decay to zero within 3 ps.
The atomic dynamics is temperature-dependent; see Figure S5 for a = 9000 K pulse. While the overall trend of the population dynamics is similar, the fraction of 3 is significantly lowered after the first pulse and does not recover during the cooling process. This suggests a more pronounced deformation at higher temperature such that ablation occurs with fewer pulses.
An increase in the population of 3 atoms after = 2 ps as seen in Figure 1C is not necessarily due to reconstruction F I G U R E 1 Atomistic motions upon ablation. A, Time-dependence of defect formation by two pulses of 6000 K in an area with 3 nm radius. The duration of each pulse is 100 fs followed by a relaxation time of 10 ps. The snapshots show only the target heated area. B, Different bonding structure of carbon atoms within the ablation region: 1 indicates carbon atoms with only one bond, 2 two bonds, 3 three bonds, and 4 four bonds. C, 1 , 2 , 3 , and 4 as a function of time. (d) Snapshots of an ablated graphene layer showing deformations including Stone-Wales and ring defects of the hexagonal structure of pristine graphene. Instead, defects such as the Stone-Wales [73] can be formed in the cooling process, resulting in the formation of pentagonal and heptagonal carbon rings; see Figure 1D. To better understand the process of ring-like defect formation, we plot temporal snapshots of this deformation in Figure 2A.
With application of the laser pulse, the 3 carbon atoms transform into 2 atoms by rotation of a C-C bond. During the cooling period, the initial number of 3 is recovered to some extent but the pristine hexagonal structure is not reconstructed. Other ring-like defects that may be formed are shown in Figure 2B. These defects form because the thermal energy within the ablation region increases the out-of-plane motion of carbon atoms; when the kinetic energy of atoms is high enough to break C-C bonds, but not high enough to liberate any atom, the atoms form extended rings. For the example of Figure 2B, over the course of ∼ 200 fs, a single defect with a radius of ≈ 0.5 nm is formed, resulting in a large number of 2 atoms. We anticipate that the presence of these ring-like defects can provide an experimental signature of laser ablation on the nano-scale since they can be probed via vibrational spectroscopy. [74][75][76][77]

Shape and size characteristics of the equilibrated ablated region
Many factors such as pulse energy, spot size of the laser beam, folding of graphene, or preexisting defects in the  Figure S6 for more snapshots). As shown in Figure 3A, by changing the pulse temperature and the beam size, the size and shape of the ablated area can be modified. To quantify the size of the ablation region formed under these different conditions, we plot the ratio between the ablated ( ) and the initially heated ( 0 ) area, which we label as in Figure 3B-C, = ∕ 0 . Figure 3B plots as a function of the radius of the heated area ( ) with ranging from 0.5 to 9.0 nm. For heated regions with a radius of heating, ≤ 1.0 nm, we observe no bond-breaking until a temperature > 6000 K. As the heating radius is increased, larger cracks form in the graphene even at temperature as low as 6000 K. At higher temperatures, the relative ablation area rapidly increases until reaching a plateau at ∼ 0.8. This indicates that a certain temperature distribution is needed and that at the edges of the pulse, no ablation occurs. Figure 3C presents the same results plotted against temperature, highlighting the temperature-dependence of . For a small beam size ≤ 1.0 nm, grows linearly with temperature, while at high temperatures increases until a plateau is reached. This plateau in terms of beam size for a given temperature suggests that further increase of the beam radius will not increase the size of the ablated region. As discussed below, this plateau of the ablated area at less than 100% of the illuminated area is also observed experimentally.
The shape and size of the laser beam is also highly dependent on internal strain. Although graphene is relatively robust to external stretching and compressive forces due to the strong carbon-carbon bonds, [78][79][80] strain can influence bond-breaking and deformation. This is demonstrated in Figure 3D as a function of temperature. We find that stretching results in a larger ablation area while compression leads to lower area for two chosen beam radii ( = 1.0 and 9.0 nm). This can be understood as being due to the weakened (strengthened) bonding upon application of tensile (compressive) strain; see also Figures S7 and S8. This suggests that for graphene on a substrate, where strain may be present, it is necessary to consider such effects.
In order to identify the shape of the ablated region, we define the fractal dimension index, [81] , as where is the perimeter and is the area of the ablated spot, respectively. For a circle, the approaches 1, while for highly convoluted shapes the trends towards 2. The as a function of temperature and radius is presented in Figures 3E and S6. At higher temperatures or larger beam radii, the ablated area is more circular. However, the effect of temperature on the shape of the ablated area is more pronounced than on the beam radius.
Next, we consider the atomistic deformations formed under different laser conditions. As previously mentioned, the simulations predict deformations of the hexagonal structure around the central ablation region, with interspersed pentagons and heptagons, Figure 1D.
To understand the degree and extent of the deformation, we calculate the relative percentage of pentagon-like  Figure S10 for information about pentagon structure detection. The relative percentage of pentagons within a region that extends over a region of 1.1 × (where is the radius of heated area) is shown in Figure 4A. This percentage increases significantly with temperature and slightly with size of heated area. At a high temperature of 10, 000 K more than 20% of polygons in the selected area are pen-tagons, whereas at 5000 K there is almost none, consistent with little deformation of graphene. If we consider a large region around the ablated area, 1.5 × , the percentage of pentagons is reduced to ∼ 5% ( Figure 4B), indicating that beyond a few atomic shells around the main heated region, the graphene is pristine-like.
Additionally, we analyze the percentage of carbon atoms with three connections ( 3 ) (which occur for the sp 2 bonds of pristine graphene) compared with 1 , 2 , or 4 atoms in Figure 5. The majority of carbon atoms within the heated area are 2 or 3 , that is, having two or three neighbors (see Figure 5B and C). At the lowest temperature considered (5000 K), 100% of carbon atoms are 3 , suggesting little distortion to the C-C bonds (although as previously noted, Stone-Wales defects defects also form 3 carbons). Increasing the temperature enhances the number of 2 atoms and reduces 3 . For temperatures up to 7000 K, the presence of 2 is more or less insensitive to the size of the heated area. At higher temperatures, however, increasing the heating area slightly reduces the presence of 2 and enhances 1 . The formation of 1 indicates the presence of dangling bonds within the ablation region. There are few 4 carbons for any range of radius or temperature,

Multiple pulse enhancement
The majority of our simulations were performed with two pulses of 100 fs duration to simulate laser ablation; however, experimentally more pulses can be applied. We find that the number of applied pulses dramatically changes the ablation area specifically at the lowest temperatures considered. An example is shown in Figure 6 for up to 20 thermal pulses applied with = 6.0 nm at 5000 K (the lowest temperature considered). With < 5 pulses applied, there is no visible sign of ablation, only a few defects are formed. Cracks begin to form after application of 10 pulses and a continuous ablated area emerges after 20 pulses. This finding agrees with experimental results [66] that suggest ablation can be achieved with relatively low energy pulses at a high repetition rate. However, we note that the shape of the ablated area formed by low-energy/high-replication pulses might not necessarily be the same as those created by high-energy/low-replication pulses; we observe a more circularly ablated area for the latter ( Figure 1A).

Long-range topological modification of ablated graphene sheets
Laser pulses not only induce nano and micro cracks in the graphene sheets, but can also introduce new 3D topological structures such as folding. [82,83] This opens new pathways for manipulation of graphene properties at the nano scale without exposing graphene to chemical contamination from etching or e-beam lithography. Our simulations indicate momentary topological waves that dissipate over the course of ∼ 1 ns, and result in a temporary periodically folded layer; see Figure S11 in SI. These waves travel with a speed of ∼ 0.4 − 1.0 cm s −1 depending to the stretching force on the suspended graphene layer. Experimentally, Yoo et al. showed that such a post-ablation topology forms for graphene on substrates but not suspended graphene. [82] Our simulations, however, indicate that such long-range waves can be created by the ablation process but are temporary and do not cause a permanent folded structure on a suspended graphene layer.

3.5
The predictive capability of the simulations In Figure 7A, we present the experimentally determined values as a function of pulse power. For a focused power on the graphene layer varying between 11.8 W, clear ablated circular regions of varying diameter are induced, [66] as F I G U R E 7 A, Experimental normalized ablated area, , against laser power varying from 1.0 to 1.8 W. The optical microscope image inset shows single layer graphene ablated by ultrafast laser pulses for corresponding power levels. B, Comparison of the experimental and simulated ablated area ratio, .The red circle shows an experiment where a modification in the graphene structure is observed instead of a clear ablated spot shown in the optical microscope inset in Figure 7A. In the measurement, as the laser power is increased, the ablation area increases, plateauing at ∼ 0.8. This agrees well with the presented simulations where for the largest heated area ( = 9.0 nm) similar values were derived (see Figure 3B). We note here that the experimental pulse is in the m length-scale while the simulations are on the nm scale. However, as shown in Figure 3B, we predict that the ablated area reaches a plateau well below 10 nm for a given temperature. Thus, we expect that at high temperature, our simulations can be scaled to the experimental dimensions.
To compare the experimentally measured and simulated , we determine the equivalent temperature associated with the applied laser pulse. Different models have been used to evaluate the laser temperature at the ablation spot. [84,85] Here, we estimate the temperature induced by the laser pulse using a two-temperature model as described in the Experimental Details section. Additionally, we extrapolate the simulation to the experimental pulse size. We do this by applying a supervised machine learning (ML) technique; we applied several regressor algorithms, such as decision tree, [86] gradient boosting regressor, and multivariate adaptive regression splines, [87] which are implemented in Python scikit-learn library [88] (version 0.19.1) (see Supporting Information for parameters). We built the models based on 54 simulation data points and considered the simulation results of un-stretched system ( = 0) and excluded data with = 0. To avoid overfitting, we used a two-fold cross-validation technique to get an unbiased estimate of the model performance. [89] We applied the above model (fit to the simulated data) to the experimental conditions and compared with measurement as shown in Figure 7B. As depicted in the figure, ML can predict the trend in a wide temperature range. The estimated values show a similar trend and order of magnitude to experiment, in particular at high temperatures (> 7000 K). The simulation suggests that there are two regimes for ablation; at low temperature (or laser power), the ablation area is quite insignificant with ∼ 0.2 while at higher temperatures the ablated area increases with increasing temperature. This finding is in agreement with our experiments which determine that at low temperature (∼ 4200K), graphene is modified but not ablated; shown by the red circle in Figure 7B. Thus, simulations at the nanoscale can be utilized to approximately predict the ablation radius on the m scale.

SUMMARY
In this study, we investigated the ablation of a single layer graphene under multiple high energy thermal pulses using reactive molecular dynamics simulations. We determined that bond-breaking occurs with application of the pulse and can result in crack and defect formation upon reequilibration. The shape of the resulting ablation region is a function of pulse energy, the beam radius, and the pulse repetition. We classified the frequency and type of defects that occur with these different parameters and by introducing the fractal dimensional index, quantified the shape of the ablated region. Lastly, we used a supervised machine learning approach to extrapolate our prediction to the length scale of measurement and showed that the predicted temperature-dependent trend in ablated area agrees well with the experimental results formed by femtosecond laser beams on a micrometer scale. This work provides a thorough understanding of the atomic-scale processes associated with laser ablation and is thus a first step towards designing nano-patterned graphene.