An Atomistic-Based Nonlinear Plate Theory for Hexagonal Boron Nitride

Through the continuity of the DREIDING force field, we propose, for the first time, the finite-deformation plate theory for the single-layer hexagonal boron nitride (h-BN) to clarify the atomic source of the structure against deformations. Divergent from the classical Föppl-von Karman plate theory, our new theory shows that h-BN’s two in-plane mechanical parameters are independent of two out-of-plane mechanical parameters. The new theory reveals the relationships between the h-BN’s elastic rigidities and the atomic force field: (1) two in-plane elastic rigidities come from the bond stretching and the bond angle bending; (2) the bending rigidity comes from the inversion angle and the dihedral angle torsion; (3) the Gaussian rigidity only comes from the dihedral angle torsion. Mechanical parameters obtained by our theory align with atomic calculations. The new theory proves that two four-body terms in the DREIDING force field are necessary to model the h-BN’s mechanical properties. Overall, our theory establishes a foundation to apply the classical plate theory on the h-BN, and the approach in this paper is heuristic in modelling the mechanical properties of the other two-dimensional nanostructures.


Introduction
Hexagonal boron nitride (h-BN) and graphene have the same lattice structure, excellent mechanical and electrical properties, and many potential applications [1][2][3][4][5]. Compared with graphene, h-BN's mechanical property research is still not intensive [6,7]. The mechanical properties of graphene and carbon nanotubes have been researched for around thirty years. However, two critical mechanical problems have not yet been thoroughly clarified [1,[8][9][10][11][12][13][14]. First, how does graphene resist out-of-plane deformations? Second, is the classical elasticity theory still valid for the one-atom-thick two-dimensional (2D) nanostructures? In the classical elasticity theory, an initial flat 2D structure can model as the film or the plate. The differences between the two structures are as follows [15,16]: the film can only resist in-plane tension and shear, while the plate can resist in-plane deformations and out-of-plane bending and twisting. Atomic calculations show that one-atom-thick graphene and h-BNs can resist in-plane and out-of-plane deformations [1,7,12,14]. Therefore, it is natural to use the classic plate theory to describe their mechanical properties. Based on the above assumption, many theories have been developed [1,7]. However, applying the plate theory to single-layer atomic structures may induce two problems. First, the macroscopic plate bending resistance comes from the opposite tension and compression condition on both sides of the midplane [16]. For one-atom-thick graphene or h-BNs, it is impossible for an atom to be in both tension and compression. Second, the graphene's bending resistance also induces another widely mentioned argument, namely, the so-called Yakobson paradox: the uncertainty of the single-layer atomic structure's thickness leads to the uncertainty of the bending mechanical properties (corresponds to the uncertainty of the bending rigidity) [1,9]. Overall, the two critical problems mentioned above cause many controversies when the classic plate theory is applied to one-atom-thick nanostructures.
In the classical elasticity theory, the deformation energy density of the Föppl-von Karman plate is [15,16] Here 2J = ε xx + ε yy , Q = ε xx ε yy − ε 2 xy are the invariants of the 2D strain ε ij in the rectangular coordinate system (O; x, y). Lets the plate midplane displacement be w(x, y) in the z direction, then H ≈ ∂ 2 w/∂x 2 + ∂ 2 w/∂y 2 /2 and K ≈ ∂ 2 w/∂x 2 ∂ 2 w/∂y 2 − ∂ 2 w/∂x∂y 2 are the average curvature and the Gaussian curvature of the midplane. In Equation (1), k p b and k p g are the in-plane rigidity's parameters, k p B and k p G are the out-ofplane bending rigidity and the Gaussian curvature rigidity (the torsional rigidity). The rigidity coefficients in Equation (1) are the functions of the elastic constants and the plate thickness [15,16]: Here, E and ν are the three-dimensional (3D) Young's modulus and the Poisson's ratio, and h is the plate thickness. According to the classical macroscopic plate theory, has Equation (2) indicates that the out-of-plane bending and torsional rigidities are related to the plate's thickness. In detail, a thicker plate has a more robust ability to resist the out-ofplane deformations, and a plate with zero thickness cannot resist out-of-plane deformations. Since one-atom-thick graphene and h-BNs cannot define an explicit thickness, the Yakobson paradox is induced. Another critical question is whether the Gaussian curvature affects the mechanical properties of h-BNs and graphene [11,[17][18][19] (this corresponds to whether the two materials have out-of-plane twisting resistance). For graphene, current researchers have presented contradictory results on this issue. Wu [17] and Atalaya [18] showed that the graphene's deformation energy does not include the Gaussian curvature term through the continuity of the first-generation Brenner atomic potential and valence bond force field. Hence the two theories mean that graphene has not the out-of-plane twisting resistance. Nonetheless, Tu [19] and Davini [20] suggested that graphene has the negative Gauss curvature rigidity by continuity of the Lenosky force field and the 2nd-generation Brenner potential. The density functional tight-binding (DFTB) calculations also showed that graphene has a negative Gaussian curvature rigidity [1,21]. To solve this inconsistency, researchers have put in a great deal of effort. Recently, Huang and co-workers obtained a modified Föppl-von Karman plate theory for graphene using the tight-binding (TB) theory of covalent bond [11]. For the first time, this theory gives the relationships between all graphene's mechanical parameters and chemical bonds under finite deformations. In detail, the rigidities against the in-plane extension and shear, which associates with 2D Young's modulus and Poisson's ratio, come from bond length changes and the misalignments of the angular of σ-orbitals after deformations. The bending resistance, which associates with the mean curvature rigidity, comes from the off-plane pyramidalization of the three adjacent π-orbitals and the twist of the adjacent π-orbitals; and the out-of-plane twisting rigidity, the Gaussian rigidity, only connects with the twist of the adjacent π-orbitals. Huang's theory provides the physical explanation of the Yakobson paradox and theoretically demonstrates that graphene can be modeled as a modified elasticity plate. This theory also provides an instructive example to model the other 2D nanomaterials' mechanical properties.
Existing experiments and atomic calculations show that the in-plane Young's modulus and bending rigidity of h-NB are slightly lower than graphene [1]. The h-BN's chemical bonds are similar to graphene: the hexagonal symmetry structure is constructed by the sp 2 -hybridized σ-bonds between nitrogen and boron, and the remaining p-orbitals form π-bonds. Therefore, a straightforward deduction is that h-BN and graphene have a similar mechanical model. Song [22] established a continuous model (CM) of h-BN nanotubes by the Brenner potential. This theory shows that h-BN has no Gaussian curvature term and cannot resist torsion. Through the UFF atom potential, Genoese obtained the h-BN's in-plane Young's modulus and Poisson's ratio that are consistent with the atomic calculations [23]. However, it is confusing that Genoese's theory leads to a negative out-ofplane Poisson's ratio if the model's bending and the Gaussian stiffnesses are compared with the classical plate theory. Current research indicates that the h-BN's mechanical parameters are closely related to the used force fields. For example, Rajan has developed a multi-body force field to consider the polar nature of the B-N bond [24]. Molecular dynamics (MD) simulations show that in-plane mechanical parameters obtained by the Rajan's force field are very close to the results of experiments and DFTB calculations. However, based on the Rajan's force field, the out-of-plane bending rigidity is 2.13 eV. That is much greater than 1.34 eV obtained by DFTB calculations [1,14]. It is necessary to use a simple atomic force field of h-BNs that can capture the primary deformation energy and the explicit physical mechanism. The TB theory has demonstrated that the Föppl-von Karman plate can describe the graphene's mechanical properties with four independent parameters [11]. However, modelling the h-BN's mechanical properties is still an open question. In this paper, we will use the DREIGING force field, which is widely used in the MD's simulations of 2D nanomaterials, to validate that h-BN can model as the modified Föppl-von Karman plate. Our research may help clarify the long-standing controversies over the mechanical properties of h-BN and graphene. The approach in this paper can be a reference for research concerning other 2D nanomaterials.

Method
If the B-N bond's heteropolarity is removed, h-BN and graphene have identical lattice structures and similar chemical bonds [25]. Hence an intuitive hypothesis is that they have similar force fields. Huang and co-workers have obtained the TB force field of graphene that includes bond stretching, the bond angle bending, the angle of off-plane pyramidalization, and the two adjacent π-orbitals twisting angle [11]. The off-plane pyramidalization angle and the inversion angle represent the four adjacent atoms' deviation from the initial plane. Under the first-order approximation, we will illustrate later that the dihedral angle change is equivalent to the π-orbital twisting angle of adjacent two atoms. Therefore, the atomic force field of h-BNs should include the bond stretching, bond bending, inversion angle, and dihedral angle because h-BN's force field is similar to graphene. The above four items can be captured by the DREIDING force field without nonbonding effect [26]: in which the four components represent respectively the potential energy induced by the bond stretching, bond angle bend, inversion angle and dihedral angle, as shown in Figure 1. By means of MD calculations with the modified parameters of the DREIDING force field, the in-plane elastic parameters and the bending stiffness are obtained. They are in good agreement with the DFTB calculations [27]. The present paper will utilize the DREIDING force field to study the h-BN's mechanical properties. The potential energy expression of the DREIDING force field is [26]: Here l i , θ i , χ i and ω i are the deformed bond length, bond angle, inversion angle, and dihedral angle. The summation ∑ is taken over all lengths or angles. θ 0 = 2π/3 and l 0 = 0.145 nm is the initial bond angle and bond length, as shown in Figure 1. Ignoring the difference between boron and nitrogen atoms, the potential energy of one atom, for example the boron atom at the point O, rewrite as Different researchers adopted different values in Equation (5). Table 1 lists several widely used data sets of h-BN and graphene (in DREIGING, the graphene's bond order is 1.5 [28], and the h-BN's bond order is 1.09 [29]). Within the existing research, the bond stretching and bond angle bending terms were adopted by most researchers. However, how to decide the inversion angle and the dihedral angle terms is a controversial question. For example, Li developed a structural mechanics approach to study the carbon nanotube's mechanical properties by using bond stretching, bond bending, and dihedral angle terms [30]. This method has been widely applied to study the mechanical properties of the h-BN and graphene [31][32][33]. However, another opposite pathway uses only the inversion angle term to represent the change in the offplane potential energy [18,24]. We will show later that both the inversion angles and the dihedral angles are necessary, and their absence may induce incorrect results.
To obtain the continuum limit form of the h-BN's deformation energy, we make two geometric presumptions: first, all atoms in the undeformed plane distribute on the continuous deformed curved surface S, as shown in Figure 1; second, we ignore the impact of the in-plane deformations on the out-of-plane deformations. The second presumption is necessary to derive the Föppl-von Karman plate theory from the classical nonlinear elastic theory [15]. In the present research, it is not necessary to distinguish boron and nitrogen atoms. Thus, only the central boron atom, the three first-neighboring nitrogen atoms, and the six second-neighboring boron atoms are considered, as shown in Figure 1.
Here we have In Equation (6), the bold letters n pqr represents the unit vector perpendicular to the plane formed by the three adjacent atoms p, q and r; d Ok is the unit vector of r Ok . We decompose n pqr into two components at point q: n τ pqr in the tangent plane of surface S, and n n pqr in the normal direction of the surface S. Due to n τ pqr n n pqr under the small curvature, we have n iOj ≈ n O and n Ojk ≈ n j . Here n O and n j represent the unit normal vectors at O and j positions. (5) is rewritten as Considering the finite displacement theory [34], has Here d 0 Oi and d 0 Oj denote the unit vector of the undeformed r 0 Oi and r 0 Oj correspondingly. Substituting Equations (8) and (9) into the first two terms of Equation (5), and through complicated but straightforward calculations, gets (calculation details refer to Refs. [11,17,22] In order to obtain the continuous limit of the two off-plane four-body terms, we approximate the deformed h-BN's surface as the Monge patch [15,35]. We supposes that the atom O stays in the initial position during deformation, and let n Ox = r x (0, 0)/|r Ox |, and n Oy = r y (0, 0)/ r Oy , here r x (0, 0) = ∂r ∂x (0, 0), and r y (0, 0) = ∂r ∂y (0, 0). According to the classic finite deformation plate theory [15], approximately gets n Ox ⊥n Oy , n Ox and n Oy are parallel to x and y axis correspondingly. Under the coordinate system O; n Ox , n Oy , n O , as shown in Figure 1, the position vectors expand as [35] r 0j (x, y) = r x x + r y y + 1 2 According to the geometric presumption (2), we ignore the influence of the in-plane deformations on the out-of-plane deformations. Therefore, the deformed position vector of the point j(x, y) is simplified as According to Equation (12), the unit vectors in Equation (7) are Here, Substituting Equations (13) and (14) into inversion angle terms in Equation (7), and through complicated but straightforward calculations, gets Considering [11,19] and substituting Equations (10), (15) and (16) into Equation (7), we get the continuity potential energy density F as Here, A 0 = 3 √ 3l 2 0 /4 is the area occupied by an atom. The stiffness coefficients in Equation (17) are Through an analogy with the classical Föppl-von Karman plate theory, the Poisson's ratio ν and the 2D Young's modulus Y can written as Equation (17) shows that the continuous limit deformation energy of h-BNs has the same format as the classical Föppl-von Karman plate theory. Our new theory can also be used to model the mechanical properties of graphene if the graphene's force field parameters are used. Since different researchers use different force field parameters, we calculate the mechanical parameters of h-BN and graphene for the force field parameters of Table 1, with the results as shown in Table 2. In order to compare the present results and the existing ones, Table 2 also lists the typical values obtained by DFTB [1], the bond orbital theory (BOT) [11], CM [19,20,23], and Experiments [36,37]. By comparing our theoretical results and pre-exciting studies, it is found that the in-plane elastic parameters of graphene and h-BN obtained by our approach are in good agreement with the pre-exciting results. However, the out-of-plane bending stiffness and torsion stiffness are different with different force field parameters. This result reveals a key connection: the out-of-plane mechanical properties significantly depend on the atomic force field. For example, the torsional stiffness comes from the dihedral angle terms. Therefore, if the dihedral angle terms are omitted in the MD calculations, the torsional properties of graphene and h-BN will be missed. In Table 2, "Present with [*]" means the values are calculated by Equations (18) and (19) with the force field parameters in Ref.
[*]; DFTB means the density functional tight-binding calculations; BOT means the bond orbital theory; and CM means continuous model.

Discussion
Equation (17) indicates that the h-BN's deformation energy density is consistent with the classical Föppl-von Karman plate theory. However, Equation (18) displays that the h-NB's in-plane parameters, k b and k g , are independent of the out-of-plane parameters, k B and k G . This independence stems from the different atomic mechanisms between the in-plane stiffness and the out-of-plane stiffness. The changes in bond length and bond bending angle allow the ability to resist the in-plane tension and shear. The changes in the inversion angles and dihedral angles provide the ability to resist the out-of-plane bending. Furthermore, the anti-twisted ability only comes from the changes in the dihedral angles. In the macroscopic plates and shells, the in-plane stiffness comes from the potential energy increase induced by the midplane stretching and shear, while the out-of-plane stiffness comes from the opposite stretching states in the two sides of the midplane [16]. This makes the macroscopic plate's in-plane and out-of-plane stiffness meet with Equation (3). However, From Table 2 and Equation (18), gets Equation (20) shows that the h-BN does not meet Equation (3). Therefore, it is impossible to relate the in-plane and out-of-plane mechanical parameters through the equivalent thickness. This indicates that the Yakobson paradox comes from the independence between the in-plane stiffness and the out-of-plane stiffness. From the bond orbital theory (BOT) based on the TB theory, the characteristics of the h-BN's chemical bond are similar to graphene if the polarity of B and N atoms is ignored [11,25]. Therefore, it can be inferred from the graphene's BOT that the h-BN's in-plane stiffness may come from the bond stretching and bending, while the out-of-plane stiffness may come from the changes in the π-bonds.
Up to now, researchers have widely accepted that the 2D Young's modulus and the Poisson's ratio can correctly describe the h-BN's in-plane mechanical properties under finite deformations. Nevertheless, how to describe the out-of-plane mechanical properties is still a controversial open question. Here we will use the present theory to clarify some arguments. Following the classic plate theory, we introduce three auxiliary functions as κ xx = −∂ 2 w/∂x 2 , κ yy = −∂ 2 w/∂y 2 and κ xy = −∂ 2 w/∂x∂y (known as bending curvatures [15,17]), so two bending moments and one twisting moment are Equations (18) and (21) reveal: (1) the Gaussian stiffness, k G , describes the h-BN's ability to resist twisting, and this ability only comes from the dihedral angle term of the force field; (2) the bending stiffness k B comes from both the dihedral angle and inversion angle terms; (3) the Gaussian stiffness influences on two bending moments. There are two different options for the out-of-plane terms of the used force fields. The first option considers only the inversion angle terms, and this option makes k G = 0, M xx − M yy = 0. It can be found from Equation (2) that this situation is impossible in classical plate theory because k G = 0 means |ν| = ∞. Moreover, this situation also leads to two confusing results: two bending moments are always equal, and the equivalent thickness of the h-BN depends on the type of loads [17,22]. The other option considers only the dihedral angle terms and can get k B = −k G from Equation (18). This means that the bending stiffness is identical to torsional stiffness. Furthermore M xx = k B κ xx and M yy = k B κ yy , namely, a bending moment is only related to a bending curvature. Comparing Equations (21) and (2), one can find that k B = −k G corresponds to the classical plate with Poisson's ratio ν = 0. According to the differential geometry theory, the Gauss-Bonnet integral formula has [35]: the Gaussian curvature K does not appear in the balance equations and their boundary conditions, which it does not affect the deformation results of the structures. However, if the structures have free boundaries, the Gaussian curvature appears in the boundary conditions and affects the deformations [15,16]. Therefore, calculations on an h-BN shell with free boundaries may make mistakes if the force field's four-body terms only include the inversion term or dihedral angle term. In addition, from the geometric point of view, the Gaussian curvature K = 0 means that the surface is developable. If a continuous plate theory of h-BNs does not include Gauss curvature, the theory only describes special cases in which an h-BN is curled as a developable surface.
The above discussions are based on our analytical theory. Compared with the DFT calculations, the main advantages of our new theory are as follows. First, the new theory demonstrates that the modified Föppl-von Karman plate theory can describe h-BN's mechanical properties. Second, the new theory provides analytical formulas of h-BN's mechanical parameters with the DREIDING force field. DFT calculations cannot give these relationships. Third, this theory discovers that the bending and torsion stiffnesses do not depend on the h-BN's thickness. This gives a clear physical explanation of the Yakobson paradox revealed by the DFT calculations. In fact, the Yakobson paradox exists in h-BNs and graphene, and other 2D nanomaterials [1,38,39], for example the single layer MoS 2 .
Recently, heterogeneous structures composed of the h-BN and other 2D nanomaterials have been widely considered [40,41]. The nonbonding van der Waals interaction is essential in heterogeneous structures. The van der Waals interaction between the nano curved surface and the in-surface or out-surface particles can be expressed as a unified function of the mean curvature and Gaussian curvature of the curved surface [42,43]. How to consider the van der Waals interaction in the present theory needs further research.

Conclusions
In this research, through the continuity of the DREIGING force field, we prove that the h-BN's mechanical properties can model as a Föppl-von Karman plate theory with independent in-plane and out-of-plane mechanical parameters under finite deformations. The in-plane stiffness depends on the bond stretching and bond bending; the out-of-plane bending stiffness depends on the changes in the inversion angles and the dihedral angles; the torsional stiffness (the Gaussian stiffness) only depends on the changes in the dihedral angles. Since the force field's sources of the in-plane and out-of-plane stiffness are independent, it is unnecessary to relate the in-plane and the out-of-plane mechanical parameters through the equivalent thickness of h-BNs. Our present theory provides the Yakobson paradox with an atomic explanation of the one-atom-thick 2D nanomaterials. This research is heuristic in modelling other 2D nanomaterials through the classical macroscopic elastic theory.