Molecular Dynamics Analysis of Graphene Nanoelectromechanical Resonators Based on Vacancy Defects

Due to the limitation of graphene processing technology, the prepared graphene inevitably contains various defects. The defects will have a particular influence on the macroscopic characteristics of the graphene. In this paper, the defect-based graphene nanoresonators are studied. In this study, the resonant properties of graphene were investigated via molecular dynamic simulations. The effect of vacancy defects and hole defects at different positions, numbers, and concentrations on the resonance frequency of graphene nanoribbons was studied. The results indicated that single monatomic vacancy has no effect on graphene resonant frequency, and the concentration of the resonant frequency of graphene decreases almost linearly with the increase of single-atom vacancy concentration. When the vacancy concentration is 5%, the resonance frequency is reduced by 12.77% compared to the perfect graphene. Holes on the graphene cause the resonance frequency to decrease. As the circular hole defect is closer to the center of the graphene nanoribbon, not only does its resonant frequency increase, but the tuning range is also expanded accordingly. Under the external force of 10.715 nN, the resonant frequency of graphene reaches 429.57 GHz when the circular hole is located at the center of the graphene nanoribbon, which is 40 GHz lower than that of single vacancy defect graphene. When the circular hole is close to the fixed end of graphene, the resonant frequency is 379.62 GHz, which is 90 GHz lower than that of single vacancy graphene. When the hole defect is at the center of nanoribbon, the frequency tunable range of graphene reaches 120 GHz. The tunable frequency range of graphene is 100.12 GHz when the hole defect is near the fixed ends of the graphene nanoribbon. This work is of great significance for design and performance optimization of graphene-based nanoelectro-mechanical system (NEMS) resonators.


Introduction
As a novel two-dimensional nanomaterial composed of sp 2 hybrid carbon atoms [1], graphene has excellent properties, such as superior thermal conductivity, high electronic conductivity, and outstanding mechanical properties [2][3][4][5]. Therefore, graphene is considered to be an ideal material for nanoelectro-mechanical system (NEMS) resonators [6]. Graphene-based nanomechanical resonators have attracted a lot of attention [7]. It has high frequency, strong mechanical nonlinearities, high impact factor, small volume and mass, and ultra-high sensitivity [8]. Graphene-based nanomechanical resonators have a broad prospect in many fields, such as radio frequency (RF) communication [9], microwave devices [10], high sensitivity sensors [11], single molecule detectors [12], and macroscopic quantum effect detection [13]. Bunch et al. first fabricated a graphene mechanical resonator with the length and width in micrometer size. The fundamental frequency is on the order of MHz [14]. Jiang et al. studied the free vibration of single-layered graphene-based mass sensor by Galerkin strip DFT method [15]. Numerous experimental and theoretical studies have thus been carried out to utilize graphene in nanomechanical resonators.
However, graphene is not perfect during the production process and structural defects are inevitable [16]. Graphene defects can be divided into three types according to how the defects are created: thermally dynamically resultant, deformation-introduced, and artificially-induced defects [17]. The formation energy of thermally activated defects are generally less, such as point vacancy defects and Stone-Wales defects (SW) [18]. Chemical vapor deposition (CVD) method enables mass production of large-area graphene. However, graphene fabricated by CVD method and epitaxial growth unavoidably contain line defects, such as gain boundary [19]. Adatom defects can be obtained by various chemical and heat treatments [20]. Defects also can be artificially introduced by electron and ion beam irradiation, such as nanopores and holes. Furthermore, the size of nanopores can be controlled by focused ion beam (FIB). Deletion of a carbon atom in the six-membered ring lattice of graphene leads to a single-atom vacancy defect. Figure 1 shows the transmission electron microscope image and atomic structure diagram of a single vacancy defect [21]. The presence of defects breaks the symmetric lattice structure of graphene, which has a significant impact on the mechanical, thermal, and electrical properties. Due to this reason, a lot of experimental and theoretical research has concentrated on investigating the effects of these defects. For instance, Li et al. found that the mechanical properties were significantly reduced with the existence of defect by the molecular dynamic (MD) simulations, and the thermal conductivity of graphene containing SV, DV, and SW decreased 57.6%,~43.4%, and~31% [22]. Ahangari et al. indicated that the point (Stone-Wales (SW) and atom vacancies) and shape defects had a negative effect on Young's modulus using ab-initio-based density functional theory (DFT) calculations [23].
In the literature, abundant experimental and theoretical simulation reports are available on the effects of various defects on the thermal conductivity [24], electron density [25], and mechanical strength [26][27][28] of pristine graphene. The effects of the type, location, and concentration of defects on the resonant properties of graphene have been rarely reported, and the behavior of defective graphene NEMS resonators is not well understood. Thus, it is of great significance to investigate the resonate properties of graphene and effects of various defects. As the device develops towards miniaturization, the mechanical model based on the traditional macroscopic continuum is no longer applicable when the scale enters the nanometer level. Molecular dynamics simulation can be used to study the configuration and energy changes of graphene materials at resonance with atomic accuracy [29]. In this paper, we have conducted a series of MD simulations to explore the resonance law of defective graphene. The effect of vacancy defects and hole defects at different positions, numbers, and concentrations on the resonance frequency of graphene nanoribbons is studied via MD simulation. The results provide necessary theoretical guidance for the design and performance optimization of graphene-based NEMS resonators. The findings can provide a comprehensive understanding in the resonante properties of defective graphene.

Materials and Methods
Molecular dynamics simulation is one of the effective methods to study the characteristics of micro-nano structures [30]. It deeply describes the complex mechanical mechanism of discrete atomic and molecular systems, shows the characteristics that cannot be achieved by experiments and continuum theory, and can accurately carry out numerical calculation.
The molecular models were first constructed. The double-ended fixed-branch graphene nanoribbon established in this simulation is shown in Figure 2, which is armchair-shaped in the length direction. The length l is 10.224 nm and the width w is 3.19736 nm. The molecular model of graphene nanoribbon consisted of 1248 carbon atoms. During the actual experiment, the graphene has to be placed on a SiO 2 substrate, therefore the two ends of the graphene nanoribbon structure are fixed in the simulation model. The red area at both ends is the fixed part, and the fixed atomic number is 104 × 2. The middle blue region is the vibrating part, which contains 1020 carbon atoms. The distance from the graphene resonant beam to the substrate is 5 nm. In this paper, Large-scale Atomic/Molecular Parallel Simulation (LAMMPS) software (Version 24Dec2020, Sandia National Laboratories, Albuquerque, NM, USA) is used for simulation. The time step is 1 × 10 −3 ps (1 fs). The Boltzmann constant in the system is 1.38 × 10 −23 J/K. An Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) potential function is used to simulate the C-C interaction. The AIREBO potential function includes not only the interaction force generated by the potential but also the long-range interaction force. The C-C bond length is 3.4 Å. The cut-off radius is set to 10.2 Å. In the system numerical solution uses the Velocity-Verlet algorithm with high computational efficiency. The boundary condition of setting the nanoribbon X direction in the horizontal plane is periodic boundary condition. The boundary condition of the nanoribbon in the Y direction is periodic boundary condition. An infinite parallel plate capacitor is simulated. A boundary condition is set as a fixed boundary condition in the normal direction of the capacitor, that is, in the Z direction. Nosé-Hoover thermostat was used to control the temperature in the simulation.
The external force applied in this paper is calculated based on the electrostatic force of an infinite parallel plate capacitor. The electrostatic force between the substrate and the resonant beam can be calculated by the following formula without considering the edge effect [31]: Among them, W is the potential energy of the resonator, ε is the dielectric constant, L is the length of the resonant beam, w is the deflection of the resonant beam, V is the voltage between the upper and lower plates, including direct current (DC) bias voltage and small alternating voltage. Derivation of potential energy in displacement direction is the magnitude of electrostatic force between two plates: The minus sign represents the direction of electrostatic attraction between the parallel plates of the capacitor. The external force in molecular dynamics simulation is the electrostatic force calculated by Formula (2), which is uniformly distributed on each atom in the effective region of graphene nanoribbons.

Effect of Single Void Defects on the Resonance Frequency of Graphene
Considering the effects of the location, size, and concentration of vacancy defects on the resonant properties of graphene, single-atom vacancies were set on the graphene nanoribbon. Atoms D1 (51.119, 14.757, 0) at the center position, D2 (51.119, 0, 0) at the edge position, and D3 (17.040, 14.757, 0) at the other positions were removed from the model, respectively. The corresponding structures of graphene nanoribbon models SV-1, SV-2, and SV3 containing single vacancies were obtained, as shown in Figure 3. In the molecular dynamics simulation, the potential function selects the AIREBO potential. The time step is set to 1 fs, and relaxation of 50,000 time steps is performed first. After the relaxation was completed, external forces of 0.297 nN, 1.191 nN, 4.762 nN, and 10.715 nN were loaded, respectively. After that, the energy was minimized again. Then the external force was removed, and the NVT ensemble was changed to the NVE ensemble. At this moment, the kinetic energy of the system and the potential energy between the bonding atoms transform to each other, and graphene nanoribbons vibrate freely. Selecting the kinetic energy or potential energy change of the system as the timedomain signal, the corresponding frequency response curve can be obtained by Fourier transform. The graphene resonant frequency of the different vacancy positions under the action of 10.715 nN external force is shown in Figure 4.
As can be seen in Figure 5, the first-order resonant frequencies of single-atom vacancy nanoribbons at different locations basically overlap. Point M in the figure indicates that the resonant frequency magnitude of graphene is 337.12 GHz under the external force of 0.297 nN, point N indicates that its resonant frequency is 429.9 GHz under the external force of 1.191 nN, point P indicates that the resonant frequency slowly increases to 433.65 GHz when the external force continues to increase to 4.762 nN. Point Q indicates that when the external force is 10.715 nN, the single-vacancy defect graphene at three positions and perfect graphene shows the same resonant frequency, which is 469.53 GHz. The increase in resonant frequency with different applied stresses is consistent with that of the perfect graphene. The resonant frequency of the graphene structure shows an increasing trend with increasing loading forces, which is nonlinear. For the case of consistent frequencies of single-atom vacancy graphene in different arrangements, probably because the loss of a single atom has little effect on the whole graphene nanoribbon system. There is not enough precision to distinguish the minor differences in the statistical physics used in this paper. The loss of a single atom does not change the resonant properties of the graphene nanoribbon in molecular dynamics simulations.  Since the effect of a single atomic vacancy on the resonant frequency of graphene nanoribbons is small, the number of defects is increased to consider the impact of different defect concentrations on the resonant frequency of graphene. In fact, vacancy defects are the most typical defects that can exist stably in graphene and will inevitably be generated in the processing of graphene. For example, when graphene is treated by ion irradiation technology, the number of defects increase with the increase of irradiation dose, and the defects are mainly vacancies. In this paper, the concentration of vacancy defects is expressed as the ratio of missing carbon atoms to the total number of intact graphene atoms. The graphene nanoribbons with vacancy defect rates of 0.294%, 1.4705%, 2.647%, 3.8235%, and 5% were modelled, and the distribution of vacancies on the graphene nanoribbons remained symmetrical and uniform. Figure 6 shows the structure of the vacancy graphene model with different concentrations. Firstly, the graphene nanoribbon model was subjected to a 50,000-step relaxation. The atomic displacement deformation clouds chart during the relaxation process are shown in Figure 7. It can be seen that the carbon atoms near the single-atom vacancy have a larger displacement than carbon atoms in other positions. This is because the carbon atoms near the vacancy have lost a covalent bond between carbon and carbon atoms. Graphene nanoribbons have dangling bonds and are structurally unstable, so these atoms fluctuate more during relaxation.   Figure 8 shows the effect of different vacancy concentrations on the resonant frequency of graphene. It can be seen that the resonant frequency of graphene nanoribbons decreases almost linearly with the increase of vacancy defect concentration. The vibration amplitude of resonant bands also shows a general decreasing trend. When the concentration of vacancy defects is 0.294%, the number of single-atom vacancy is 3. It is the same as the resonant frequency of graphene with one single-atom vacancy defect. When the concentration increases to 5%, the resonant frequency of graphene decreases by 12.77%. In contrast, the concentration of defects in the graphene nanoribbon significantly affects its resonant frequency. In the Euler-Bernoulli beam theory [14], the resonant frequency of the double-ended solidly supported beam when the stress σ acts uniformly on it is given by where E, D, ρ, l, w, and t are Young's modulus, third-order modulus of elasticity, density, length, width, and beam thickness, respectively. Obviously, from Equation (3), the fundamental frequency of graphene nanoribbons decreases with decreasing elastic modulus without changing the size of the model and without strain. The simulation results of Han et al. showed that single atomic and diatomic vacancy defects have essentially no effect on the elastic modulus of graphene films [32]. Therefore, in the simulation of graphene nanoribbons containing only one vacancy defect, the reason why different vacancy positions did not change the resonant frequency may be because a single-atom vacancy did not change the Young's modulus of graphene. For graphene models with different defect concentrations, Zandiatashbar et al. investigated the effect of defects on the mechanical properties of graphene. They showed that the elastic modulus did not change even with a high concentration of hybridized carbon atoms, while vacancy defects significantly reduced the resonant frequency of grapheme [33]. According to Equation (3), the resonance frequency of the nanoribbon is proportional to the Young's modulus. Therefore, when the number of defects increases, that is, the concentration increases, the decrease in the resonant frequency may be caused by the decrease in Young's modulus.

Effect of Circular Hole Defects on the Resonant Frequency of Graphene
Circular hole defects are also one of the common defects in graphene. For example, in the preparation process of reduced graphene oxide there is a strong oxidation reaction, which leads to the loss of atoms and causes defects such as holes. Continuous electron beam irradiation will also cause small holes in the graphene film. Figure 9 shows the spherical aberration correction transmission electron microscope (TEM) image of monolayer graphene with small holes [34].  The parameters used are the same as those used for the simulation of single-atom defects. First, the relaxation process is analyzed, and the two ends of the graphene structure are fixed. Then, the full relaxation is carried out to minimize the energy of the system. The stress cloud chart in the relaxation process is shown in Figure 11. It can be seen that the atoms around the hole are less stressed. This is because there are only a few C-C bonds around these carbon atoms, which are bound by small interatomic forces. The stress near the edge of the nanoribbon is relatively large, which may have a certain impact on the resonant frequency of graphene. The relationship between the resonant frequencies of graphene and different external loads, as well as different circular hole defect locations, is shown in Figure 12. It can be seen that the closer the circular hole position is to the graphene center, the resonant frequency of graphene increases. This increase trend is almost linear. Under the external force of 10.715 nN, the resonant frequency of graphene reaches 429.57 GHz when the circular hole is located at the center of the graphene nanoribbon, i.e., H5 position, which is 40 GHz lower than that of single vacancy defect graphene. When the circular hole is close to the fixed end of graphene, that is, the H1 position, the resonant frequency is 379.62 GHz, which is 90 GHz lower than that of perfect graphene. The change trend of the resonance frequency with the position of the circular hole is consistent under different external force loads. Figure 13 shows the effect of the location of the circular hole defect on the resonant frequency of graphene. It can be seen that as the circular hole gets closer to the center of the graphene nanoribbon, the tunable range of graphene gradually increases, and the rate of increase becomes faster and faster. That is, the variation range of resonant frequency caused by different external forces is larger. For electrostatically driven resonators, a large tuning range can achieve higher sensitivity. Point E indicates that when the center of the hole defect is at the center of the nanobelt, that is, the H5 position, the frequency tunable range of graphene reaches 120 GHz. The tunable frequency range of graphene is 100.12 GHz when the center of the circular hole defect is near the fixed ends of the graphene nanoribbon. As the circular hole defect is closer to the center of the graphene nanoribbon, not only does its resonant frequency increase, but the tuning range is also expanded accordingly. Therefore, the circular hole defect should be located in the center of the graphene nanoribbon.

Conclusions
Studying the effect of defects on nanomechanical resonators based on graphene is essential for obtaining a high-quality graphene resonator. In this paper, we employed molecular dynamic method to investigate the effects of vacancy defects and hole defects on vibration characteristics of graphene-based nanomechanical resonator. First, the vacancy defects of different positions and concentrations are modeled in detail. The LAMMPS software was used to simulate the resonance characteristics of the models. The relaxation performance and the frequency response characteristics of graphene with vacancies under different external loads were emphatically analyzed. The results show that the resonant frequency of graphene nanoribbon increases in a step-like manner as the external force increases. The non-linear deformation dominates the dynamic behavior of graphene. By Fourier transform calculation, we found that single monatomic vacancy has no effect on graphene resonant frequency, and the concentration of the resonant frequency of graphene decreases almost linearly with the increase of single-atom vacancy concentration. In the next stage, nanohole defects were imposed to the graphene nanoribbon. The effect of hole location on the resonant frequency of graphene was analyzed and the simulation results show that holes obviously reduce the resonant frequency of graphene. As the circular hole defect is closer to the center of the graphene nanoribbon, not only does its resonant frequency increase, but the tuning range is also expanded accordingly. The above findings can provide a comprehensive understanding in the resonant properties of defective graphene. Our work provides a necessary theoretical foundation for the following research on the design and performance optimization of graphene nanoresonators.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.