Graphene Origami with Highly Tunable Coefficient of Thermal Expansion

The coefficient of thermal expansion, which measures the change in length, area, or volume of a material upon heating, is a fundamental parameter with great relevance for many applications. Although there are various routes to design materials with targeted coefficient of thermal expansion at the macroscale, no approaches exist to achieve a wide range of values in graphene-based structures. Here, we use molecular dynamics simulations to show that graphene origami structures obtained through pattern-based surface functionalization provide tunable coefficients of thermal expansion from large negative to large positive. We show that the mechanisms giving rise to this property are exclusive to graphene origami structures, emerging from a combination of surface functionalization, large out-of-plane thermal fluctuations, and the three-dimensional geometry of origami structures.

W hereas most materials expand upon heating due to a positive coefficient of thermal expansion, graphene can have, as another interesting property of this famous material, a negative coefficient of thermal expansion 1 owing to out-of-plane thermal fluctuations. 2,3 However, this does not provide tunability of the coefficient of thermal expansion to a target value (negative or positive). On the other hand, this tunability is critical for many applications including construction, 4 fuel cells, 5 and electronic devices 6 because mismatch between the coefficients of thermal expansion of constituent materials can cause large thermal stresses, leading to reduced reliability 7 or even fracture. 8 There are three main approaches to control the coefficient of thermal expansion of a material. 9−11 The first approach exploits supramolecular mechanisms such as transverse vibration, 12 the rigid-unit model, 13 and phase transformation 14 to obtain a negative coefficient of thermal expansion. The second approach introduces multimaterial lattices to achieve a negative or nearly zero effective coefficient of thermal expansion, originating from different thermal deformations of the constituents (due to different coefficients of thermal expansion) and the geometrical constraints of the multimaterial. 15,16 The third approach consists of designing composites consisting of materials with negative (using the first approach) and positive coefficients of thermal expansion. 17 There have been a few works tailoring the coefficient of thermal expansion of graphene by crystal defects 18 and doping, 19 but the thermal expansion is only tunable within a narrow range. Therefore, approaches that can provide coefficients of thermal expansion in a wide range, including not only negative but also zero and positive values, must be developed to broaden thermal applications of graphene-based materials.
The origami approach, or the art of paper folding, which originated in East Asia, is used to generate 3D structures with extraordinary properties such as bistability in a square twist origami 20 and negative Poisson's ratio in the well-known Miura origami. 21 It has been applied in various fields from solar cells 22 and robotics 23 to biomedical applications 24 and electronic devices. 25 In particular, the Miura origami can be used as a platform for controlling the coefficient of thermal expansion by a mechanism similar to the multimaterial lattices mentioned above. 26 However, this approach requires special arrangement of the faces of materials with different stiffnesses, which is impossible when graphene is involved. In addition, despite significant progress in 2D material synthesis, most 3D graphenebased structures are cellular, such as foams, sponges, and aerogels (see ref 27 and references therein). Only simple origami structures (basic polyhedra) have been obtained so far, 28 whereas the realization of complex 3D origami structures, such as the graphene Miura origami structure, is still a great challenge.

RESULTS AND DISCUSSION
In this study, we use molecular dynamics (MD) simulations to demonstrate that graphene Miura origami structures (GMSs) can exhibit large and tunable coefficients of thermal expansion, both positive and negative. In the temperature range from 250 to 350 K, a positive coefficient of thermal expansion of (33 ± 1) × 10 −6 K −1 is achieved, outperforming aluminum and magnesium alloys, 29 and a negative coefficient of thermal expansion of (−465 ± 28) × 10 −6 K −1 , which is significantly more negative than previously reported for any noncellular material. 30,31 Our study leverages recent work using surface functionalization to resolve difficulties in controlling the folding angle and direction of GMSs, which was a principal bottleneck for forming complex 3D origami structures. Specifically, the surface is functionalized in predefined areas of a graphene sheet with a pattern characterized by the unit cell size L x,0 × L y,0 , fold width w, and density of adatoms ρ (ratio of the numbers of hydrogen and carbon atoms in the red and blue areas in Figure 1a). The pseudosurface stress induced by the functionalization drives folding of the graphene sheet at the functionalized areas to form a GMS. We show that the coefficient of thermal expansion of such a GMS is determined by two simultaneous mechanisms: on the one hand, an increasing bending stiffness of graphene for increasing temperature (due to out-of-plane thermal fluctuations) reduces the folding angle induced by the pseudosurface stress to yield a positive coefficient of thermal expansion. On the other hand, interplay of the out-of-plane thermal fluctuations with asymmetry of the in-plane stiffness of the GMS increases the folding angle to yield a negative coefficient of thermal expansion. It turns out that a wide range of coefficients of thermal expansion, from large negative to large positive, can be obtained by changing L x,0 × L y,0 , w, and/or ρ.
We employ MD simulations to simulate the formation of GMSs due to pattern-based surface functionalization and to calculate the coefficients of thermal expansion. Figure 1a shows a schematic of the surface functionalization approach. Carbon atoms are randomly selected on top and bottom of the graphene sheet in the red and blue areas, respectively, with density ρ, and hydrogen atoms then are placed at the top sites of these carbon atoms at a distance of 1.1 Å. Periodic boundary conditions are applied in the zigzag (x) and armchair (y) directions to eliminate edge effects. The interaction between the atoms is modeled by the second generation REBO (REBO-II) potential. 32 Initially, a molecular statics simulation is conducted using the conjugate gradient method. Then, at 300 K, the system is equilibrated using a canonical (NVT) ensemble for 100 ps and afterward fully relaxed (until the potential energy converges) using an isothermal−isobaric (NPT) ensemble with the stress components along the in-plane directions controlled to be zero ( Figure  S1 in the Supporting Information). After equilibration, the GMS is either heated to 350 K with a heating rate of 25 K/ns or cooled to 250 K with a cooling rate of −25 K/ns using a NPT ensemble. Figure 1b presents effective lengths L x in the x-direction (normalized with respect to the length at 300 K) as functions of temperature (T) for GMSs with the same w = 3.7 nm and ρ = 15% but different L x,0 = L y,0 = 25 nm (GMS #1), 31 nm (GMS #2), and 37 nm (GMS #3). As we find for all cases L x ∝ aT + b, the coefficient of thermal expansion in the x-direction at 300 K is approximated as α x = a. GMS #1 expands upon heating with a relatively large positive value of (33 ± 2) × 10 −6 K −1 , clearly exceeding our result of (−1 ± 1) × 10 −6 K −1 for pristine graphene (which agrees with a previous MD simulation 2 and is close to the experimental value of (−3.5 ± 0.5) × 10 −6 K −118 ). In contrast, GMSs #2 and #3 exhibit large negative values of (−15 ± 4) × 10 −6 and (−86 ± 6) × 10 −6 K −1 , respectively. Similarly, we obtain for the coefficient of thermal expansion in the ydirection values of (12 ± 1) × 10 −6 K −1 (GMS #1), (−2 ± 1) × 10 −6 K −1 (GMS #2), and (−15 ± 2) × 10 −6 K−1 (GMS #3). Therefore, modifying L x,0 = L y,0 enables us to realize both in the x-and y-directions coefficients of thermal expansion in a wide range from large negative to large positive.
A GMS can be regarded as a set of graphene pieces connected at folds characterized by the folding angle θ, which is zero when the structure is flat (Figure 2a). The change of L x with temperature is controlled by two deformation modes, the dilation and flapping of connected graphene pieces (Figure 2a). To separate the two effects by excluding flapping, we study the coefficient of thermal expansion for hydrogen atoms randomly distributed on both sides of a graphene sheet (size 50 nm × 50 nm with periodic boundary conditions), using total densities of hydrogenation (ratio of the numbers of hydrogen and carbon atoms) from 0 to 10%. The value increases only slightly when the total density of hydrogenation increases ( Figure 2b); that is, dilation cannot explain the results for the GMSs. Consequently, flapping (change of θ) is the main cause for the observed variation in the coefficient of thermal expansion between the three GMSs. We will address the mechanisms controlling the flapping in the following.
To describe the folding angle θ caused by hydrogenation, we consider a thin plate under bending due to a pseudo surface stress ρs (where s is the pseudo surface stress at ρ = 100% hydrogenation) and a folding pattern consisting only of a single To understand thermal effects on θ through s and D, we derive these quantities at different temperatures by MD simulations. We obtain D following ref 33. To obtain s, we consider a graphene sheet with 50% hydrogenation on each side (graphane chair 34 ), that is, ρ = 100%, as in this case there is pure stretching (no bending), so that s = hσ can be determined from the in-plane normal stress σ accessible by MD simulation of pristine graphene stretched biaxially with the same in-plane normal strain. For increasing temperature, s decreases slightly ( Figure  3a) and D increases significantly (Figure 3b; consistent with previous simulations 35,36 ); that is, eq 1 implies that θ decreases (compare the inset of Figure 3b). As the change of D is more significant than that of s (for example, s decreases by 0.8% and D inceases by 3.5% from 300 to 350 K), the increasing bending stiffness of graphene upon heating is the main mechanism behind the increasing coefficient of thermal expansion of the GMSs. This mechanism, which we call the "positive mechanism" in the following, implies that the coefficient of thermal expansion of any GMS is always larger than that of pristine graphene. Specifically, GMS #1 has a positive coefficient of thermal expansion, though that of pristine graphene is negative. However, the positive mechanism cannot explain the large negative values obtained for GMSs #2 and #3.
The second mechanism is based on the large out-of-plane thermal fluctuations of graphene combined with the specific geometry of the Miura origami. The folds of a GMS act as mechanical constraints against these fluctuations, which grow for increasing temperature and serve as an external force to bend the graphene pieces, resulting in a change of θ and, thus, of L x (Figure 4a). To better understand the deformations of the three GMSs, we plot their energy density versus strain curves under uniaxial stress in the x-direction at 1 K in Figure 4b. All curves turn out to be asymmetric, with a smaller slope under compression than under tension; that is, the GMSs are less stiff under compression than under tension. The increasing energy of the out-of-plane thermal fluctuations for increasing temperature results in a reduction of L x due to the asymmetry of the energy density versus strain curves, corresponding to a negative coefficient of thermal expansion. In the following, we call this the "negative mechanism". The reduction of L x is larger  for GMSs with smaller Young's modulus; that is, the negative mechanism is stronger. We fit each curve in Figure 4b as u E xx xx in the strain interval from −3 to +3% (where E xx and ε xx are Young's modulus and the strain in the x-direction, respectively), yielding E xx = 157, 76, and 39 MPa for GMS #1, #2, and #3, respectively, that is, a growing contribution of the negative mechanism. As the contribution of the positive mechanism is the same for the three GMSs (same w and ρ), the small contribution of the negative mechanism for GMS #1 explains its large positive coefficient of thermal expansion. Correspondingly, growing contributions of the negative mechanism explain the increasingly negative coefficients of thermal expansion of GMSs #2 and #3.
We now focus on the effects of L x,0 = L y,0 , w, and ρ on the coefficient of thermal expansion. The results in Figure 5 indicate a decrease for increasing L x,0 = L y,0 , consistent with Figure 1b. Note that for GMSs with the same w and ρ the contribution of the positive mechanism is the same regardless of L x,0 = L y,0 , as can be concluded from eq 1. On the other hand, the contribution of the negative mechanism is larger for GMSs with larger L x,0 = L y,0 due to the decreasing Young's modulus, which explains the observed trend. Very large negative coefficients of thermal expansion can be achieved (see Figure 5). We find for L x,0 = L y,0 = 25 nm (positive mechanism dominates) an increase of the coefficient of thermal expansion for increasing w and ρ, which can be explained as follows. We extract from When the positive mechanism dominates, we can use eqs 1, 2, and 3 to rewrite which increases for increasing w, ρ, and θ. Therefore, for GMSs with larger w and ρ (larger θ, see eq 1), the coefficient of thermal expansion is larger. For L x,0 = L y,0 = 50 nm (negative mechanism dominates), increasing w and ρ lead to a more negative coefficient of thermal expansion due to a smaller in-plane stiffness ( Figure S2 in the Supporting Information). Moreover, for intermediate L x,0 = L y,0 = 37 nm, we find that the positive  ACS Nano www.acsnano.org Article mechanism dominates for small w and ρ, whereas the negative mechanism dominates for large w and ρ. Finally, we note that material properties similar to those described above can be achieved by single-side surface functionalization or/and larger fold widths. For example, we obtain α x = (−220 ± 28) × 10 −6 K −1 in the case of single-side surface functionalization for L x,0 = L y,0 = 62 nm, w = 12.3 nm, and ρ = 15% (see Figure S3 in the Supporting Information).

CONCLUSIONS
We use MD simulations to demonstrate that GMSs obtained by hydrogenation can exhibit a wide range of coefficients of thermal expansion from large negative to large positive by controlling the folding pattern, in particular, in terms of the unit cell size, the folding width, and the density of hydrogenation. This property turns out to be a consequence of the interplay between two mechanisms (the positive and negative mechanisms) that we explain by microscopic considerations. Whereas many unusual properties of conventional metamaterials arise from specific geometrical features, 21,38 the wide range of coefficients of thermal expansion accessible in GMSs turns out to be the result of an inseparable combination of surface functionalization, large out-of-plane thermal fluctuations (intrinsic property of graphene), and the Miura origami geometry (extrinsic property). First, surface functionalization is the basis for forming the GMSs and provides the pseudosurface stress for the positive mechanism. Second, the large out-of-plane thermal fluctuations due to the single-atomic thickness of graphene contribute by increasing the bending stiffness for increasing temperature (in the case of the positive mechanism) and by causing flapping (in the case of the negative mechanism). The fact that the pseudosurface stress and thermal fluctuations cannot be effective mechanisms in three-dimensional materials and are unlikely to be effective in other two-dimensional materials highlights the importance of graphene in enabling the demonstrated wide range of coefficients of thermal expansion. Third, the Miura origami geometry has the key role to provide a three-dimensional structure. Our work provides fundamental insights into mechanisms that can be utilized to engineer the thermal properties of 2D materials.

METHODS
We use the open-source LAMMPS code 39 to perform molecular statics and MD simulations and the OVITO software 40 for visualizing the simulation results. The zigzag and armchair directions of graphene correspond to the x-and y-directions, respectively. We consider the 2 × 2 supercell shown in Figure 1a after confirming that for a GMS with L x,0 = L y,0 = 25 nm, w = 2.5 nm, and ρ = 15% the difference in the obtained coefficients of thermal expansion with respect to a 4 × 4 supercell is negligible. An initial relaxation is performed by a molecular statics simulation with fixed size of the supercell. The conjugate gradient minimization is deemed to be converged when the relative change of the total energy in successive iterations falls below 10 −16 . Then a NVT ensemble is used for relaxation at 300 K for 100 ps. All other MD simulations (equilibration as well as heating and cooling) are executed in a NPT ensemble to ensure zero in-plane stress components. A time step of 1 fs is chosen. A heating/cooling rate of ±25 K/nm is used after confirming for a GMS with L x,0 = L y,0 = 25 nm, w = 2.5 nm, and ρ = 15% that the results do not change relevantly as compared to a heating/ cooling rate of ±2.5 K/ns. We assign the initial velocities of the atoms five times randomly and consider five random distributions of hydrogen atoms for each GMS and the graphene sheet to estimate the statistical error by means of the corrected sample standard deviation. To obtain the energy density versus strain curve, the GMS is initially fully relaxed at 1 K using the procedure mentioned above. It is then stretched/ compressed in the x-direction with a strain rate of ±10 8 s −1 using a NPT ensemble in that the stress component in the y-direction is controlled to be zero during the loading process.