Interaction Energy between Graphene and a Silicon Substrate Using Pairwise Summation of the Lennard-Jones Potential

We consider the computation of a long-range interaction energy between a single graphene sheet and a silicon substrate, which arises from vacuum fluctuations. The interaction energy obtained by summation of the LennardJones potential between a carbon atom in a single graphene sheet and a silicon atom is compared with the dispersion energy (Casimir energy) obtained by combining the Lifshitz theory and the Dirac model for graphene. Deviation of the pairwise summation of the Lennard-Jones potential from the Casimir energy is corrected by adding a powerfunction term, whose coefficient depends on the distance between atoms. We also consider the interaction between a graphene sheet and a silicon dioxide substrate. [DOI: 10.1380/ejssnt.2017.40]


I. INTRODUCTION
The excellent mechanical properties of graphene [1,2] have led to many promising applications in nanoelectromechanical systems (NEMS), such as a nanoelectromechanical contact switch [3] and resonator [4][5][6].These new devices may open a novel way to reduce energy consumption in information appliance [7] and realize highly sensitive sensors [8].To understand the mechanical structure and dynamical properties of graphene-based NEMS, molecular dynamic (MD) simulations are often used.The accuracy of MD simulations mainly depends on the degree of accuracy of the force between molecules used in the MD simulations.In particular, since the adhesion between graphene and a substrate is a key issue in the graphenebased NEMS [9], the modelling of the interaction between a graphene sheet and substrate is very important.
The Lennard-Jones (6,12) potential is frequently used to express the interaction energy between a graphene sheet and substrate [10][11][12], which originates from the fluctuation of an electromagnetic field in a vacuum.The greatest advantage of the Lennard-Jones potential over other interatomic potentials may be the simplicity of the expression, and this strong point is helpful to reduce the computation time.However, the accuracy of the Lennard-Jones potential usually decreases as the distance between atoms increases because of the finiteness of the speed of light.For large separations, a retardationation effect must be taken into account.In addition, mutlibody effects must be taken into account.Recently, the interaction between a graphene sheet and a flat substrate using the Lifshitz theory has been studied [13][14][15][16][17], and a good agreement between the theoretical value calculated by the combination of the Lifshitz theory and Dirac model and the measurement using a dynamical atomic force microscope (AFM) was reported [18].
In order to consider the modelling of the graphenebased NEMS using the Lennard-Jones potential, we examine the difference between the interaction energy obtained by a pairwise summation of the Lennard-Jones potential [19,20] and that obtained by the Lifshitz theory (Casimir energy) [21,22].The multi-scale simulation is one possible scheme to correctly simulate graphenebased NEMS.The interaction between carbon atoms in a graphene sheet is described by atomic interaction energies such as the Tersoff potential, and the local pressure between a graphene sheet and substrate is calculated by the Lifshitz theory.However, in this study, we consider another approach to improve the MD simulation.We generalize the Lennard-Jones potential to agree with the Casimir energy between a graphene sheet and a substrate arranged in parallel without significantly losing the simplicity of the Lennard-Jones potential.
This paper is structured as follows.In Sec.II, we explain briefly the Lennard-Jones potential.In Sec.III, the interaction energy between a graphene sheet and silicon substrate is calculated using the Lennard-Jones potential, and the asymptotic behavior for large separations is considered by employing the continuum approximation.In Sec.IV, to calculate the interaction energy between dielectric plates, the Lifshitz theory is introduced.In Sec.V, the interaction energy between two parallel silicon plates is calculated using the Lifshitz theory and compared with that calculated using Lennard-Jones potential.To characterize the difference between the results obtained by these methods, we introduce the local power exponent that represents the degree to which the interaction energy decreases by increasing the separation distance.In Sec.VI, the local power exponent of the interaction energy between a graphene sheet and a silicon substrate is calculated, and the generalization of the Lennard-Jones potential is considered.In Sec.VII, the representation of the interaction energy between a graphene sheet and silicon dioxide using the generalized Lennard-Jones interaction is described.In Sec.VIII, we summarize our results and present the method of improving the MD simulation for the graphene-based NEMS.

II. LENNARD-JONES POTENTIAL
To express the interaction energy between a pair of neutral atoms, we often employ the Lennard-Jones (6-12) potential: e-Journal of Surface Science and Nanotechnology  where r is a distance between atoms.The first term represents repulsive interactions for small separations, and the second term represents the van der Waals interaction.The Lennard-Jones parameters ϵ and σ depend on the combination of elements.For the combination of identical atoms, Rappé et al. [23] presented a table of Lennard-Jones parameters for various elements.The Lennard-Jones parameters between different elements A and B can be approximately determined by the following relationship where σ X−X and ϵ X−X are the Lennard-Jones parameters between a pair of element X.The Lennard-Jones parameters used in this study are given in Table I, where the parameters are calculated based on molecular mechanics force field given by Rappé et al. [23].We note that V (r) takes the minimum −ϵ at r = 2 1/6 σ.

III. PAIRWISE SUMMATION OF THE LENNARD-JONES POTENTIAL
We consider the interaction energy between suspended graphene and a silicon substrate as shown in Fig. 1 using the Lennard-Jones potential.We assume that the suspended graphene is located parallel to the substrate with a separation gap a.The configuration of carbon and silicon atoms is shown in Fig. 2. Silicon atoms are arranged in a diamond cubic crystal structure with a lattice constant d Si of 5.431 Å.The origin of our chosen coordinate system is located at a center of a six-membered ring, where the z-axis is parallel to the (0,0,1) axis of silicon.We set the x and y-axes to be perpendicular to the z-axis.Carbon atoms are arranged in a honeycomb structure, and the minimum distance between carbon d C is 1.421 Å.In general, the interaction energy between two bodies A and B is given by the summation of the Lennard-Jones potential over all possible combinations of atoms in two bodies.We assign positive integers i and j to atoms of bodies A and B, respectively.Then, the interaction energy is given by where r ij is the distance between atoms labelled by i and j.
Let us consider the interaction energy between a single carbon atom that is located at (x, y, a) and a substrate.If the silicon substrate is semi-infinite, the interaction energy depends only on a and the relative position of a carbon atom measured from the nearest silicon atom on the top surface of the substrate.When the position of the nearest atom is (X, Y, 0), the relative position of the carbon atom on the x − y plane is determined by (ξ, η) ≡ (x − X, y − Y ), where −d Si /2 ≤ ξ, η ≤ d Si /2.Thus, the interaction energy is written as U (ξ, η, a).The number density function of carbon atoms observed at (ξ, η) is uniform due to the lattice mismatch between graphene and the silicon substrate.Thus, the interaction energy Volume 15 (2017) Inui and Iwasaki between a graphene sheet and a silicon substrate can be expressed using the interaction energy between a carbon atom and a silicon substrate U (ξ, η, a) as U (ξ, η, a)dξdη. (5) In the practical computation of U (ξ, η, a), the size of the substrate is limited.We introduce l max as a parameter that characterize the size of the substrate and sum over all atoms in the substrate within a domain −l max /2 ≤ x, y ≤ l max /2 and −l max ≤ z ≤ 0 in the calculation below.
If the separation distance between the graphene sheet and silicon substrate is much larger than the lattice constant of the crystal silicon, the distance between the carbon and silicon atoms changes almost continuously from a to ∞.Thus, we may use the continuum approximation, and the potential between the graphene sheet and silicon substrate is given by Here, r 1 ≡ {x 1 , y 1 , z 1 } denotes the position in the silicon substrate, and the integral for r 1 is calculated over the domain −∞ ≤ x 1 , y 1 ≤ ∞ and z 1 ≤ 0. The position in the graphene sheet is denoted by r 2 ≡ {x 2 , y 2 , a}, and the integral for r 2 is calculated over the area The distance between the carbon and silicon atoms R(r 1 , r 2 ) is expressed by We denote the number density of silicon and carbon atoms as ρ Si = 4.994× 10 28 /m 3 and ρ C = 5.246×10 20 /m 2 , respectively.By employing cylindrical coordinates, the integrals in ( 5) are analytically calculated, and the interaction energy is given by For large separations, the interaction energy can be asymptotically expressed by Figure 3 shows |U (a)| calculated by ( 5) on a log-log scale, where the integral is evaluated with a grid interval length of 0.01d Si for l max /d Si = 100, 200, and 300.The solid line shows |U a≫σ (a)|, whose slope is −3.We confirm that U (a) approaches U a≫σ (a) as a increases.However, a low convergence is also observed.For a large separation, a large system size l max is necessary.The difference between U (a) and U a≫σ (a) for small separations is due to a lack of validity of the continuum approximation.

IV. LIFSHITZ THEORY
In graphene-based NEMS such as nanoelectromechanical contact switches, the separation distance between the graphene sheet and substrate is much larger than the length σ characterizing the Lennard-Jones potential.If expression (1) was valid for an arbitrary distance a, and the interaction energy between the graphene sheet and substrate was given by the summation of the Lennard-Jones potential between a carbon atom and an atom in the substrate, U a≫σ (a) gives a good approximate potential.However, in actuality, these condition are not satisfied.The retardation and many-body effects must be taken into account.The retardation effect between dielectric bodies was studied by Casimir [21] and generalized by Lifshitz [13,22].According to Lifshitz theory, the interaction energy between dielectric plates at temperature T , which is the Casimir energy, is given by Here, the normalized Matsubara frequency ζ l is defined by 4πk B T al/cℏ where c is the speed of light in vacuum, k B is the Boltzmann constant, and ℏ is the reduced Planck constant.The variable y is defined by where k ⊥ is the modulus of the wave-vector projection of light on the plate.The summation is taken over two contributions of transverse magnetic (TM) and transverse electric (TE) modes, and the prime on the summation symbol denotes that 1/2 should be inserted if l = 0.The function r (n) σ (iζ l , y) in (10) is the reflection coefficient of the n-th plate defined by where ϵ (n) ≡ ϵ (n) (icζ l /(2a)) is the permittivity of the nth plate along the imaginary axis.Using the Kramers-Kronig relation [22], the permittivity along the imaginary axis is calculated from the imaginary component of the permittivity as The Casimir force per unit area P C (a, T ) ≡ −∂U C (a, T )/∂a is given by If the permittivities are one for arbitrary frequencies, which they are for perfectly conductive plates, the Casimir energy at absolute zero is given by

V. CASIMIR ENERGY BETWEEN SILICON PLATES
The electric properties of graphene are very unique, and the expression of the reflection coefficient is complicated in comparison with other materials.Thus, we will begin by considering the Casimir energy between silicon plates.The material dependence of the Casimir energy is mainly determined through the permittivity along the imaginary frequency axis.We obtained the permittivity of silicon along the imaginary frequency axis by substituting the imaginary component of the permittivity, which is taken from the optical data of silicon [24], into Kramers-Kronig relation (13).
Figures 4(a) and (b) show the Casimir energy and force, respectively between silicon plates at 300 K on a log-log scale.For small separation distances, the Casimir force is approximately expressed by a potential energy  This formula is rewritten as where H is the Hamaker constant [20].Similarly, the Casimir energy for small separation distances is given by The Hamaker constant for silicon plates is 2.634 × 10 −19 J.
Let us estimate the Hamaker constant for silicon plates using the Lennard-Jones potential.Using the continuum approximation, the interaction energy is given by Thus, the term representing the van der Waals interaction is By comparing ( 20) and ( 22) both methods successfully predict that the interaction decays proportional to a −2 for small separation distances.Equation ( 22) yields the Hamaker constant for the interaction energy: For silicon plates, the Hamaker constant obtained by the summation of Lennard-Jones interaction is 2.747 × 10 −19 J, and is slightly different from that obtained by the Lifshitz theory.In addition, it is found that the slope in Fig. 4(a) increases as the separation distance increases.To observe this change, we assume that the Casimir energy can be expressed in the following form: Although the local power exponent α and coefficient C(a) depend not only on the separation distance but also on the temperature, we omit the temperature dependence for simplicity and concentrate our attention on the dependence of the separation distance.The exponent α is determined from the Casimir energy and force by Figure 5 shows the dependence of the local power exponent α on the separation distance at 300 K. We note that the thermal radiation and permittivity of silicon substrate depend on temperature [25].Thus, α also depends on temperature.As expected, α approaches 2 as the separation distance decreases.Conversely, as the separation distance increases, α approaches 3 and decreases as the distance decreases further.This change in the power exponent is mainly due to the finiteness of the speed of light.As seen in (15), the Casimir energy between perfectly conductive plates decreases proportional to a −3 as the separation distance increases.This implies that the asymptotic behavior of the Casimir energy between silicon plates is similar to that between perfectly conductive plates at large separations.
To examine the deviation between the interaction energy between silicon plates by summing the Lennard-Jones potential from that obtained by the Lifshitz theory, we introduce the relative error defined by Figure 6 shows the relative error of the potential energies between obtained based on the Lifshitz formula and the pairwise summation method.Since the Casimir energy depends on the temperature, the relative error also depends on temperature.Fig. 6 shows the relative error at 300 K.The deviation increases with the separation distance.
Here, we assume that the interaction energy between two atoms is expressed by the following power function where C β > 0. For example, the van der Waals interaction can be expressed as a special case of β = 6.We consider http://www.sssj.org/ejssnt(J-Stage: http://www.jstage.jst.go.jp/browse/ejssnt/) e-Journal of Surface Science and Nanotechnology Volume 15 (2017) the interaction energy between semi-infinite plates, whose atoms interact with the potential in (27).If β > 4, the interaction energy between the semi-infinite plates obtained by the continue approximation is given by For β = 6 and 7, we have The local exponent of U β (r) is β − 4, and that for silicon plates changes between 2 and 3. Accordingly, the local exponent β in (27) changes between 6 and 7.Because the power function with an exponent of 7 is necessary to express the interaction energy for large distances, we generalize the term representing the attractive interaction in the Lennard-Jones potential between silicon plates separated with a gap a in the following form: Using the generalized Lennard-Jones potential, the attractive interaction energy between two plates including atoms with number densities ρ 1 and ρ 2 is given by .

(32)
As our focus is on the interaction energy for a ≫ σ, the repulsive term in the Lennard-Jones potential, which rapidly decreases as the distance between atoms increases can be neglected.The coefficients C6 (a) and C7 (a) are determined by the solution of the following simultaneous equations: Figure 7 shows the coefficients C6 and C7 for silicon plates, which are normalized by ϵ Si .For small separation distances, C7 is very small.This means that (1) is a good approximate expression for small separation distances.On the other hand, C7 increases in comparison with C6 .This means that the retardation effect becomes more important as the separation distance increases.

VI. CASIMIR ENERGY BETWEEN A GRAPHENE SHEET AND SILICON SUBSTRATE
We consider the interaction between suspended graphene and a silicon substrate in the framework of the Lifshitz theory.The main task involving the calculation of the Casimir energy given by ( 10) is to express the reflection coefficient as a function of the frequency and the modulus of the wave-vector projection of light on the plate.The electric properties of graphene are quite different from ordinary materials, and quantum field theory is required to derive the reflection coefficient of graphene.The one of the important parameters characterizing the reflection coefficient of a two-dimensional material is a mass gap of quasiparticle excitations ∆.This value for graphene is conjectured to be very small, and thus we assume that ∆ = 0.In addition, the temperature dependence is expected to be a small except term with l = 0. Thus, we used the refection coefficient at absolute zero except for the term with l = 0.
If the interaction energy between a silicon and carbon atom in a graphene sheet is expressed by (27), the interaction energy between a graphene sheet and silicon substrate obtained by the continuum approximation is given by By comparing with (28), we find the exponent of a is decreased by one.This is because graphene is twodimensional material.Figure 8 shows that the local exponent α(a), which is defined in the previous section, changes between 2 and 3, which is similar to Fig. 5.The remarkable difference between Fig. 5 and Fig. 8 is that the local exponent for the combination of a graphene sheet and silicon substrate is considerably larger than 2 at small separations, and it differs from the asymptotic behavior in (20).Figure 9 shows the relative error, which is defined similarly with (26).The relative error is small in comparison with that for the silicon plates.The reason is explained below.The exponent obtained by the continuum approximation in (34) is β − 3. Thus, if we assume the form in (27) for the interaction between a carbon and silicon atom, the exponent β must change between 5 and 6.Thus, we define the generalized Lennard-Jones potential between a carbon atom in graphene and a silicon atom by Volume 15 (2017) This artificial interatomic potential yields the interaction energy between a graphene sheet and silicon substrate: 3 Figure 10 shows the changes in C5 (a) and C6 (a) normalized by ϵ Si The coefficient C5 (a) decreases as the separation increases.Interestingly, the extinction of C5 (a) at large separations implies that the interaction energy between a graphene sheet and silicon substrate can be approximately expressed by modifying only ϵ Si-C in the Lennard-Jones potential (6,12).For instance, C6 (a) = 2.0 ϵ Si-C and C5 (a) = 3.8 ×10 −3 ϵ Si-C at a = 100 nm.Thus, the interaction energy at a = 100 nm and T = 300 K can be expressed as a summation over the Lennard-Jones (6,12) potential with 2ϵ Si-C .

VII. CASIMIR ENERGY BETWEEN A GRAPHENE SHEET AND SILICON DIOXIDE SUBSTRATE
Suspended graphene is often fabricated above a silicon dioxide substrate.Thus, we consider the Casimir energy between a graphene sheet and silicon dioxide substrate.To calculate the Casimir energy based on the Lifshitz theory, the permittivity of silicon dioxide is necessary.The structure of a silicon dioxide layer on a silicon wafer produced by thermal oxidation is amorphous.Thus, we use the permittivity of fused quartz described by where C UV = 1.098,C IR = 1.703, ω UV = 2.033×10 16 rad/s, and ω IR = 1.88×10 14 rad/s [22].Figure 11 shows the local power exponent of the interaction energy between a graphene sheet and SiO 2 substrate.When calculating the interaction energy between the graphene sheet and SiO 2 substrate as a pairwise summation of the Lennard-Jones potential, two interatomic interactions must be considered.The first is the interaction between a carbon atom in the graphene sheet and a http://www.sssj.org/ejssnt( silicon atom.The second is the interaction between a carbon and an oxygen atom.Thus, we express the general-ized Lennard-Jones potential (attractive terms) for these interactions as Unlike the generalized Lennard-Jones potential between the graphene sheet and silicon substrate, there are four unknown parameters.These parameters can be determined by introducing higher derivatives of the Casimir energy.However, as an alternate method, we simply introduce the following constraint for β = 5 and 6: By this constraint, the number of independent parameters decreases to 2. The interaction energy between the graphene sheet and a silicon dioxide substrate obtained by the continuum approximation is given by ŨSiO where ρ Si = 2.641 ×10 28 /m 3 is the number density of silicon atoms.Here we note that the factor of two before ϵ O-C /ϵ Si-C (= 0.386) originates in the assumption that the number density of oxygen atoms is two times larger than that of silicon atoms.We calculated the interaction energy between a graphene sheet and a flat substrate by two methods and examined the difference between them.In the first method, the interaction energy is calculated by summation of the Lennard-Jones potential over all possible combinations of carbon atoms in a graphene sheet and atoms in the substrate.In the second method, both the graphene sheet and substrate are regarded as dielectric bodies, and the Lifshitz formula is applied.For large separation distances between the graphene sheet and substrate, the second method yields a more accurate interaction energy.However, the first method is more simple and suitable for MD simulations.
To clarify the difference in the interaction energies obtained by the two methods, we introduced the local power exponent.The interaction energy obtained by the first method with continuum approximation decays, obeying a power law, as the separation distance increases, and its power exponent is independent of the separation distance.On the other hand, the interaction energy obtained by the second method cannot be expressed by a power func-tion with a constant power exponent.Thus, we expressed the interaction energy by a power function whose exponent depends on the separation distance.The difference between the two methods is mainly caused by the presence or absence of the retardation effect.Accordingly, the change in the local power exponent on the separation distance provides the degree of the retardation effect.
We presented the generalized formula of the Lennard-Jones potential to accurately express the interaction energy between a graphene sheet and silicon substrate.The Lennard-Jones potential includes a term describing the van der Waals interaction, which is proportional to r −6 .By summing this term, the van der Waals interaction between semi-infinite plates is correctly described for small separation distances within the continuum approximation for combinations of semi-infinite slabs.By adding the retardation term r −7 , the coverage of the generalized Lennard-Jones potential can be extended over a range of separation distance.Using the pairwise summation method is that the interaction energy between curved surfaces.On the other hand, the Lifshitz formula in (10) can be used only for parallel surface.
The interaction energy between a two-dimensional material and semi-infinite substrate, which are arranged in parallel cannot be described by summation of the r −6 term for small separations.A term proportional to r −5 is required.This additional term originates in a lack of integration over the axis perpendicular to a surface of a graphene.We show that the linear summation of r −5 and r −6 can be used to express the interaction energy between a graphene sheet and silicon substrate.This does not mean that the new term is needed to express the interaction energy between a carbon and silicon atom for short distances between atoms.The generalized Lennard-Jones potential presented here is an effective potential to express the macroscopic bodies.The procedure of calculating the interaction energy between graphene and a substrate is summarized as follows.First, interaction energy is calculated using Lifshitz theory and Dirac model as a function of the separation distance.Second, the coefficients of the generalized Lennard-Jones potential are determined.Finally, the pairwise summation is performed.
Although the Lennard-Jones potential is widely used to simulate graphene-based NEMS, the marginal interaction distance must be determined.It is well known that the Lennard-Jones (6,12) potential is inadequate for calculating the interaction energy at large separations.The pairwise summation method of the Lennard-Jones (6,12) potential between a graphene sheet and silicon substrate unexpected gave a good approximate interaction energy near a ≈ 100 nm.However, it should be noted that the Lennard-Jones (6,12) potential is not the correct expression for the interaction of a carbon atom in graphene and a silicon atom in the substrate for large distances between them.

FIG. 2 .
FIG. 2. Configuration of carbon and silicon atoms ((a) side view and (b) top view).The large solid circles denote the positions of the silicon atoms, and the small solid circles denote the positions of the carbon atoms.The region inside square dashed line indicates the domain of integration in (5).

FIG. 3 .
FIG. 3. Potential energy between a silicon plate and a single graphene sheet obtained by summation of the Lennard-Jones potential between a silicon and carbon atom for lmax/dSi = 100, 200, and 300.The solid line shows |Ua≫σ(a)|, whose slope is −3.

F a N m 2 FIG. 4 .
FIG. 4. (a) Potential energy and (b) force between silicon plates calculated by the Lifshitz formula.The dashed lines are the straight lines with a slope of −2 and −3 for the potential energy and the force, respectively.

FIG. 5 .
FIG.5.Dependence of the local power-exponent of the potential between silicon plates on the separation distance.

FIG. 6 .
FIG.6.Plot of the relative error of the interaction energy between silicon plates obtained by summation of the Lennard-Jones potential to that obtained by the Lifshitz theory.

FIG. 7 .
FIG. 7. Dependence of the coefficients C6 (solid line) and C7 (dashed line) normalized by ϵSi-Si on the separation distance for Lennard-Jones potential between silicon atoms.

FIG. 8 .FIG. 9 .
FIG.8.Dependence of the local power-exponent of the potential between a silicon plate and graphene sheet on the separation distance.

Figure 12
shows C5,Si-C (a) and C6,Si-C (a) normalized by ϵ Si-C .By comparing with Fig. 10, similar dependences of C5,Si-C (a) and C6,Si-C (a) on the separation distance are observed VIII.CONCLUSION

TABLE I .
Lennard-Jones parameters.