# Adhesion Behavior between Multilayer Graphene and Semiconductor Substrates
## 2 Method
### Description of the Potentials and Running Simulation  
In the simulation, to model the interaction of carbon atoms in graphene, Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) potential [32] was adopted. It has been shown that the AIREBO potential is capable of calculating the interlayer interaction of graphene very well due to  
Figure 1: Molecular dynamics (MD) model of three-layer graphene on an Si substrate (front view and oblique view). Green atoms represent the topmost layer of graphene, violet atoms represent the second layer of graphene, blue atoms represent the bottom layer of graphene; the silicon surface is an Si (100) crystal plane.  
the inclusion of the long-range Lennard-Jones (LJ) potential [32]. The AIREBO potential has been widely used for graphene in MD simulations [19; 21; 33]. It is in a formation similar to a pairwise dispersion-repulsion potential.  
\[E=\frac{1}{2}\sum_{i}\sum_{j\
eq i}\left[E_{ij}^{REBO}+E_{ij}^{LJ}+\sum_{k\
eq i }\sum_{l\
eq i,j,k}E_{kijl}^{TORSION}\right] \tag{1}\]  
\(E_{ij}^{REBO}\) is the REBO potential [34], which describes the interactions between intralayer graphene atoms, and its formulation is:  
\[E_{ij}^{REBO}=\sum_{i}\sum_{j(>i)}\left[V_{ij}^{R}\left(r_{ij}\right)-b_{ij}V _{ij}^{A}\left(r_{ij}\right)\right] \tag{2}\]  
\(V^{R}\left(r_{ij}\right)\) is the function of all interatomic repulsive interactions while \(V^{A}\left(r_{ij}\right)\) represents the attractive interactions. The distance between pairs of nearest-neighbor atoms, \(i\) and \(j\), is expressed by \(r_{ij}\). The bond-order function, \(b_{ij}\), can deal with many-body effects and is capable of calculating the formation and cracking of covalent bonds. This is because the bond-order function, \(b_{ij}\), not only contains a dihedral angle interaction weighing the bond strength, but also describes radicals and the conjugate or non-conjugate structure of carbon, as well as the rotation barrier, which prevents unrealistic bond rotation from occurring. The long-range interactions between graphene atoms are determined by the LJ potential, which is expressed as \(E_{ij}^{LJ}\). \(E_{kijl}^{TORSION}\) depends on the four-body potential torsion term of the dihedral angle.  
The Tersoff potential was applied to calculate the interactions among Si atoms [35]. This potential has been widely used for modeling Si [36; 37]. The interaction between graphene and the Si substrate was described by the van der Waals interaction using the LJ potential function, and the interaction strength can be modified by simply changing the energy depth [38]. The validity of these potential functions used in our model was verified by the fact that the simulated physical parameters, such as the lattice constant, Young's modulus, and shear modulus, etc., were very near to the experimental results [32; 37; 38; 39].  
Our MD simulation was performed with the LAMMPS code [40], which is widely used. The system was consistent with an NVE microcanonical ensemble. Within this ensemble, the number of atoms, \(N\), the system volume, \(V\), and the energy, \(E\), were conserved to ensure that the total energy was constant. The equations of particle motion were solved using the Verlet algorithm, and the simulation timestep was 0.001 fs, which was adequate for system relaxation by examining the stability through the root mean square deviations of the atoms. All the simulations were performed at 300 K. A constant vertical upward velocity was applied to graphene during the exfoliation process. In addition, we systematically studied the whole range of velocities: First, for each model where the number of graphene layers was one to ten, the calculations were executed 20 times with velocities ranging from 1 A/ps to 20 A/ps with an interval of 1 A/ps; second, according to different exfoliation phenomena from the simulated results obtained in the first step, we tried to narrow down the range of velocities to determine the critical V\({}_{\text{up}}\) with an interval of 0.01 A/ps.