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.

Whereas 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 graphene-based 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 × L, fold width w, and density of adatoms ρ (ratio of the numbers of hydrogen and carbon atoms in the red and blue areas in Figure a). 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 × L, w, and/or ρ.

Figure 1

(a) Formation of a GMS from a graphene sheet. Hydrogen atoms are located on the top side of the sheet in the areas colored in red and on the bottom side of the sheet in the areas colored in blue. (b) Effective lengths L in the x-direction (normalized with respect to the length at 300 K) as functions of temperature for GMSs with w = 3.7 nm and ρ = 15%; L = L = 25 nm for GMS #1, 31 nm for GMS #2, and 37 nm for GMS #3. The black solid, dashed, and dash-dotted lines are linear fits for GMS #1, #2, and #3, respectively.

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 a 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 b presents effective lengths L 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 = L = 25 nm (GMS #1), 31 nm (GMS #2), and 37 nm (GMS #3). As we find for all cases L ∝ aT + b, the coefficient of thermal expansion in the x-direction at 300 K is approximated as α = 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–1[18]). 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 y-direction 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 = L 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 a). The change of L with temperature is controlled by two deformation modes, the dilation and flapping of connected graphene pieces (Figure a). 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 b); 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.

Figure 2

(a) Dilation and flapping of connected graphene pieces in a GMS. (b) Coefficient of thermal expansion in the x-direction as a function of the total density of hydrogenation: GMSs #1, #2, and #3 compared to a graphene sheet.

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 fold of width w. The bending moment density exerted on the thin plate is M = hsρ/2, where h is the effective thickness and D the bending stiffness of graphene. Kirchhoff–Love plate theory then impliesTo 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 a) and D increases significantly (Figure b; consistent with previous simulations[35,36]); that is, eq implies that θ decreases (compare the inset of Figure b). 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.

Figure 3

Mechanism of enhancement of the coefficient of thermal expansion (“positive mechanism”). (a) Pseudosurface stress induced by hydrogenation of graphene as a function of temperature. (b) Bending stiffness of graphene as a function of temperature. The results are obtained for a 22 nm × 22 nm graphene sheet. Inset: Schematic of the decreasing folding angle. As the temperature increases, the pseudosurface stress decreases and the bending stiffness increases, causing the fold angle to decrease, which results in a positive coefficient of thermal expansion.

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 (Figure a). 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 b. 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 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 is larger for GMSs with smaller Young’s modulus; that is, the negative mechanism is stronger. We fit each curve in Figure b as in the strain interval from −3 to +3% (where E and ε are Young’s modulus and the strain in the x-direction, respectively), yielding E = 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.

Figure 4

Mechanism of reduction of the coefficient of thermal expansion (“negative mechanism”). (a) Schematic of the change in the folding angle due to out-of-plane thermal fluctuations. The dashed lines indicate the original effective length of the GMS. (b) Energy density versus strain curves for the three GMSs of Figure b. For the same energy density, the compressive strain is larger than the tensile strain (asymmetric curves); that is, the GMSs are stiffer under tension than under compression (stiffness asymmetry).

We now focus on the effects of L = L, w, and ρ on the coefficient of thermal expansion. The results in Figure indicate a decrease for increasing L = L, consistent with Figure b. Note that for GMSs with the same w and ρ the contribution of the positive mechanism is the same regardless of L = L, as can be concluded from eq . On the other hand, the contribution of the negative mechanism is larger for GMSs with larger L = L due to the decreasing Young’s modulus, which explains the observed trend. Very large negative coefficients of thermal expansion can be achieved (see Figure ). We find for L = L = 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 Figure a the effective length of the unit cell in the x-direction and from ref (37) we obtain θ = cos–1(2sin2(ϕ/2)/sin2β – 1). As L = L implies β = π/4, we haveandWhen the positive mechanism dominates, we can use eqs , 2, and 3 to rewriteaswhich increases for increasing w, ρ, and θ. Therefore, for GMSs with larger w and ρ (larger θ, see eq ), the coefficient of thermal expansion is larger. For L = L = 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 = L = 37 nm, we find that the positive mechanism dominates for small w and ρ, whereas the negative mechanism dominates for large w and ρ.

Figure 5

Effect of L = L, w, and ρ on the coefficient of thermal expansion of GMSs at 300 K.

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 α = (−220 ± 28) × 10–6 K–1 in the case of single-side surface functionalization for L = L = 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 a after confirming that for a GMS with L = L = 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 = L = 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 ±108 s–1 using a NPT ensemble in that the stress component in the y-direction is controlled to be zero during the loading process.