## Abstract

We consider how the bending stiffness of a multilayer graphene sheet relies on its bending geometry, including the in-plane length L and the curvature κ. We use an interlayer shear model to characterize the periodic interlayer tractions due to the lattice structure. The bending stiffness for the sheet bent along a cylindrical surface is extracted via an energetic consideration. Our discussion mainly focuses on trilayer sheets, particularly the complex geometry-dependency of their interlayer stress transfer behavior and the overall bending stiffness. We find that L and κ dominate the bending stiffness, respectively, in different stable regions. These results show good quantitative agreement with recent experiments where the stiffness was found to be a non-monotonic function of the bending angle (i.e., Lκ). Besides, for a given in-plane length, the trilayer graphene in the flat state (κ → 0) is found to have the maximum bending stiffness. According to our analytical solution to the flat state, the bending stiffness of trilayer graphene sheet can vary by two orders of magnitude. Furthermore, once multilayer graphene sheets are bent along a cylindrical surface with small curvature, the sheets perform similar characteristics. Though the discussion mainly focuses on the trilayer graphene, the theoretical framework presented here can be readily extended for various van der Waals materials beyond graphene of arbitrary layer numbers.

## 1 Introduction

Two-dimensional (2D) materials such as graphene are a class of atomically thin sheets that have promised extraordinary electrical, thermal, and mechanical properties [1–3]. Multilayer 2D materials, also known as van der Waals (vdW) materials, have been recently developed by vertically stacking individual 2D materials, which have provided a useful routine for the unique properties of these layers being well-exploited in electronic devices [4–6]. The vdW materials are thin and very susceptible to bending deformations, just like a few pieces of paper [7,8]. As a result, they have shown various out-of-plane deformation forms in applications, such as ripples, buckles, scrolls, folds, bubbles, tents, and so on [9–14]. Though sometimes undesired for device performances, these forms have found exciting applications for mechanics metrology, novel strained electronics, and in situ transmission electron imaging [15,16] as well as enabled a number of anomalous phenomena, including flexoelectricity, pseudomagnetic fields, and photon emission [17–19]. In these scenarios, bending stiffness of the vdW materials plays a crucial role in controlling out-of-plane geometry as well as the associated functionalities.

As for multilayer sheets, on the assumption that the layers are perfectly bonded without interlayer slip (i.e., plane sections remain plane upon bending), bending stiffness would have a cubic relationship with thickness and be linearly proportional to Young’s modulus [20,21]. However, multilayer vdW materials feature strong in-plane covalent bonding and weak interlayer vdW interactions [21,22], so that layers are readily subjected to relative slips, as evidenced by recent observations on the structural super-lubrication, self-tearing, and self-rotation [23–25]. Such slip would make the bending stiffness of vdW materials significantly overestimated by traditional plate theory [21,26–29]. More intriguingly, recent reports have revealed that the bending stiffness of multilayer graphene sheets is a variable depending on their bending geometry, rather than a constant material property [27,29]. On the one hand, for multilayer graphene sheets with prescribed layer numbers, the theoretical prediction showed that the bending stiffness increases with the increasing in-plane length [29]. On the other hand, experiments and density functional theory (DFT) calculations indicated that the bending angle has a decisive effect on the bending stiffness [27]. These findings not only support the invalid cubic relationship between bending stiffness and thickness claimed in recent reports [21,27–29] but also may explain the discrepancies among different measurements on even the same type of vdW materials in these reports [8,27,30–34].

These advances have posed a question about how the bending stiffness of a vdW material depends on its two geometric descriptors—the in-plane length *L* and the bending angle *θ*. Alternatively, to characterize the bending geometry, we may use the in-plane and out-of-plane length scale: the length *L* and the radius of bending curvature 1/*κ* (since *θ* = *Lκ*). The essence of this question is how much the bending geometry of a vdW material modifies the classical “plane sections remain plane” assumption via interlayer shear. Apparently, besides the bending geometry, modifications are expected due to the intrinsic material properties, including the bending stiffness of a single layer, the lattice constant, the in-plane stiffness, the interlayer spacing, and the interlayer shear traction-separation law. This has been implied in a recent experimental work where among graphene, hexagonal boron nitride (hBN), and molybdenum disulfide (MoS_{2}) multilayers with comparable thickness, the graphene system shows the largest Young’s modulus but the smallest bending stiffness because of its weak interlayer shear resistance [21]. In this work, we adopt the material properties of graphene sheets (which have been relatively well-studied) and investigate the collective effect of *L* and *κ* on their bending stiffness that has been unclear previously. However, one may find that the proposed theoretical framework applies to various vdW materials, where the application can be achieved simply by the substitution of specific material parameters.

The rest of the paper is organized as follows. In Sec. 2, we analyze the free body diagram of each layer in a multilayer, which is bent along a cylindrical surface. We assume fixed interlayer distance and constant amplitude for the interlayer shear traction as a sinusoidal function of the interlayer slip. The equilibrium equations for these layers are derived, and the effective bending stiffness is expressed based on the stored energy. In Sec. 3, we discuss the stability of interlayer slip for trilayer graphene as an exemplary case. In Sec. 4, we predict specific values of bending stiffness for trilayer graphene sheets and compare these with experimental measurements. A more detailed bending stiffness map is also provided based on the two geometric factors. We further explore the multilayer graphene of three or more layers in Sec. 5. Finally, the concluding remarks are presented in Sec. 6.

## 2 A Bending Model for Multilayer Graphene Sheets

When a multilayer sample is bent along a cylindrical surface, and the interlayer is perfectly bonded without slip, the deflection curve will only occur in one plane, just like a pure bending of a beam with a rectangular cross section. There exist three characteristics or hypotheses [35]: (1) The neutral surface exists, in which longitudinal fibers do not change in length, and the neutral axis (longitudinal) passes through the centroid of the cross-sectional area. As for the straight sample with a rectangular cross section, the neutral surface is the middle plane. (2) All cross sections of a straight beam remain plane and perpendicular to the longitudinal axis during deformation. (3) Any deformation of the cross section within its own plane can be neglected. However, if the interlayer action is limited, the second hypothesis above will no longer hold [21]. For example, when we bend a book along a cylindrical surface, whose interlayer may be deemed as ideal lubrication with no resistance to slip, the cross section at the free end which is a plane before deformation becomes a surface and not perpendicular to the longitudinal axis (as seen in more detailed descriptions as the Appendix). As for a multilayer graphene sheet whose interlayer is neither perfectly bonded nor ideally lubricated, the interlayer shear and slip would coexist, as shown in Fig. 1(a). Based on the observation of the multilayer sample bending along a cylindrical surface, we have the following four hypotheses in order to simplify the problem:

The material is homogeneous and linear elastic.

The middle plane of the whole multilayer sample is a neutral surface.

Any deformation of the cross section within its own plane, including the change of interlayer distance, can be neglected.

The interlayer shear response is isotropic, represented by the response in the zigzag direction.

On account of our third hypothesis, we only take interlayer shearing, one of the interlayer failure modes, into consideration but ignore two other failure modes, rippling and kink buckling [29].

Besides, we only discuss the condition of odd layers, though the scenarios of even layers may be analyzed in a similar way. Suppose that the total layer number is *N* = 2*n* + 1, *n* ∈ N^{+}. We denote the middle layer by *i* = 0, and *z*_{0} denotes the corresponding thickness direction coordinate value of the middle plane. From the middle layer to the curved convex side, *i* gradually increases, and to the concave side, *i* decreases. Thus, *i* belongs to {− *n*, …, −1, 0, 1, …, *n*}, and the corresponding thickness direction coordinate value is *z*_{i} = *id*, where *d* denotes the interlayer distance and equals 0.34 nm for multilayer graphene.

In most experiments, the thickness of multilayer graphene sheets is small relative to its in-plane size. Thus, when the sample is bent along a cylindrical surface, we treat it as a multilayer beam. Besides, for the symmetry, we only take the half into consideration, and the symmetrical cross section without relative interlayer displacement can be treated as one end clamped.

We then consider the case of the perfectly bonded interlayer, where the second hypothesis in traditional theory still holds, as shown in Fig. 1(b). The length of the longitudinal axis in *i*th layer changes from *L* to *L*(1 + *z*_{i}/*R*), scilicet *L*(1 + *z*_{i}*κ*), where *R* is the curvature radius of the middle surface after bending, and *κ* = 1/*R* is the corresponding curvature.

Having considered the two limiting cases, we move on to the real bending stiffness of a multilayer sheet. To do so, we consider the energy stored in the multilayer sheet when it is bent along a cylindrical surface. There are various ways to enforce this specific displacement field for each layer in the sheet, including sandwiching the sheet between two rigid concentric circles, applying distributed transverse pressure, and so on. Here, however, we use an alternative way that the middle layer is mapped to a cylindrical surface without change in the in-plane length [29], by which the problem becomes to solve ordinary differential equations (ODEs) for merely the in-plane displacement field of each layer. The out-of-plane displacement fields would be prescribed by the curvature of the cylindrical surface and our hypothesis (3).

*n*+ 1 individual, separated layers. Such loads guarantee that the assembly of these layers into a multilayer behaves exactly the same as the perfectly bonded sheet with a curvature of

*κ*for the middle layer. We then confine the assembled multilayer between two rigid, frictionless concentric circles with the inner radius of (1/

*κ*−

*nd*), and the outer radius of (1/

*κ*+

*nd*). Finally, we convert the interlayer from the perfectly bonded case into realistic cases simply by applying a relaxation load

*P*

_{i}at the end of each layer

*E*

_{2D}is 2D in-plane stiffness of graphene, which is taken as 340 N · m

^{−1}[3]. The use of the relaxation load cancels out the mechanical response (stretching or compression) originating from the initially applied linearly distributed load. Note that when the problem is proceeded using the initial undeformed configuration, the magnitude of relaxation load just equals to that of the initially applied load. However, we will use the deformed configuration assuming no interlayer slip as the reference, so that the denominator in Eq. (1) takes (1 +

*κz*

_{i}). This would stretch/compress the sheet into a length of

*L*if the interlayer is ideally lubricated. This process ensures that the strain energy stored in the multilayer sheet purely results from the bending of the sheet onto a cylinder. Note that we do not apply relaxation moment at the free end or relaxation pressure on the surface because these terms would not contribute to the stored energy in the deformed multilayer sheet under the assumption of rigid confinement and our hypothesis (3). In other words, the out-of-plane displacements (work conjugate of pressure), including the bending angle at the edge (work conjugate of edge moment), have been prescribed.

*τ*and relative displacement [29,36]. The sinusoidal function reflects the periodicity of the lattice, which is similar to the classic Peierls–Nabarro dislocation model [37]. For A–B stacking graphene sheets, according to our fourth assumption, we consider a case of one layer sliding on the other static layer and describe the interlayer shear traction as

*u*is the displacement of the sliding layer, $l0=3lc\u2212c$ is the periodical length where

*l*

_{c−c}is the C–C bond length and taken as 0.14 nm, and the amplitude

*τ*

_{max}can be deduced from molecular dynamics (MD) simulations on interlayer shear response [29]. $\tau max=Fmaxsydc$, where $Fmaxsy$ is the maximum average shear force per atom, and $dc=4/(33lc\u2212c2)\u224539.3Atom/nm2$ [29]. According to Fig. 11(c) in the literature [29], we can approximatively obtain $Fmaxsy=0.015nN/atom$. Thus,

*τ*

_{max}≅ 590 MPa. Besides, we ignore the edge effect, which should be trivial for typical experimental samples whose bending length is much larger than

*l*

_{0}[29].

Equation (2), along with the line of our fourth hypothesis, actually dictates that samples are bent along the zigzag direction. We shall show shortly that this simple form allows for a number of explicit understandings on the important role played by the bending geometry in the bending stiffness of a multilayer sheet. The consideration of an arbitrary bending direction, however, requires a more complicated interlayer shear traction law, which should present preciser predictions but may preclude any analytical efforts. We thus leave this in a future study.

In the following modeling process, when considering the state after bending with imaginary perfectly bonded interlayer as the reference configuration, we add a bar above the symbol of all quantities, such as $u\xaf$, except for those that do not change with the coordinate system transformation. In particular, the in-plane displacement of the middle plane of the whole sample is always zero. However, for the sake of unity, it is still denoted by $u\xaf0$.

*l*

_{0}has been considered to match the geometry of deformed lattice [29].

*σ*refers to the stress resultant rather than stress. According to the geometrical relationship on the assumption of small deformation, the strain is

For the unity of symbols, we also introduce $u~i(x~)=u\xafi(x\xafi)$ and $\tau ~i=\tau \xafi$.

*B*is the bending stiffness of monolayer graphene and taken as 1.44 eV in our theoretical analyses [38]. In particular, for the middle layer

## 3 Stability Analysis of Interlayer Slip Within a Trilayer Graphene Sheet

The trilayer sheet is a good starting point because of not only its simple formulation but also that the associated discussion could provide a basis for sheets of more than three layers. In this section, we, therefore, analyze the interfacial shear stress transfer (in particular, the stability of interlayer slip) within a trilayer graphene sheet.

### 3.1 Qualitative Theory.

*l*

_{s}= (1 −

*κd*)

*l*

_{0}/(2

*π*), and dimensionless quantities, $\u03d1=u\xaf\u22121/ls$, $\eta =x\xaf\u22121/ls$, and $\omega 1=(ls\tau max)/E2D$. Then, the governing equation can be written as

*T*=

*ω*

_{1}

*η*so that the governing equation becomes

*T*

_{p}=

*T*, $\u03d1p=\u03d1+\pi $, and $\alpha p=d\u03d1p/dTp=\alpha $, Eq. (23) is able to transformed into the equation of nonlinear pendulum [39], which may help us to understand our differential system. In the classical pendulum system,

*T*

_{p}, $\u03d1p$, and

*α*

_{p}correspond to the time, the angle of pendulum, and the angular velocity, and $\omega 12$ corresponds to the ratio of the acceleration of gravity to the length of the frictionless pivoted massless rod. A particular feature of Eq. (24) is the critical points at (

*k*π, 0),

*k*∈

*Z*, which are located at the center when

*k*is odd and the saddle otherwise.

*T*

_{b}=

*ω*

_{1}

*L*(1 −

*κd*)/

*l*

_{s}, and the subscript ‘b’ refers to the free end.

*C*

_{0}is a constant for a specific trajectory in the phase diagram, which corresponds to the total energy of the pendulum, including the kinetic and potential energy.

As shown in Fig. 2, where the reference configuration is temporarily supposed as a stress-free state, the relaxation load is applied as an active external load depending on the bending curvature in the model. Thus, the sign of the stress at the clamped end should be the same as that of the relaxation load despite the presence of interlayer shear force. In addition, the overall deformation response of the studied layer should be consistent with the relaxation load. These two judgments can also be obtained from the phase diagram of the differential equation, as shown in Fig. 3(a). For example, for the case of the curved concave side we consider, where the state after bending without interlayer slip is used as the reference configuration and the relaxation load is a tensile load, the overall displacement and strain should be positive. This is equivalent to that any place on the concave side of the trilayer sample except the clamped end correspondingly locates in the first quadrant of the phase diagram, where $\u03d1$ and *α* are corresponding to the displacement and strain, respectively. Thus, for the case of the curved concave side, *C*_{0} > 0^{2}/2 + cos 0 = 1. In the pendulum system, this critical *C*_{0} = 1 corresponds to the maximum possible potential energy, which is obtained from the highest point of the system. When the energy is larger than this critical value, the pendulum can swing over; otherwise, it can only swing back and forth [39]. In our system, the clamped end of the sample exists. In the phase diagram, as illustrated in Figs. 2(b) and 2(c), the boundary condition at the clamped end corresponds to the axis, $\u03d1a=\u03d1|T=0=0$, where “a” refers to the clamped end. This corresponds to that the pendulum is able to swing over, and the total energy must over the maximum possible potential energy. In this way, we can also obtain the result *C*_{0} > 1.

*T*increases. The effect of curvature may first appear in the relaxation stress and then the boundary condition at the free end. The second formula of Eq. (9) further determines

The combination of the boundary condition at the clamped end, $\u03d1a=0$, and Eq. (27) implies that different from the pendulum system with initial boundary conditions, our system has moving boundary conditions.

In the phase diagram, two different trajectories cannot intersect. We take two arbitrary trajectories with *T*_{b1} = *T*_{b2} but $\u03d1b1>\u03d1b2$, where the subscripts “1” and “2” refer to these two trajectories. As shown in Figs. 3(b) and 3(c), *α*_{a} increases with increasing $\u03d1b$, i.e., *α*_{a1} > *α*_{a2}. Because these two trajectories cannot intersect, $\alpha 1|\u03d1=\u03d1b2>\alpha 2|\u03d1=\u03d1b2=\alpha b2$. According to Eq. (25), $d\alpha /d\u03d1=sin\u03d1/\alpha \u22650$ when $\u03d1\u2264\pi $, while $d\alpha /d\u03d1=0$ when $\u03d1=\pi $. In the process of mechanical deformation, *κ* is the direct variable. Correspondingly, we use *α*_{b} for analyses. $\u03d1b=\pi $ corresponds to *α*_{b} ≥ 2. Thus, when *α*_{b} ≤ 2 (and $\u03d1\u2264\pi $), $d\alpha /d\u03d1>0$ (i.e., $\alpha b1>\alpha 1|\u03d1=\u03d1b2>\alpha b2$) can be ensured. This implies that *α*_{b} increases monotonically with the increasing $\u03d1b$. We suppose that the curvature increases from zero. Thus, when *α*_{b} increases from zero but lower than 2 and $\u03d1b\u2208(0,\pi )$, the system is stable without the change of monotonicity. Introduce another set of dimensional quantities, $x^=x\xaf\u22121/l0$, $L^=L/l0$, $u^=u\xaf\u22121/l0$, $\kappa ^=\kappa d$, then *α*_{b} = 2 corresponds to $\kappa ^=\kappa ^c\u22450.015976$. Besides, when *T*_{b} is less than a certain critical value, *α*_{b} becomes a monotonically increasing function of $\u03d1b$ in the first quadrant of the entire phase diagram, as shown in Fig. 3(b). Finding this critical value of *T*_{b} is an inverse problem, and we will obtain this value analytically in Sec. 3.2. It is the possible transformation of monotonicity that requires us to use $\u03d1b$ rather than *α*_{b} in boundary conditions to obtain the results in Figs. 3(b) and 3(c). When the function of *α*_{b} and $\u03d1b$ is monotonic, it is reasonable to choose either *α*_{b} or $\u03d1b$ as a boundary condition, while the choice of *α*_{b} may be more straightforward.

Furthermore, we can obtain the distribution of the relative slip displacement between layers through the phase diagram, as shown by the colored solid line in Figs. 3(b) and 3(c). Although the displacement of the interlayer slip first increases monotonically from the clamped end, we may observe the phenomenon of oscillation when the length and curvature are large enough to ensure that $\u03d1$ stably reaches π, which we attribute to the sinusoidal relationship between interlayer shear traction and interlayer relative slip displacement. This oscillation phenomenon does not exist when the interlayer shear constitutive relationship is elastic or elastoplastic [40].

### 3.2 Critical Bending Length.

*i*= −1, the boundary condition, corresponding to the governing equation, Eq. (22), is

*α*> 0 (i.e., $d\u03d1/d\eta >0$). Combining Eqs. (26) and (30), we have

*η*changes from 0 to

*L*(1 −

*κd*)/

*l*

_{s}while accordingly $\u03d1$ changes from 0 to $\u03d1b$. We take

*η*= 0 as the initial state and integrate Eq. (32), obtaining

*F*is the first kind of incomplete elliptic integral, which is defined as $F(\phi |m)=\u222b0\phi (1\u2212msin2\zeta )\u221212d\zeta $ [41]. At the clamped end

*κ*

*T*

_{b}discussed at the end of Sec. 3.1 corresponds to the case that the inflection point discussed previously coincides with the stationary point, so let

*L*, from which we may further deduce the corresponding dimensionless quantities introduced before, $\u03d1b*\u22454.188$, $\alpha b*\u22452.320$, $Tb*\u22451.994$, and $L^*\u224539.14$. Besides, we may consider the physical implication of these parameters by taking $L^*$ as 39.

According to the analyses on stability in Sec. 3.1, when $Tb<Tb*$, that is, $L^\u2264L^*$, $\u03d1b$ and *α*_{b} are mapped one on one, while when a non-monotonic function relationship between $\u03d1b$ and *α*_{b} appears, a bifurcation of saddle point appears, which corresponds to a snap-through sliding in terms of mechanical response [29].

Now, it is clear that when $L^\u2264L^c=39$ or $\kappa ^\u2264\kappa ^c=0.015976$, the interlayer slip behavior is stable. For the convenience of discussion, the phase diagram is roughly divided into four regions ($[0,\kappa ^c]\xd7[0,L^c]$, $[0,\kappa ^c]\xd7(L^c,+\u221e)$, $(\kappa ^c,+\u221e)\xd7[0,L^c]$, $(\kappa ^c,+\u221e)\xd7(L^c,+\u221e)$), and the bifurcation phenomenon can only arise in the fourth region ($(\kappa ^c,+\u221e)\xd7(L^c,+\u221e)$), as shown in Fig. 4. Because the typical experimental measurement results are in stable regions, which will be shown in Sec. 4, we focus on the stable regions. It should be pointed out that there is a small stable part in the fourth region. However, in order to briefly discuss the continuous shear behavior, i.e., the stable part, we do not discuss this stable but small region.

In the stable regions, $\u03d1=\u03d1(\eta ,\kappa )$ can be obtained from the implicit function, Eq. (37), which is equivalent to that the displacement solution *u* = *u*(*x*; *κ*, *L*) is obtained. Then, following the process of Sec. 2.2, we can deduce the effective bending stiffness *D*^{eff} = *D*^{eff}(*κ*, *L*), which implies that the bending stiffness is dominated by both curvature and bending length.

## 4 Theoretical Predictions and Comparison With Experimental Results

The theoretical analysis in Sec. 3.2 has suggested that the bending stiffness is dominated by both curvature and bending length. As for the case that the sample is bent along a cylindrical surface, the product of the curvature and the bending length is the bending angle. A question arises about whether it is possible that the bending stiffness is controlled solely by the bending angle, as previously concluded from the general trend of DFT and experimental results [27].

*l*

_{p}, which includes a straight part and two curved parts with curvature radius

*R*

_{p}and bending angle

*θ*(Fig. 5(b)). Besides,

*H*

_{p}denotes the height of the step and

*L*

_{p}denotes the length of the projection of the detached part onto the hBN substrate.

*H*

_{p},

*R*

_{p}, and

*θ*have been viewed as three independent quantities. As a result,

*l*

_{p}and

*L*

_{p}can be, respectively, described as

Obviously, the curved part in the experiment does not have two free ends (or the case that one end is clamped, and the other is free), but have adjacent straight parts. We first denote the length of one curved part as *L*_{I}, which is equal to *R*_{p}*θ*, and the half-length of the straight detached part as *L*_{II}. Then, further taking effect from the straight part on hBN into consideration, we consider an additional three times *L*_{II} and average the curvature, as shown in Fig. 5(c). Due to the symmetry, the bending length in experiments is adopted as *L*_{I} + 4*L*_{II} in our theoretical analyses, and the corresponding curvature is (*L*_{I} + 4*L*_{II})/*θ*.

We compare the theoretical prediction with the experimental measurement in Fig. 6(a). Our model agrees with the experimental results well and reveals the phenomenon that it is not necessary that a smaller bending angle corresponds to a smaller stiffness, which may conflict with the whole trend implied by other experimental data [27]. For example, as shown in Fig. 6(a), the sample with the smallest bending angle $(12deg)$ does not have the largest bending stiffness. We attribute this novel phenomenon to the effect of bending length. To further elucidate the mechanism behind this phenomenon, we consider both curvature and bending length to deduce the bending stiffness based on our model, as shown in Fig. 6(b). We also mark the corresponding experimental points with white cross-shaped symbols and plot contour of bending angle with white dashed curves. For a specific bending length, bending stiffness indeed decreases initially with the increase of curvature (or bending angle), but may also experience a slight increase when the bending angle is about $30deg\u223c40deg$. Along the contour curve of $20deg$, we observe that for a specific bending angle, the bending stiffness increases solely with the increase of bending length. Moreover, as for the sample with the smallest bending angle $(12deg)$, we reveal that it is the smallest bending length that makes the corresponding bending stiffness smaller than that of the sample, which owns a larger bending angle $(16deg)$. Further theoretical analyses on the bending stiffness of trilayer graphene based on our model will be shown in Sec. 5.1.

In our analyses, an additional 3*L*_{II} is included to consider the effect from the straight part on hBN substrate. If we choose other additional length other than 3*L*_{II}, such as 2*L*_{II} or 4*L*_{II}, the value of the predicted bending stiffness shall be different because of the different bending length. However, the whole trends based on different choices are all similar in a robust way, and agree well with the whole trend of experimental data. Furthermore, to match the experimental data on values, we choose 3*L*_{II}.

We attribute the difference between the theoretical prediction and the experimental measurement, especially when the bending angle of the sample is relatively small, to the following two factors.

First of all, in the experiment, the graphene may not be bent along the zigzag direction; that is, the interlayer shear action may differ from our fourth hypothesis to some extent. This also contributes to the difference between the theoretical prediction results and the experimental measurement results.

Despite some differences, our theoretical predictions are in good agreement with the experimental measurement results. One of our key findings is presented: the bending stiffness is not a monotonic function of the bending angle due to the fact that the bending length also plays a vital role.

## 5 Further Results and Discussion

### 5.1 Trilayer Graphene

#### 5.1.1 Theoretical Results.

In stable zones shown in Fig. 4, that is, $\kappa ^<\kappa ^c$ or $L^<L^c$, we solve our model, revealing the mechanism clearly that the bending stiffness is determined by both the bending length and curvature, as shown in Fig. 7.

We find that when the bending length is relatively large ($\kappa ^<\kappa ^c$ and $L^\u226bL^c$), or the curvature is relatively small ($L^<L^c$ and $\kappa ^\u226a\kappa ^c$), the bending stiffness is dominated by the bending length, increasing monotonically as bending length increases and only decreases insensitively with the increase of curvature. This find is consistent with the previous conclusion that the bending stiffness increases with the increase of length to some extent but extends the conclusion to a wider region, where the bending length can exceed 9.7 nm ($L^\u223c40$) [29].

Meanwhile, we can observe from Figs. 7(a) and 7(b) that when the curvature gradually decreases, regardless of the bending length, the bending stiffness always converges to an initial value, and this initial value is a monotonically increasing function of the bending length, which will be further discussed in Sec. 5.1.2.

The effect from curvature on the bending stiffness gradually becomes apparent as we further consider the condition of relatively large curvature (but $L^<L^c$) in the stable region. As shown in Fig. 7(b), when the curvature is close to the critical curvature $(\kappa ^\u223c\kappa ^c)$, it is interesting that the increase of the bending length has a small effect on the increase of the bending stiffness, and even there exists an area where as the bending length increases, the bending stiffness decreases under the same bending curvature. In these regions with relatively large curvature, the curvature dominates the change of the bending stiffness.

Here, we take microscale and nanoscale bending behaviors as two typical examples to illustrate the implication of the above-mentioned results. On the one hand, when we consider the case of a large bending length, such as $L^=104$, i.e., 2.4 *μ*m, which corresponds to a microscale bending behavior, for $\kappa ^=\kappa ^c$, the bending angle is up to 6528 deg. In other words, even if this trilayer graphene sample is wound for more than 18 turns around the cylindrical surface, the sample can still amazingly maintain its bending stiffness approximately the same as its initial value. On the other hand, when the bending length is a few nanometers, even if the bending angle is relatively small, the corresponding curvature can be considerable, and there may be a great difference between the bending stiffness and the initial value. For example, for a sample with bending length of 5 nm, when $\kappa ^$ reaches $\kappa ^c$, the bending angle is only about $14deg$. It should be noted that the bending angle mentioned earlier is considered for the half sample when the sample has two free ends, or the whole sample when the sample has one clamped end and one free end. These discussions imply that the possible changes of bending stiffness ought to be considered when we design microscale or nanoscale out-of-plane deformation of vdW materials in order to precisely induce the anomalous physical phenomena or meet the functional requirements of the devices.

#### 5.1.2 A Solution to the Flat State.

*κ*→ 0, the governing equation, Eq. (22), can degenerate into

*κd*≈ 1. Use dimensionless quantities, $u^=u\xaf\u22121/l0$, $x^=x\xaf/l0$, $L^=L/l0$, and $\kappa ^=\kappa d$, and introduce $\omega 2=2\pi \tau maxl0/E2D$. Equation (43) then becomes

Here, we take $L^\u2208[1,104]$ as an example to illustrate the bending stiffness of trilayer graphene in the flat state, as shown in Fig. 8(a). On the one hand, the bending stiffness increases monotonically with the bending length in the flat state. One the other hand, the bending stiffness seems to have an upper bound and a lower bound. The upper bound can be analytically deduced from Eq. (48) as (3*B* + 2*E*_{2D}*d*^{2}), and the lower bound is easily known as 3*B*. In other words, the bending stiffness of trilayer graphene predicted by our model can range from 4.32 eV to 495 eV, which we suppose that can partially explain the wide divergence on the measurements of bending stiffness [8,27,30–34]. Revisiting the single-plate continuum model which is modified for the discrete nature [27,31], $Dupeff=NB+E2Dd2(N3\u2212N)/12$, we find the upper bound we predict for trilayer graphene is just the case of *N* = 3. Similarly, the lower bound we predict for trilayer graphene can also be seen as in the case of *N* = 3 for $Dlowereff=NB$ under the ideal lubricated interlayer condition. What is more, on account of *B* ≪ *E*_{2D}*d*^{2}/12 [38], the ratio between the upper bound and the lower bound can significantly exceed *N*^{3}.

For the sample with a specific bending length, as the curvature increases from the flat state, the bending stiffness generally decreases, and tends to the lower bound, as illustrated in Fig. 8(b). We may combine Figs. 7 and 8 to obtain a comprehensive understanding of how the bending length and curvature jointly determine the bending stiffness.

It should be pointed out that for a specific bending length, the bending stiffness is not a monotonic decreasing function of the bending curvature in the strict sense. Here, we use the sample length $L^$ equal to 20 as an example to show the transition of monotonicity. As shown in Fig. 9(a), the effective stiffness *D*^{eff} decreases with the increase of the bending curvature $\kappa ^$ in general, but when the bending curvature $\kappa ^$ is about 0.04, there exists a small area where the monotonicity changes, which is magnified in the inset figure of Fig. 9(a).

If we multiply the curvature by the bending length, we can also obtain the figure that shows the relationship between the bending stiffness and the bending angle. The format of the figure remains unchanged, but the abscissa is linearly mapped. We additionally take the sample length $L^$ equal to 13.5 as an example and consider the change of bending stiffness with respect to the bending angle, as shown in Fig. 9(b). Because of the reduction of bending length, the bending stiffness of the sample with bending length $L^=13.5$ is generally smaller than that of the sample with bending length $L^=20$. Furthermore, in DFT calculation, the relationship between the bending stiffness and the bending angle can also be obtained by compressing the sample to form a ruck [27], and the bending curve must experience the change of concavity, which implies that the corresponding bending length in our model is smaller than the real length of sample along the bending direction in DFT calculation. By comparison, it can be found that the DFT calculation result with length $L^=20$ is similar to the result in Fig. 9(b) [27].

### 5.2 Multilayer Graphene With More Than Three Layers.

According to the discussion in Sec. 2.2, we numerically solve our model for multilayer graphene with more than three layers under the condition of small curvature. Figure 10(a) shows the relationship between the bending stiffness and bending length for graphene with five layers, revealing that the bending stiffness is insensitive to curvature when the curvature is small, which is similar to the case of trilayer—considering that the larger layer number, the greater the overall slip under the same curvature, we use $\kappa ^=0.01/N$ to represent the small curvature bending state of different layers. We observe that the bending stiffness always has an upper bound as the bending length increases. Thus, for samples with more than three layers, we approximate the bending stiffness of $L^=104$ as the upper bound and normalize the bending stiffness, as shown in Fig. 10(b). It can be observed that for the cases with different layers under small curvature bending, the increasing trends of the bending stiffness with the bending length are generally similar to that of trilayer in the flat state, and merely the bending length which is corresponding to a rapid increase of the bending stiffness increases with the increase of the layer number.

## 6 Conclusion

In this work, the mechanism of how bending stiffness of multilayer graphene with a specific layer number is jointly determined by the bending length and curvature is investigated theoretically. We establish a model to describe the multilayer graphene sheet bent along the cylindrical surface and obtain the energy expression form of bending stiffness. For trilayer graphene, our model theoretically predicts the experimental measurement results well and reveals that the bending stiffness is dominated by the bending length and curvature, respectively, in different stable regions, rather than solely by the bending angle.

As for the implications, we call for consideration of the possible change of bending stiffness, when we design the various microscale or nanoscale out-of-plane deformation patterns of vdW materials. On the one hand, when the bending length is relatively large ($L^\u226bL^c$ and $\kappa ^<\kappa ^c$), or the curvature is relatively small ($\kappa ^\u226a\kappa ^c$ and $L^<L^c$), the bending stiffness is dominated by the bending length. In this region, the bending stiffness increases monotonically as bending length increases, and only decreases insensitively with the increase of curvature. On the other hand, the effect of curvature on the bending stiffness becomes apparent when the curvature is relatively large ($\kappa ^\u223c\kappa ^c$ and $L^<L^c$). These two results imply that while the microscale trilayer graphene is able to maintain its initial bending stiffness well even if it is rolled a dozen times, the bending stiffness of nanoscale trilayer graphene has to experience a significant change when the bending angle reaches only a dozen degrees.

The flat state, in which the bending stiffness reached its maximum for a specific bending length, is analytically explored for the trilayer graphene sheet, which indicates that the bending stiffness of trilayer graphene can range from 3*B* to (3*B* + 2*E*_{2D}*d*^{2}), i.e., from 4.32 eV to 495 eV. We suppose that this large range of two magnitude orders is able to partially explain the wide divergence on the measurements of bending stiffness before. Furthermore, when multilayer graphene samples are bent along a cylindrical surface with small curvature, samples with more than three layers share similar characteristics with trilayer graphene. Although we have focused on multilayer graphene, our results should be extended to other vdW materials, giving fundamental guidance on the rational design of functional vdW structures with various out-of-plane deformation patterns.

**Acknowledgment**

This work was supported by the NSF of China through Grants Nos. 11890681, 11672301, 11521202, 11890682, and 11832010, and the Strategic Priority Research Program of Chinese Academy of Sciences through Grant No. XDB36000000. X. M. thanks Fei Pan, Pinshane Y. Huang, Guorui Wang, Yuan Hou, Hao Long, and Hongtian Qiu for their helpful discussions. We thank Zhaohe Dai for their useful comments and edits on the manuscript. We also thank Pinshane Y. Huang for data from the literature [27].

### Appendix. Bending a Book Along a Cylindrical Surface

In order to explore the basic geometric relationship of the behavior that the sample is bent along a cylindrical surface, we bend a book, as shown in Fig. 11. In the experiment, we try our best to ensure that there are no obvious gaps between the sample layers; that is, the interlayer distance (the distance between the middle planes of adjacent layers) is equal to the thickness of each layer, so as to meet the third hypothesis of our model. If we choose some points in each layer and draw straight lines perpendicular to the deflection curve of the corresponding layer through these points, then these straight lines will converge to the same point, the center of curvature, which we denote as the point *O*′.

*L*and

*t*, respectively, denote the width and thickness of the book, and

*R*denotes the curvature radius of the middle plane of book, whose value can be obtained from the relation that the radius of cylindrical top surface equals

*R*−

*t*/2. The coordinates of any point in the middle plane of the

*i*th layer after bending is

*l*denotes the length of the arc along the path of the corresponding middle plane from the point to the clamped end, and

*L*

_{i}denotes the length of the corresponding middle plane after bending. For the case that the interlayer is ideally lubricated, the middle plane of each layer remains neutral during the bending deformation, that is,

*L*

_{i}=

*L*.

As for Eq. (A1), by setting *l* = *L*, we can theoretically obtain the parametric equation of the cross section at the free end. Similarly, the parametric equations of other cross sections can be obtained, and some representative curves of profiles are drawn in Fig. 11 with dashed lines. These profiles, which are theoretically predicted, are in good agreement with those of experiments. While the cross section at the clamped end remains perpendicular to each layer, the cross section at the free end, whose profile corresponds to the dashed curve in the northwest part of Fig. 11, is no longer a plane and not perpendicular to each layer.

We derive Eq. (A1) solely from the geometry, and the properties of the material, including interlayer properties, do not enter into the discussion. Therefore, Eq. (A1) is valid irrespective of the properties of material, which we can use to consider the bending behavior with finite interlayer shear effect.

## References

_{2}