Bending modulus of the rippled graphene: the role of thickness

Bending modulus is a key parameter to characterize the stiffness of materials. Commonly, it is believed that the bending modulus is closely related to the thickness as described by the thin plate theory. However, the thin plate theory fails in multilayer van der Waals materials like multilayer graphene, suggesting a more complex relationship between the bending modulus and thickness. Here, rippled graphene structures containing non-hexagonal carbon rings with different thicknesses are constructed to study the thickness-dependent bending modulus by the first-principles calculations. It is found that the bending modulus of rippled graphene depends on several factors, such as geometry, bending curvature, and thickness. Particularly, for the egg-tray graphene structures with similar structural pattern and bending curvature, i.e., eliminating the effects of structural pattern and bending curvature, the bending modulus could show a linear relationship to the thickness. Moreover, this linear relationship is very robust even in the case of changing the thickness through heteroatom doping.


Introduction
Flexible electronics technology, as an emerging technology with great potential, is promising in the fields of display [1], sensing [2][3][4], and energy storage [5,6]. It requires materials to be light, thin, and flexible. Graphene, as an attractive twodimensional (2D) material, is a competitive candidate for flexible materials due to its excellent mechanical flexibility.
Commonly, the flexibility of a material is characterized by bending modulus. However, it is difficult to accurately measure the bending modulus of graphene. The reported bending modulus of graphene varies in a wide range from 0.8 to 10,000 eV [7][8][9][10]. The large difference in bending modulus of graphene can be ascribed to both the sample quality and measurements.
So far, it is common sense that the mechanical strength of graphene is closely related to its structural details. For example, graphene grain boundaries (GBs) can affect the mechanical strength of graphene [11][12][13][14]. Depending on the density and detailed arrangement of defects, the intrinsic strengths of graphene GBs vary from 46 to 95 GPa [12], lower than those (102 GPa in zigzag direction and 113 GPa in armchair direction) of pristine graphene. Meanwhile, the residual functional groups in graphene fabricated from reduced graphene oxide (GO) can also decrease the mechanical strength of graphene [15][16][17]. Generally, the larger the functional group coverage is, the lower the mechanical strength will be. Compared with the ~ 1 TPa of graphene, Young's modulus of GO can be degraded to 290 GPa with increasing the functional group coverage to 70% [15].
In addition to the mechanical strength, flexibility of graphene also depends closely on geometry. When incorporating topological defects such as pentagons and heptagons, various graphene allotropes can be formed [18]. Particularly, planar graphene allotropes with pentagons and heptagons show lower bending moduli than graphene, exhibiting superior flexibility [19]. Besides, for the graphene functionalized with groups like hydroxyls and epoxides, the bending modulus is quite different from that of pristine graphene. The measured bending modulus of GO can be as low as ~ 0.026 eV [20], far below that of graphene. Further theoretical calculations show that the bending modulus of GO can be either lower or higher than that of graphene, depending on the type and coverage of the functional groups, as well as the bending direction [21,22]. For example, hydroxyls can obviously enhance the bending modulus of GO up to 17.5 eV, but epoxides can reduce the bending modulus of GO to 1.2 eV, lower than that of graphene (1.5 eV) [21]. In addition, external strain can effectively tailor the bending modulus of graphene. When graphene is stretched, the bending modulus can be significantly degraded due to the weakening of carbon bonds [10].
Traditionally, the bending modulus is described by the classical thin plate theory. According to the thin plate theory, the bending modulus is calculated by D = Eh 3 ∕12(1 − 2 ) , where D is bending modulus, E is Young's modulus, h is thickness, and υ is Poisson's ratio. From the thin plate theory, thickness is a key parameter that determines the bending modulus of a material. However, this formula applies to a perfectly glued multilayer with N identical layers [23]. When the interlayer interaction changes, this formula can be unsuccessful, such as failing in multilayer graphene [23,24]. Furthermore, lots of monolayer 2D materials are not flat and show surface corrugations, which also have a certain thickness (amplitude of corrugation). So how about the relationship between the bending modulus and thickness in a monolayer 2D material with surface corrugation since the thin plate theory could not be suitable in this case?
Previously, it was reported that for the wrinkled graphene formed of hexagonal carbon rings, the bending modulus increases by the square power of wrinkling amplitude when using cosine functions to model the wrinkling of the surface [25]. This square power of rippling amplitude-related bending modulus is obtained based on the von Karman theory considering the nonlocal stress. In addition to modelling wrinkling of graphene by using cosine functions, incorporating topological defects such as pentagons and heptagons can also induce wrinkling as reported in rippled graphene allotropes. Through intensive literature research, plenty of rippled graphene allotropes with varied thickness have been reported through introducing topological defects in graphene, such as graphene with Stone-Wales defects [26], CG568 [27], egg-tray graphene [28], fused C26 polyhedra [29], Hexa-C20 [30], and penta-graphene [31]. Commonly, these rippled graphene allotropes have different carbon rings and structural patterns. The thickness of rippled graphene allotrope comes from the positive and negative curvature in geometry. For example, pentagon and heptagon will introduce positive and negative curvature separately in graphene [32]. In addition, the thicknesses and bending moduli of these rippled graphene allotropes with different topological defects and structural patterns will be significantly different.
Meanwhile, due to the high symmetry and stability of egg-tray graphene (formation energy is only ~ 0.2 eV/atom), it is easy to dope and construct derivative configurations with similar topological defect structures, which is a typical and ideal theoretical model of rippled graphene allotropes. Therefore, in this work, mainly taking the egg-tray graphene as an example, we studied their relationship between the bending modulus and thickness using the first-principles calculations. The calculation results demonstrate that the bending modulus of egg-tray graphene depends on several factors, such as geometry, bending curvature, and thickness. Particularly, when eliminating the effects of structural pattern and bending curvature, the bending modulus could be linearly related to thickness. Moreover, this linear relationship is still kept under doping modification.

Computational methods
The calculations are performed employing the plane-wave pseudopotential technique as implemented in the Vienna ab initio simulation package (VASP) [33]. The exchange and correlation interactions are described by the Perdew-Burke-Ernzerhof (PBE) parameterization of the generalized gradient approximation (GGA) [34,35]. To ensure convergence of the total energy and stress, values for kinetic energy cutoff are 400 eV for undoped structures and 600 eV for doped structures. A uniform spacing of 0.03 Å −1 of the Monkhorst-Pack K-point grid is adopted to sample the first Brillouin zone for each structure shown in Table 1. To avoid the interaction between the structure and its periodic images, a vacuum thickness of 15 Å is adopted. All the structures are optimized using the conjugate gradient algorithm until the force on each atom is less than 0.05 eV/Å and the energy is converged to 10 −5 eV. Here, we use a tube model to simulate the bent egg-tray graphene. Then, the bending modulus D is calculated by the equation W = 1 2 D 2 [10], where W is strain energy density calculated by the energy difference between the tube model and the unbent model divided by the surface area of the structure, and κ is bending curvature calculated by the reciprocal of tube radius. To ensure that the strain energy is just the bending energy, the non-self-consistent calculations are performed for the nanotube models. All calculations are carried out under the standard conditions of zero temperature and zero pressure.

Role of structural pattern
As mentioned above, many rippled graphene allotropes have been reported. However, the relationship between the bending modulus and thickness is very complex for the rippled graphene allotropes since they show different structural patterns, such as carbon hybridization, topological defect types, and positions. The bending modulus of rippled graphene allotrope is related to the factors including structural pattern, thickness, and bending curvature. For example, for the planar graphene allotropes with different structural patterns, even they have the same thickness (0 Å), their bending moduli are obviously different [19]. To further demonstrate the effect of structural pattern, we directly calculate the bending moduli of some reported graphene allotropes with different thicknesses, as listed in Table 1.
The bending modulus as a function of thickness for the rippled graphene allotropes listed in Table 1 is plotted in Fig. 1. It can be noticed that the rippled graphene allotropes with different thicknesses and structural patterns have obviously different bending moduli. Moreover, the relationship between the bending modulus and thickness is very complex, which cannot be fitted with a simple function such as linear or quadratic function. Thus, how about the case if the rippled graphene changes the thickness without altering the structural patterns? One possible rippled graphene to study this question is the egg-tray graphene, which was predicted and named by Liu et al. in 2019 and might be fabricated via bottom-up chemical synthesis methods with precursors containing pentagonal and heptagonal rings [28].

Egg-tray graphene models with different thicknesses
Here, two kinds of rippled graphene allotropes, egg-tray graphene t 1 and t 2 , are chosen to study the dependence of bending modulus on thickness. Both of the egg-tray graphene t 1 and t 2 belong to tetragonal lattice. Moreover, they have the similar pattern of non-hexagonal rings with one pair of nonadjacent pentagons (highlighted in gold) and the other pair of nonadjacent heptagons (highlighted in red) in a unit cell, as shown in Fig. 2a and f. The pentagons and heptagons induce curvature and contribute to the thickness of egg-tray graphene. The egg-tray graphene t 1 and t 2 only show slight differences in the arrangement of hexagons.
The egg-tray graphene is a periodically rippled structure with a certain pitch length. Then, the thickness can be changed by varying the pitch length through adding graphene nanoribbons. Figure 2b-e show the structures with different thicknesses by adding one to four graphene nanoribbons along x direction in t 1 configuration, as labeled t 1 -1, t 1 -2, t 1 -3, and t 1 -4. The same operation is applied to t 2 configuration, as shown in Fig. 2g-j. Table 2 shows the change of thickness during adding graphene nanoribbons. The thickness of egg-tray graphene increases with the number of graphene nanoribbons, but the increment is gradually decreased. Based on this limitation Fig. 1 The relationship between the bending modulus and thickness for the graphene allotropes listed in Table 1. The bending modulus of graphene (blue star) is also presented for comparison and computational cost, here the maximum number of added graphene nanoribbons is four. But even so, adding four graphene nanoribbons can lead to more than 150% increase in thickness, such as from 1.638 Å in t 1 to 2.506 Å in t 1 -4, and from 1.998 Å in t 2 to 3.348 Å in t 2 -4. Besides, adding graphene nanoribbons will increase the stability of the egg-tray graphene, as indicated by the decreased formation energy in Table 2, which is calculated by the energy difference between the egg-tray graphene and flat graphene.

Bending moduli of the egg-tray graphene models
To focus on the role of thickness, in addition to eliminating the influence of structural differences, it is necessary to avoid other factors that can affect the bending modulus, such as bending curvature [24]. Fortunately, for the eggtray graphene supercells constructed in Fig. 2, the lattice parameter b that determines the bending curvature is almost the same since the graphene nanoribbons are added in x direction. From the optimized lattice parameters in Table 2, it can be noticed that for both the t 1 and t 2 cases, the difference in lattice parameter b is less than 5% despite different thicknesses. Moreover, all the ten structures shown in Fig. 2 have a close lattice parameter b due to their highly similar structural patterns. Therefore, these egg-tray graphene models simultaneously have the highly similar structure pattern and bending curvature, which are good candidates to explore the relationship between bending modulus and thickness. The calculated bending moduli of the egg-tray graphene supercells are plotted in Fig. 3. As shown in Fig. 3a and b, for both the t 1 and t 2 cases, the bending modulus increases linearly with the thickness. Particularly, when combining the t 1 and t 2 cases, i.e., considering the bending modulus of all the ten structures presented in Fig. 1, the bending modulus is still linearly related to the thickness (Fig. 3c) since they have similar structural pattern and bending curvature.
In the classical thin plate theory, there is a limiting case. When the interface is frictionless, each layer would bend independently so that the bending modulus of the multilayer would scale linearly with the layer number. Inspired by this limiting case, we suppose a possible explanation for the linear relationship between the bending modulus and thickness. For the egg-tray graphene structures studied here, they almost have the same pattern of pentagons and heptagons. The change of thickness is due to adding different numbers of identical graphene nanoribbons. Under bending, the added graphene nanoribbons contribute equally to deformation in regard to the bending strain at the same bending curvature. Table 2 Optimized structural parameters and formation energies of the egg-tray graphene models shown in Fig. 2. The N C is the total number of carbon atoms in the unit cell, a and b are the lattice parameters along x and y directions separately, h is the thickness, and E f is the formation energy per atom  . 3 Linear relationship between the bending modulus and thickness of the egg-tray graphene shown in Fig. 2: (a) t 1 cases, (b) t 2 cases, and (c) t 1 and t 2 cases. Here, r 2 is the square of the correlation coefficient Therefore, the strain energy shows a linear relationship to the volume of the egg-tray graphene, leading to a linearly increased bending modulus to the thickness. However, it is easy to judge that this linear relationship will be greatly limited, and it cannot be directly applied to other rippled graphene allotropes with different structural patterns. Even for egg-tray graphene, if the structural pattern is changed by randomly adding graphene nanoribbons ( Figure S1), this linear relationship will also break, as shown in Figure S2.

Doping effect
To further confirm and explore this linear relationship, egg-tray graphene structures are doped because doping can also change the thickness but maintains the geometrical structure. Here, doping of boron (B) and nitrogen (N) and their co-doping are considered since B prefers to be doped at pentagonal sites and N prefers to be doped at heptagonal sites [44]. We first investigate the optimal doping sites for B and N doping and their co-doping by calculating the total energy, which are presented in Figure S3 and Tables S1 and S2. It was found that both the B and N doping will decrease the thickness of egg-tray graphene, as shown in Fig. 4a and b. However, co-doping of B and N atoms almost does not change the thickness of egg-tray graphene due to the bonding of B and N atoms. Then, the bending modulus of doped egg-tray graphene is calculated. Particularly, for both the single-atom doping and co-doping, although their bending moduli and thicknesses changed almost simultaneously, the bending moduli of doped egg-tray graphene still show a linear increase to the thicknesses, as shown in Fig. 4c. Therefore, the linear relationship between the bending modulus and thickness for egg-tray graphene is well maintained under doping due to their almost unchanged lattice parameters and structural patterns.

Validity of tube models
In our calculation, the bending moduli of all rippled graphene allotropes are calculated by the common tube model [19,45]. The advantage of the tube model is to ensure uniform bending because a circular tube has a uniform bending curvature. However, since the rippled graphene allotrope has a rippled structure, the tube model will show a polygonal cross section, as shown in Fig. 5. Considering that our tube model has a polygonal cross section that is different from the commonly circular shape, does this tube model affect the relationship between the bending modulus and thickness? To further verify this relationship, we consider more cases by changing the supercell size and calculation program. The 1 × 5 × 1 and 1 × 6 × 1 supercells are built and curled into tube models to calculate the bending modulus.
The results are plotted in Fig. 6 and compared with those of 1 × 4 × 1 supercells. It can be noticed that for both the t 1 (Fig. 6a) and t 2 (Fig. 6b) cases, the linear relationship between the bending modulus and thickness is well kept for different supercell sizes. Moreover, the different bending moduli indicate that bending curvature can also play a role in determining the bending modulus of a material. Another problem is that small supercells such as 1 × 4 × 1, 1 × 5 × 1, and 1 × 6 × 1 supercells could not ensure a uniform bending curvature and may undergo a large partial deformation, making the thickness of a bent tube quite different from its unbent 2D counterpart. To make a reasonable tube model, a large supercell should be adopted, such as 1 × 18 × 1 Fig. 4 Thicknesses of the doped egg-tray graphene: (a) t 1 cases, and (b) t 2 cases; and (c) the relationship between the bending modulus and thickness. Here, r 2 is the square of the correlation coefficient and 1 × 20 × 1 supercells, as shown in Fig. 5. However, large supercells have large carbon atoms (even up to thousands), which are beyond our DFT computational capabilities.
Here, the bending moduli of large egg-tray graphene supercells are calculated using the GULP program [46,47] with Brenner potentials [48]. Taking the egg-tray graphene t 1 cases (t 1 -1, t 1 -2, t 1 -3, and t 1 -4) as examples, the calculated bending modulus as a function of thickness is plotted in Fig. 7. It can be noticed that when using the large supercell to curl into a nearly circular tube, the bending modulus also shows a linear relationship to the thickness.
Based on the data from DFT and classical force field method, the tube model can well describe the relationship between the bending modulus and thickness for the eggtray graphene with similar structural pattern and curvature. Due to the role of bending curvature, the value of bending modulus at a certain thickness can be different when using different dimensions of supercells. However, the linear relationship is well kept in different dimensions of supercells.
Generally, although the bending modulus will increase with the thickness of a material, the relationship between the bending modulus and thickness can be very complex. The classical   6 Linear relationship between the bending modulus and thickness calculated by DFT for the egg-tray graphene models shown in Fig. 2 with different supercell sizes: (a) t 1 cases and (b) t 2 cases. Here, r 2 is the square of the correlation coefficient thin plate theory gives that bending modulus is proportional to the three power of thickness for perfectly glued multilayer with N identical layers, but it fails in multilayer van der Waals materials like multilayer graphene [23,24]. For wrinkled graphene formed of hexagonal carbon rings, the bending modulus increases with the square power of wrinkling amplitude when using cosine functions to model the wrinkling of the surface [25]. Here, we demonstrate that for the rippled graphene allotropes with pentagons and heptagons, the bending modulus depends on the structural patterns and bending curvature. When eliminating the effects of structural pattern and bending curvature, the bending modulus of egg-tray graphene could show a linear relationship to the thickness.

Conclusions
We study the relationship between the bending modulus and thickness of the rippled graphene. Generally, the bending modulus of rippled graphene not only depends on the thickness, but also is related to geometric details and bending curvature. Therefore, the relationship between the bending modulus and thickness can be very complex. However, taking the egg-tray graphene as an example, the bending modulus of rippled graphene could show a linear relationship to the thickness when the effects of structural pattern and bending curvature are eliminated. Moreover, this linear relationship is very robust even in the case of changing thicknesses through slightly heteroatom doping and different bending curvatures.