Universal description of potential energy surface of interlayer interaction in two-dimensional materials by first spatial Fourier harmonics

We propose a hypothesis that the potential energy surface (PES) of interlayer interaction in diverse 2D materials can be universally described by the first spatial Fourier harmonics. This statement (checked previously for the interactions between graphene and hexagonal boron nitride layers in different combinations) is verified in the present paper for the case of hydrofluorinated graphene (HFG) bilayer with hydrogen bonding between fluorine and hydrogen at the interlayer interface. The PES for HFG bilayer is obtained through density functional theory calculations with van der Waals corrections. An analytical expression based on the first Fourier harmonics describing the PES which corresponds to the symmetry of HFG layers is derived. It is found that the calculated PES can be described by the first Fourier harmonics with the accuracy of 3\% relative to the PES corrugation. The shear mode frequency, shear modulus and barrier for relative rotation of the layers to incommensurate states of HFG bilayer are estimated. Additionally it is shown that HFG bilayer is stable relative to the formation of HF molecules as a result of chemical reactions between the layers.


I. INTRODUCTION
Since the discovery of graphene [1] this material has attracted considerable attention due to its unique physical properties. The interaction between graphene layers is responsible for the tunable band gap [2] and such phenomena as superconductivity in twisted graphene bilayers [3], commensurateincommensurate phase transition [4] manifested through formation of a network of domain walls [5] with topologically protected helical states [6,7] in bilayer graphene, selfretraction of graphene layers [8,9] and so on. Graphene applications based on interlayer interaction such as nanoelectromechanical systems (NEMS) composed of graphene layers which slide with respect to each other have been proposed [10][11][12]. In addition to graphene, a wide family of other 2D materials has been recently synthesized including hexagonal boron nitride (h-BN, see Ref. [13] for a review), graphane [14], various transition metal dichalcogenides (see Ref. [15] for a review), phosphorene [16], borophene [17], germanene [18], etc. Heterostructures consisting of layers of different 2D materials should be also mentioned (see, e.g., Ref. [19] on graphene/h-BN nanoscrolls and Ref. [20] for a review). An important characteristic of interlayer interaction is the potential energy surface (PES) that is the interlayer interaction energy as a function of the coordinates describing the relative displacement of the layers. Particularly the PES determines the commensurate-incommensurate phase transition in 2D bilayers, self-retraction of 2D layers and operation of NEMS based on relative motion of such layers.
The first-principles calculations for graphene [21][22][23][24][25] and * allexandrleb@gmail.com † liv ira@hotmail.com ‡ popov-isan@mail.ru § andrey.knizhnik@gmail.com ¶ poklonski@bsu.by h-BN [26] bilayers and graphene/h-BN heterostructure [27][28][29] have shown that the PESs of interlayer interaction in such systems can be described using the first spatial Fourier harmonics. It is interesting to note that while the amplitude of corrugations of the PES for graphene bilayer computed using the simple Lennard-Jones potential is an order of magnitude less than for the surface that follows from the density functional theory (DFT) calculations, it is described excellently by the same expression [22]. The approximation by the first Fourier harmonics also works well for PESs of interwall interaction of infinite and commensurate carbon nanotube walls [30][31][32][33][34][35] and nanotube walls with edges [36] and defects [31], both obtained from first principles [32][33][34][35][36] and using classical potentials [30,31]. Based on the results listed above for PESs of interlayer interaction between 2D layers and nanotube walls, we propose here the hypothesis that the possibility of approximation of the PES by the first Fourier harmonics is a universal property for diverse 2D materials.
This hypothesis leads to the important conclusion that physical properties determined by different regions of the PES are interrelated, that is the measurement of any physical property determined by the PES gives the information about the whole PES and, therefore, can be used to estimate other properties determined by this PES [23]. In the case of interaction between graphene layers, the proposed hypothesis is not only based on the DFT calculations but also confirmed by the following experimental data. The barrier for relative motion of graphene layers derived through measurements of the shear mode frequency (related with the PES near its minimum) and the stacking dislocation width (related with the PES along the path between neighbouring minima through the saddle point) equals 1.5-1.8 meV [23] and 2.4 meV [37] per carbon atom of one layer (per carbon atom of the upper/adsorbed layer), respectively. The difference between these two estimates of the barrier is less than the scatter of the values that follow from the DFT calculations ranging from 0.5 to 2.4 meV per carbon atom of one layer, see Ref. [38] and references therein. Approximations by Fourier harmonics beyond the first ones have arXiv:2007.11289v1 [cond-mat.mes-hall] 22 Jul 2020 been also considered for graphene [39], h-BN [39] and MoS 2 [40] bilayers and graphene/h-BN heterostructure [39].
The examples where the PES is described using the first Fourier harmonics listed above correspond to the simplest 2D materials, graphene and h-BN. In the present paper, we consider this approximation for the PES of chemically modified graphene layers by the example of interaction between hydrofluorinated graphene (HFG) layers with hydrogen bonding between fluorine and hydrogen atoms at their interface. Bilayer and multilayer systems consisting of chemically modified graphene layers such as graphane, fluorographene and HFG in different combinations have been studied recently because of their interesting electronic properties and possible applications in nanoelectronics [41][42][43][44][45]. HFG is a so-called janus nanostructure with piezoelectricity within one layer [46][47][48]. This 2D material has been synthesized [49] and the structure and electronic properties of the monolayer have been investigated theoretically in Refs. [45][46][47][48][50][51][52]. The piezoelectric enhancement has been predicted for HFG bilayer [45]. While the energies of interlayer interaction have been found for some symmetric stackings of hydrofluorinated graphene bilayer [45,53] and graphane/fluorographene heterostructure [43], the whole PES has not been yet considered for the systems of chemically modified graphene layers. Based on the PES approximation, we also estimate properties of HFG bilayer associated with relative sliding of the layers: shear mode frequency, shear modulus and barrier for relative rotation of the layers to incommensurate states.
The DFT study of the interaction between HFG layers allows us to consider an additional problem. The chemical modification of 2D layers makes possible chemical reactions at the interface between the layers. However, such a possibility has not yet been addressed. Here we demonstrate that HFG bilayer is stable relative to formation of hydrogen fluoride (HF) molecules as a result of chemical reactions between the layers.
In the following, we first give the details of our DFT calculations. In Section IIIA, the results on structure of functionalized graphene monolayers and HFG bilayer in different symmetric stackings are presented. Then we consider the PES of HFG bilayer and its approximation by the first Fourier harmonics. The characteristics of HFG bilayer associated with relative in-plane motion of the layers are estimated in Section IIIC. In Section IIID, the stability of HFG bilayer relative to formation of HF molecules is addressed. Finally, we discuss PESs for different 2D bilayers.

II. COMPUTATIONAL DETAILS
The DFT calculations are performed using the VASP code [54]. The projector augmented-wave method (PAW) [55] is applied to describe the interactions of valence electrons with atomic cores. The trigonal unit cell including 4 atoms of each layer (2 C atoms, 1 H atom and 1 F atom in the case of HFG) and having the height of 25 Å is considered under periodic boundary conditions. A dipole correction [56] is used in the direction perpendicular to the layers to cancel out interactions between periodic images. Integration over the Brillouin zone is performed using the Monkhorst-Pack method [57] with the 36 × 36 × 1 k-point grid. The maximum kinetic energy of plane waves is 600 eV. The convergence threshold of the selfconsistent field is 10 −8 eV. The second version of the van der Waals density functional (vdW-DF2) [58] is used. As follows from the previous comparison with the experimental data [38], this functional in general provides better results for the properties of bilayer graphene, graphite and h-BN related to interlayer interaction, such as shear and bulk moduli, shear mode frequencies, etc. than other functionals corrected for van der Waals interactions (PBE-D2, PBE-D3, PBED3(BJ), PBE-TS and optPBE-vdW). The structure of the monolayers has been also optimized with the exchange-correlation functional of Perdew, Burke and Ernzerhof (PBE) [59] for comparison.
To get the structure of the monolayers, both the positions of the atoms and unit cell are changed. To compute the equilibrium interlayer distances for HFG bilayer in different stackings, the positions of the atoms in the xy-plane parallel to the layers are fixed, while the size of the unit cell and positions of the atoms in the perpendicular z-direction are optimized. The geometry optimization is performed till the maximum residual force reaches 0.0003 eV/Å. The binding energy per carbon atom of one layer is computed as E b = (E bi − 2E mono )/n C , where E bi is the bilayer energy per unit cell, E mono is the monolayer energy per unit cell and n C = 2 is the number of carbon atoms in each layer per unit cell. The 2D polarization P per HFG layer is found as P = µ z /σN, where µ z is the electric dipole moment in the z-direction perpendicular to the layers per unit cell, σ is the area of the unit cell in the xy-plane parallel to the layers and N is the number of the HFG layers considered.
To obtain the PES for the coaligned HFG layers, they are placed at the interlayer distance corresponding to the groundstate AA stacking (in which the hydrogen and fluorine atoms of one layer are on top of the similar atoms of the second layer) and then rigidly shifted with respect to each other. The PES for the counteraligned layers is computed at the interlayer distance optimal for the AB1 stacking (in which the hydrogen atoms of one layer are on top of the hydrogen atoms of the second layer and which is the most energetically favourable for the counteraligned layers). The calculations are performed on the grid of 24 × 12 points, with the step of 0.187 Å and 0.216 Å in the armchair and zigzag directions, respectively. Using the PES symmetry in the zigzag direction, the grid of 24 × 24 points is finally reconstructed.
The structure and energy of the HF molecule are computed using one Γ point in the simulation box with the side of 12 Å = 1.2 nm and dipole correction in the direction of the bond.

A. Structure
First we have calculated the structure of graphane, fluorographene and HFG monolayers and compared the geometrical parameters obtained with literature data. The structure of graphene monolayers hydrogenated or fluorinated from only one side (required to study the stability of HFG bilayer) has been also computed. The chair conformation has been considered for all monolayers since it is the most favourable for graphane [60][61][62][63][64][65][66][67][68][69][70], fluorographene [60,64,65,69,71,72], and HFG [47]. Experimental observations [73][74][75][76][77][78] for fluorographite are consistent with this conformation (see also [79] for a review). Different from graphene, where all atoms lie in the same plane (disregarding long-range ripples), carbon atoms of hydrogenated and fluorinated graphene layers belong to two planes of upper and bottom carbon atoms with the distance δ between these planes (see Fig. 1a). The distance δ is referred to here as out-of-plane buckling of carbon atoms.
The geometrical parameters for graphane and fluorographene monolayers are summarized in Table I and for HFG monolayer in Table II. As seen from Table I, the computed structures agree well with the results of previous firstprinciples calculations using different functionals and experimental data. The geometrical parameters obtained with the PBE functional are exactly the same as in previous calculations using the similar approach. The account of van der Waals interactions through the vdW-DF2 functional leads to a small (within 2%) increase in the lattice constant, carboncarbon and carbon-fluorine bond lengths and a decrease in the carbon-hydrogen bond length. The out-of-plane buckling δ in this case increases for graphane, decreases for fluorographene and is almost unchanged for HFG.
The lattice constant of fluorographene is 2.5% greater than that of graphane (Table I) and the lattice constant of HFG monolayer lies in between (Table II). These differences are mostly related to changes in the carbon-carbon bond length. The carbon-hydrogen(fluorine) bond lengths are virtually the same in graphane(fluorographene) and HFG monolayer. Oneside functionalization leads to a decrease in the lattice constant and carbon-carbon lengths and an increase in the carbonfluorine and carbon-hydrogen bonds as compared to the monolayers functionalized from the both sides. The out-of-plane buckling δ is 0.46-0.49 Å for the monolayers functionalized from the both sides. It is reduced almost twice, to δ ≈ 0.27 Å, in the case of one-side functionalization.
The geometrical parameters of HFG bilayer in different symmetric stackings are summarized in Table III. As seen  from comparison of Tables II and III, the internal structure of the layers is almost unaffected by the interlayer interaction. The changes in the bond lengths, angles and out-of-plane buckling induced by the interlayer interaction do not exceed 1%. The differences in the structures of two interacting layers lie in the same range. Virtually no change in the internal structure of the layers is observed upon changing the bilayer stacking. For the optimal interlayer distances, the changes in the bond lengths for different stackings are within 0.001 Å (Table III). The calculation for the AB2 stacking (PES maxima) at the interlayer distance optimal for the AA stacking (PES minima) demonstrates that the internal structure of the layers is also poorly affected by the changes in the interlayer distance for the distances between equivalent planes of carbon atoms of the top and bottom layers in the range of 5.17-5.37 Å. As compared to the AB2 stacking at the optimal interlayer distance, the lattice constant and carbon-carbon bond length of the structure with the smaller interlayer distance are decreased by only 0.2%. Similar differences are observed in the geometries of the AA and AB2 stackings at the interlayer distance optimal for AA. This means that the changes in the internal structure of the layers can be neglected upon relative sliding at the constant interlayer distance. The computed 2D polarization of HFG monolayer (Table II) is in good agreement with the results of previous calculations [45,47]. The 2D polarization of HFG bilayer obtained in our paper, however, clearly exceeds the result of Ref. [45] but this discrepancy can be attributed to the use of different van der Waals-corrected functionals and other differences in the computational details. As compared to the monolayer, the 2D polarization per layer in HFG bilayer is enhanced by 1.7% in the AB1 and AA stackings and reduced by 0-0.2% in the other symmetric stackings according to our calculations.

B. Potential energy surface of hydrofluorinated graphene bilayers
Our calculations show that the most energetically favourable stacking of HFG bilayer is AA (Table III), the stacking in which the layers are coaligned and all hydro-TABLE I. Calculated properties of monolayer graphene hydrogenated (X = H) or fluorinated (X = F) from two or one sides a : lattice constant a (in Å), carbon-carbon length l CC (in Å), carbon-X bond length l CX (in Å), angles θ CCC and θ CCX (in degrees), torsional angles θ CCCC and θ CCCX (in degrees) and out-of-plane buckling of carbon atoms δ (in Å).    Stacking   [45,53], where close energies were obtained for the AB1 , AB2 and AA stackings (Table III). The PES maxima for the coaligned and counteraligned HFG layers correspond to the AB2 and AA stackings with the relative energy of about 7.3 meV per carbon atom of one layer (Table III). In the AB2 stacking the fluorine and hydrogen atoms located between the layers are in the "on-top" positions, while the fluorine and hydrogen atoms at the outer sides of the layers are in the middle of the hexagons. In the AA stacking, the fluorine (hydrogen) atoms of one layer are on top of the hydrogen (fluorine) atoms of the second layer. Large relative energies were also obtained for the AB2 and AA stackings in Refs. [45,53]. The relative energies of these stackings of 8-9 meV per carbon atom of one layer reported in Ref. [53] are fairly close to our results.
As can be expected, the spacings between the layers are very close in the AB1 , AB2 , AB1 and AA stackings. The distances between the equivalent planes of carbon atoms of the layers in these stackings are 5.17-5.18 Å and the distances between the planes of hydrogen and fluorine atoms at the interface are 2.18-2.19 Å (Table III). In the AB2 and AA stackings, these distances are increased by about 0.2 Å.
Let us now derive the expression for the PES of HFG bilayer described by the first Fourier harmonics. We use that the potential energy surface of an atom adsorbed on a 2D hexagonal lattice can be approximated by the first Fourier harmonics as [87] U at = U 1 2 cos(k x u x ) cos(k y u y ) + cos(2k x u x ) + 3 2 where x and y axes are chosen in the armchair and zigzag directions, respectively, k x = 2π/ √ 3a, k y = 2π/a (a is the lattice constant), u describes the relative position of the atom with respect to the lattice (point u = 0 corresponds to the case when the atom is located on top of one of the lattice atoms) and parameters U 1 and U 0 depend on the interlayer distance. The first five terms of the Fourier expansion for the interaction of atoms with a graphene layer and a (111) face of the fcc lattice can be found in Ref. [88]. In the case of HFG, we sum up interactions of CF (carbon-fluorine) and CH (carbon-hydrogen) groups with sublattices of CF and CH groups of the second layer, U = U FF + U HH + U HF + U FH .
Let us first consider the coaligned layers. In this case, U FF/HH =U 1FF/HH 2 cos(k x u x ) cos(k y u y ) U HF = U 1HF 2 cos k x u x − 2π 3 cos(k y u y ) U FH = U 1FH 2 cos k x u x + 2π 3 cos(k y u y ) Here subscript HF corresponds to the interactions of sublattices with the hydrogen and fluorine atoms located at the outer sides of the bilayer and FH to the case when the hydrogen and fluorine atoms are located between the layers. The zero displacement, u = 0, corresponds to the AA stacking. Taking into account Eqs. (2)-(4) for U = U FF +U HH +U HF + U FH we finally arrive at the expression U = U A 2 cos(k x u x ) cos(k y u y ) + cos(2k x u x ) + 3 2 where for the coaligned layers, Here we introduced the notation U 0 = U 0FF + U 0HH + U 0HF + U 0FH . Correspondingly, the energy of the AA stacking is expressed as E(AA) = U C + 9U A /2 and the energies of the AB1 and AB2 are given by E(AB1) = U C + 9U B /2 and E(AB2) = U C − 9U B /2, respectively. Therefore, the parameters of approximation (5) can be found as As seen from these equations, U A is responsible for the energy of the AA stacking with respect to the average energy of the AB stackings, while U B corresponds to the difference between the energies of the AB stackings. The PES corrugation, i.e. the energy difference between the global maximum and minimum, is given by Here we should use the energies of different stackings at the same interlayer distance. At the interlayer distance optimal for the ground-state AA stacking, we get E(AB1) − E(AA) = 0.363 meV/atom and E(AB2) − E(AA) = U max = 10.279 meV/atom. In this way we obtain the parameters of the approximation listed in Table IV. The PES described by Eq. (5) with these parameters is shown in Fig. 3a. Note that account of relaxation of the internal structure of the HFG layers upon relative sliding at the constant interlayer distance leads to FIG. 3. Interlayer interaction energy of hydrofluorinated graphene bilayer U (in meV per carbon atom of one layer) as a function of the relative displacements u x and u y (in Å) of the layers along the armchair and zigzag directions, respectively, approximated according to Eq. (5). Panels (a) and (b) correspond to coaligned and counteraligned layers, respectively. The interlayer distance is constant and equals the optimal one for the AA and AB1 stackings (point u = 0) in the cases of the co-and counteraligned layers, respectively. The energy is also given relative to the AA and AB1 stackings, respectively. TABLE IV. Parameters and quality of PES approximation by the first Fourier harmonics for different 2D bilayers at the interlayer distance d (in Å): parameters U A , U B and U C for Eq. (5) per atom of one layer a (in meV/atom), PES corrugation b U max (in meV/atom), barrier for relative sliding of the layers b U bar (in meV/atom), standard deviation from the PES obtained in the DFT calculations δU (in meV/atom), and relative deviation with respect to the PES corrugation δU/U max (in %).  Table III). For the counteraligned layers, Note that from Eqs. (6), (8), (13) and (15), it is seen that at the same interlayer distance, the following condition should be complied for the parameters U A ,U A , U C and U C : The relations similar to Eqs. (9)-(11) hold between U A , U B and U C and the energies of the AA , AB1 and AB2 stackings. At the interlayer distance optimal for the AB1 stacking, E(AA ) − E(AB1 ) = U max = 10.075 meV/atom and E(AB2 ) − E(AB1 ) = 0.331 meV/atom. The values of the parameters that follow from these relative energies are given in Table IV. The PES described by Eq. (5) with these parameters is shown in Fig. 3b.
The standard deviation of expression (5) with the parameters from Table IV from the PES obtained by the DFT calculations (this PES is shown in Fig. 1 of Supplementary Material [89]) is 0.28 meV/atom both for co-and counteraligned layers. This corresponds to 2.7% of the PES corrugation. The maximum deviation of the approximation from the DFT results is 0.39 meV/atom (the full map of deviations is shown in Fig. 2 of Supplementary Material [89]). Note that the condition given by Eq. (16) is complied well, though the calculations for co-and counteraligned layers are performed at a slightly different interlayer distance (by 0.004 Å). The difference between the right-hand and left-hand sides of this equa-tion is within 1.6%.
It is seen from Table IV that the values of the parameters U A and U B for the HFG layers are very close. The parameter U A is twice greater and U B is small compared to U A and U B . This means that the term U 1FH , i.e. repulsion of the hydrogen and fluorine atoms between the layers when they are close, dominates over the other pairwise terms (see Eqs. (6), (7), (13) and (14)).

C. Properties associated with relative sliding of the layers
The PES at a constant interlayer distance can be used to estimate a number of properties associated with relative inplane motion of the layers that can be measured experimentally [23,26,29]. Here we consider for HFG bilayer the shear mode frequency, shear modulus and barrier for relative rotation of the layers to incommensurate states.
The frequency f of the shear mode E 2g , in which adjacent layers slide rigidly in the opposite in-plane directions, can be determined from the PES curvature in a given metastable state [23,26,29] as where a is the lattice constant, U eff = (a/2π) 2 ∂ 2 U/∂u 2 x is the second-order derivative of the energy per carbon atom of one layer in energy units and µ is the reduced mass. The latter can be computed for the HFG bilayer as µ = (2m C + m H + m F )/4, where m C , m F and m H are masses of carbon, fluorine and hydrogen atoms, respectively. From Eq. (5), it follows that the PES curvatures for the AA, AB1 and AB2 stackings correspond to U eff (AA) = −2U A , U eff (AB1) = U A − 3U B and U eff (AB2) = U A + 3U B . Similar expressions hold for the counteraligned layers. From the values of the parameters listed in Table IV, we thus get that the shear mode frequencies for the AA, AB1, AB1 and AB2 stackings are very close and lie in the range of 17.4-18.5 cm −1 (Table V). These values are smaller than those reported for graphene bilayer based on the DFT calculations of 35 cm −1 [22,24] and 21-34 cm −1 (Ref. [38], depending on the functional used) and experimental studies of 28 ± 3 cm −1 [90] and 32 cm −1 [91]. They are also smaller than the DFT results for h-BN bilayer of 33-34 cm −1 [26] and 25-47 cm −1 [38] and graphene/h-BN heterostructure of 37 cm −1 [29]. The difference with the data for graphene [22][23][24]38] and h-BN [26,38] bilayers can be explained by the smaller PES corrugation (Table IV) and larger reduced mass for HFG bilayer. As for the graphene/h-BN heterostructure, it has a completely different PES [29].
The same PES curvature also determines the shear modulus, which can be estimated as [26] where σ = √ 3a 2 /4 is the area per carbon atom in the HFG layer and d = 5.17 Å is the interlayer distance. The estimated TABLE V. Shear mode frequencies f (in cm −1 ), shear moduli C 44 (in GPa) and barriers ∆U rot (in meV per carbon atom of one layer) for relative rotation of the layers to incommensurate states estimated for different stackings of HFG bilayer corresponding to energy minima. When the HFG layers are rotated with respect to each other by an arbitrary angle that does not correspond to a moiré pattern, the PES should become smooth, similar to graphene [23,92,93]. Even in the structures corresponding to moiré patterns, the PES corrugation is known to be very small [94]. Therefore, the interaction energy in such incommensurate states can be estimated as an average over the PES in the commensurate state given by Eq. (5): The barrier ∆U rot for relative rotation of the layers to incommensurate states can be find by substracting from U rot the energy in the minimum. The values of ∆U rot estimated for the AA, AB1, AB1 and AB2 stackings lie in the range 3.1-3.6 meV per carbon atom of one layer (Table V). For the same reasons as the shear mode frequency and shear modulus, these barriers are smaller than the previous predictions for graphene bilayer of 4 meV/atom [92,93] and 5 meV/atom [23], h-BN bilayer of 6.3 meV/atom [26] and graphene/h-BN heterostructure of 7.4 meV/atom [29] and 7.0 meV/atom [95].

D. Stability of hydrofluorinated graphene bilayer
The recent progress in synthesis of various chemically functionalized graphene layers allows us to propose the possibility of chemical reactions at the interface between the layers with different chemical functionalization. Here we consider the stability of the HFG bilayer relative to decomposition into graphene monolayers hydrogenated or fluorinated from only one side and HF molecules as a result of chemical reactions between the layers. In such a decomposition, one HF molecule per one unit cell of the HFG bilayer is formed. The computed energy of the HF molecule is −6.383 eV. The computed total energies of the HFG bilayer and graphene monolayers hydrogenated or fluorinated from only one side are −40.591, −16.680 and −14.701 eV, respectively, per one unit cell. Thus, our calculations show that the total energy of the HFG bilayer is lower by 2.827 eV per one unit cell than the total binding energy of the products formed upon bilayer decomposition.
That is the HFG bilayer is stable with respect to the considered reaction between the layers. This result, however, does not allow us to exclude the possibility of chemical reactions at the interface between some other chemically functionalized 2D layers.

IV. DISCUSSION AND CONCLUSIONS
The parameters and deviations of the PES approximation by the first Fourier harmonics for different 2D bilayers are summarized in Table IV including the present results for HFG bilayers with 2D polarization within one layer and the previous results for a set of 2D bilayers without 2D polarization such as graphene bilayer [23], h-BN bilayer [26] and graphene/h-BN heterostructure [29]. Both for graphene and coaligned h-BN layers, the minima of the PES corresponding to the AB stackings are degenerate and U B = 0. Furthermore, the PES corrugations for these materials given by U max = 9U A /2 (see Eq. (12)) are close in magnitude. The very small |U B | for HFG bilayer with the counteraligned layers corresponds to the PES with a small energy difference between the AB stackings of only 9|U B | = 0.3 meV/atom. The PES corrugation in the latter case is, however, almost twice smaller than for the graphene and h-BN bilayers (Table IV).
The type of the PES for counteraligned h-BN layers is similar to that for coaligned HFG layers (Table IV). In the both cases, the parameter U A is negative, which means that the AA stacking is the local minimum and one of the AB stackings is the global maximum (see Eqs. (9) and (10)). The close values of U A and U B for the HFG bilayer correspond to the close energies of the AA and AB1 stackings with the energy difference 9(|U A | − |U B |)/2 = 0.4 meV/atom. For the h-BN bilayer, the difference between the energy minima is much more pronounced, 3.1 meV/atom, and the AB1 minimum is very shallow. The PES corrugation is also one and a half greater for the h-BN bilayer compared to the HFG bilayer (Table IV).
The PES for the graphene/h-BN heterostructure is different from the ones discussed above (Table IV). Here U A is positive and U B is comparable in magnitude to U A . In this case the PES has two inequivalent maxima which correspond to the AA stacking and one of the AB stackings.
The average relative deviation of the approximation by the first Fourier harmonics with respect to the maximal corrugation for the PESs obtained in the DFT calculations is within 1% and 3% for the bilayers without and with 2D polarization, respectively (Table IV). Thus, the hypothesis proposed here that the PES of interlayer interaction in diverse 2D materials can be universally described by the first spatial Fourier harmonics is confirmed in all the cases considered so far. Note that for graphene and h-BN bilayers and graphene/h-BN heterostructure, Fourier expansions up to the third term have been also studied [39]. According to these calculations, in the cases of h-BN bilayers and graphene/h-BN heterostructure, the parameters corresponding to the second and third terms are more than an order of magnitude smaller than those for the first term. For bilayer graphene, the difference is by a factor of five. Therefore, the results of Ref. [39] also confirm that the first Fourier harmonics are sufficient to describe the PESs of these materials.
For further confirmation of this hypothesis, the set of considered 2D materials should be extended. It would be now interesting to test transition metal dichalcogenides. DFT calculations of energy profiles for in-plane sliding pathways between the symmetric stackings [96][97][98] and PESs [96,98] have been recently performed for MoS 2 [40,[96][97][98] MoSe 2 [98], and MoTe 2 [98] bilayers. The energy profiles and PESs obtained seem to be qualitatively of the same shape as the PES of HFG bilayer calculated here and the energy profiles and PES for h-BN bilayer [26]. Indeed as long as 2D layers consist of two types of alternating units arranged in the honeycomb lattice, it can be expected that the same Eq. (5) holds for the PES. For MoS 2 , the PES has been approximated by the first three terms of the Fourier expansion [40]. According to this approximation, the parameters corresponding to the second and third terms are 6-7 times smaller than for the first one. There is a good chance that the PESs of interlayer interaction for other transition metal dichalcogenides can be also accurately approximated by the first Fourier harmonics. However, these materials are beyond the scope of the present paper and will be considered elsewhere.
Even though in the present paper we limit ourselves to consideration of the PES of HFG bilayer at the constant interlayer distance, the approximation derived can be useful for modeling of a number of properties and phenomena associated with relative sliding of the layers and involving small changes of the interlayer distance. Let us first discuss the phenomena where the changes of the interlayer distance can be neglected and then the way how the PES approximation can be extended to take into account the dependence on the interlayer distance.
It is reasonable assume that the interlayer distance is constant if the relative displacement takes place close to the PES minima. As examples of such properties, we have estimated the shear mode frequency and shear modulus of HFG bilayer (Table V). It is also known from the previous DFT calculations [39] for graphene and h-BN bilayers as well as graphene/h-BN heterostructure that the optimal interlayer distance depends on the relative displacement of the layers in the way very similar to the potential energy. This means that the interlayer distance should not change much if the relative displacement of the layers lies far away from the PES maxima. For HFG bilayer, the barrier for relative displacement of the layers between adjacent PES minima is small compared to the PES corrugation, similar to graphene and h-BN bilayers (Table IV). Thus, the same as for graphene and h-BN bilayers [39], it can be expected that the interlayer distance does not change much along the minimum energy path between adjacent energy minima. In such a case, the PES at the constant interlayer distance can be used to model formation of domain walls between commensurate domains with AB1 and AB2 stackings in the supported bilayer when the size of commensurate domains is much larger than the domain wall width [4,5,26,29,37,[99][100][101][102][103] and phenomena related to sliding of a flake on the 2D layer of the same material such as atomic-scale slip-stick motion of the flake attached to a microscope tip [104][105][106] and diffusion of the flake in the commensurate state [92,93]. We have also roughly estimated the barrier to relative rotation of the HFG layers to incommensurate states (Table V).
Our previous DFT calculations [22] revealed that relative energies of symmetric stackings of graphene bilayer depend on the interlayer distance in the same exponential way. This means that the terms corresponding to the first Fourier harmonics can be multiplied by an exponential factor to describe the PES dependence on the interlayer distance. We, however, leave verification of this fact for HFG bilayer beyond the scope of the present paper. The PES approximation with account of the dependence on the interlayer distance can be useful, for example, for modeling of structure and energetics of moiré patterns. In the limit of large spatial periods, bilayer superstructure corresponds to domain wall networks in which the size of commensurate domains is much greater than the domain wall width. Such domain wall networks can be described analytically [4,5,26,29,[101][102][103]. Atomistic models [107][108][109][110][111] are extensively used for simulations of bilayer superstructures in the opposite limit of small spatial periods. The PES approximation by the first Fourier harmonics with account of the dependence on the interlayer distance makes possible development of continuum models [27,28] adequate for studies of intermediate cases. In particular, such models can be employed to simulate the structures formed upon relative rotation of the layers by the angles of about 1 • , at which superconductivity was discovered for twisted bilayer graphene [3].
We have also shown in the present paper that HFG bilayer is stable relative to decomposition into graphene monolayers hydrogenated or fluorinated from only one side and HF molecules as a result of chemical reactions between the layers. However, reactions between other types of functionalized 2D layers cannot be excluded and may require further investigation.
The raw data required to reproduce our findings are available to download from Ref. [112].