Nonequilibrium Thermodynamics of Acoustic Phonons in Suspended Graphene

Recent theory has predicted large temperature differences between the in-plane (LA and TA) and out-of-plane (ZA) acoustic phonon baths in locally-heated suspended graphene. To verify these predictions, and their implications for understanding the nonequilibrium thermodynamics of 2D materials, experimental techniques are needed. Here, we present a method to determine the acoustic phonon bath temperatures from the frequency-dependent mechanical response of suspended graphene to a power modulated laser. The mechanical motion reveals two counteracting contributions to the thermal expansion force, that are attributed to fast positive thermal expansion by the in-plane phonons and slower negative thermal expansion by the out-of-plane phonons. The magnitude of the two forces reveals that the in-plane and flexural acoustic phonons are at very different temperatures in the steady-state, with typically observed values of the ratio $\Delta T_{\mathrm{LA+TA}}/\Delta T_{\mathrm{ZA}}$ between 0.2 and 3.7. These deviations from the generally used local thermal equilibrium assumption ($\Delta T_{\mathrm{LA+TA}}=\Delta T_{\mathrm{ZA}}$) can affect the experimental analysis of thermal properties of 2D materials.

The thermal properties of graphene [1] are unconventional, because of the large difference between its inplane and out-of-plane lattice dynamics [2,3]. Therefore, much research has focused on characterizing graphene's thermal conductivity, for example by using Raman spectroscopy or electrical heaters [4][5][6][7][8][9][10][11]. Recent theoretical work by Vallabhaneni et al. has suggested that local optical heating of suspended graphene can lead to a large temperature differences between the in-plane (longitudinal LA and transverse TA) and out-of-plane (flexural, ZA) acoustic phonon baths, which is caused by differences in the thermal conductivities of the different types of phonons, and their weak mutual interactions [12]. It has been confirmed experimentally that electrons and optical phonons can show very different temperatures compared to the acoustic phonons in 2D materials [13][14][15][16][17][18], but whether strong thermal nonequilibrium between the acoustic phonon modes themselves exists has not been established. Since it has been hypothesized that such a thermal nonequilibrium might impact the interpretation of the widely used Raman spectroscopy technique to measure the thermal conductivity of graphene [12], there is a need to characterize the temperatures of the in-plane and flexural acoustic phonon baths separately.
Recently, several optomechanical techniques to characterize the time-dependent heat transport in suspended 2D materials have been developed [19][20][21][22]. Here, we demonstrate the use of an optomechanical technique to distinguish two thermal expansion force contributions with different time-constants and opposite signs. It is argued that these contributions can be attributed to the in-plane and flexural acoustic phonons. The differences in time-constant and sign allow us to obtain information on the modal temperatures of the respective phonon baths. thermally actuate and measure the motion of suspended graphene membranes. The sample fabrication is identical to that in earlier work [19]. Single-layer graphene grown by chemical vapor deposition (CVD) is transferred over dumbbell-shaped cavities in a Si/SiO 2 substrate (300 nm deep, various diameters) using a support polymer. The polymer is dissolved and the sample is dried using critical point drying, which breaks one of the dumbbell drums while the other side survives, resulting in a circular graphene drum with a venting channel to the environment that prevents gas from being trapped between the membrane and the substrate.
To actuate the motion, the membrane is heated in a vacuum (pressure lower than 1 × 10 −5 mbar) by a sinusoidally-power-modulated blue laser. The blue laser, which is focused at the center of the drum has a wavelength of 405 nm, an average incident laser power of 0.36 W and its sinusoidal modulation amplitude is 0.24 W. Electrons in the graphene drum are photoexcited, and decay into thermal phonons in less than a picosecond [12,23,24] via electron-phonon scattering. Compared to the timescales at which phonons exchange heat (> 0.1 ns), the power transfer from light to lattice vibrations via electron-phonon scattering can thus be considered instantaneous. The out-of-plane membrane motion is read out using a photodiode (PD) that detects the reflected intensity of a 633 nm red helium-neon laser with a power of 1.2 mW focused on the center of the membrane, that is modulated by the position-dependent absorption of the graphene membrane [25,26]. The estimated waist diameter of the focal point is 0.67 µm for the red laser and 0.57 µm for the blue laser [20]. A vector network analyzer (VNA) measures the frequency-dependent amplitude and phase of the signal at the output of a photodetector relative to the modulated blue laser power. The signal is corrected for parasitic phase shifts due to delays in the optical and electronic path using a calibration measurement [19], which ensures that the voltage change from the photodiode is linearly proportional to the deflection of the membrane. All experiments are performed at room temperature.
Figures 2(a) and (b) show the real and imaginary amplitude of the membrane's motion as a function of frequency. To analyse the data, the membrane temperature response ∆T to a modulated input power P ac is modeled by the heat equation C∆Ṫ + ∆T /R = P ac e iωt , where C is the effective heat capacitance and R is the effective thermal resistance of the membrane. ∆T is the average temperature change over the suspended drum area with respect to the environmental temperature T 0 , such that the total temperature is given by: T = T 0 + ∆T . The thermal expansion force is assumed to be proportional to the change in temperature ∆T (t): F (t) = α eff ∆T (t). Furthermore, we assume that far below the mechanical resonance frequency the displacement amplitude z(t) = F (t)/k, where k is the effective membrane stiffness. The time-dependent thermal expansion force represented in the frequency domain is then [19,20]: where C slow is a constant, representative of the amplitude of the thermal expansion force at low frequencies, used for fitting. ω is the driving frequency, τ 1 = RC the thermal time constant and F ω is obtained by the Fourier transform of F (t). The imaginary part of Eq. 1 has an extremum with amplitude C slow /2 at radial frequency ω = 1/τ slow as indicated in Figs. 2(a) and (b). Only the imaginary part of Eq. 1 is fit to the data, showing good agreement with the experimentally obtained imaginary amplitude. If the real part corresponding to this fit is plotted, however, it is found that below the resonance frequency there is an additional offset C fast , between the real part of Eq. 1 and the measurement, that is almost frequency independent, as indicated in Figs. 2(a) and (b). To quantify the value of C fast experimentally, the average value of the difference between the real part of the model and the experimental data (see Figs. 2(a) and (b)) at frequencies below the resonance frequency is taken. All drums with a negative value of C slow have a positive offset in the real part C fast and drums with a positive C slow have a negative C fast (Fig. 2(c)). We deduce from this correlation between C slow and C fast that the offset C fast is not due to optical cross-talk from the blue laser [19], but related to the membrane motion, because optical cross-talk in the setup is independent of the motion of the drum and independent of the sign of C slow .
Since the motion that corresponds to the offset C fast cannot be accounted for by the force in Eq. 1, it is interpreted as evidence for the existence of an additional second contribution to the thermal expansion force, with a different time-constant τ fast . This results in a modified expression for the total thermal expansion force F ω : If ωτ fast 1, the second contribution to the actuation force produces a constant offset in the real part and does not affect the imaginary part of F ω . Furthermore, a key finding of this work is that C fast always has an opposite sign with respect to C slow (Fig. 2(c)), meaning that both these forces are counteracting at low frequencies.
To show this effect more clearly, the frequency domain response of Eq. 2 is converted to a step response in the time domain in Fig. 3, using typical values of τ slow found in this work and an estimate of τ fast based on theory [27]. Our measurement thus indicates that when a constant heating power is suddenly applied at t = 0, the membrane initially expands at short timescales τ fast and then slowly contracts at longer timescales τ slow . To explain this observation, the microscopic origin of the thermal expansion contributions in graphene is analyzed in more detail.
The force that actuates the membrane F (t) is directly proportional to the in-plane stress σ(t), which is linearly related to the in-plane strain (t) by the elasticity matrix. For a membrane clamped around its circumference, this thermally induced strain is related to the internal energies of the phonons and the mechanical properties of the lattice by the equation [28]: where 0 is the initial pre-strain at the reference temperature when ∆T = 0 K, B the bulk modulus, γ j the modedependent Grüneisen parameter and U j (t) the phonon energy per unit volume for phonon mode i. Note that (t) is the total strain with respect to the initial positive (tensile) pre-strain 0 at a reference temperature for which U j = 0, which is reduced by thermal expansion of the membrane. Thermal expansion of the substrate is neglected in this analysis, since it absorbs less laser power than the graphene and because the volume where the heat can diffuse through is much larger, resulting in negligible temperature changes of the substrate. Only the contributions of the acoustic phonon modes are included in the following analysis since the in-plane optical phonon states are not occupied at room temperature and the flexural optical phonons have a Grüneisen parameter close to zero [29]. It is well known that the Grüneisen parameter for the flexural phonons γ ZA has a negative sign in graphene, while the Grüneisen parameter for the inplane longitudinal acoustic (γ LA ) and transverse acoustic mode (γ TA ) is positive [28]. At low laser modulation frequencies, the internal energy U j (t) of all phonon modes is in-phase with the blue laser, such that the sign of the thermal expansion force only depends on the sign of the Grüneisen parameter. Based on these considerations, the most likely conclusion is that the opposite signs of C slow and C fast in the experiments in Fig. 2 can be attributed to the opposite signs of the in-plane and out-of-plane phonon mode Grüneisen parameters [28]. We furthermore hypothesize that the flexural ZA phonons have a longer thermal timescale τ slow , because they experience a large thermal interface resistance at the edge of the drum [27], while the fast timescale τ fast is attributed to the in-plane phonons. The theory that theoretically supports the correctness of this hypothesis is presented in a separate article [27]. The average internal energy U j of the suspended graphene is modulated by the blue laser with an amplitude that depends on the heat flux absorbed by each mode P j , the mode's Grüneisen parameter γ j and its thermal time constant τ j . We find expressions for the average internal energies U j in the Supplemental Information [30] and substitute these in Eq. 7 to obtain: where τ LA+TA is the fast time constant associated with the in-plane phonons (attributed to τ fast for both LA and TA phonons) and τ ZA is the slow time constant from the flexural phonons (attributed to τ slow ). Furthermore, an analytical expression for the time constants corresponding to the flexural phonons τ ZA is derived by taking only the interaction between the phonon modes at the boundary into account. Using the model in Ref. [27], it is found that the time constant τ ZA can be approximated by the expression: where a is the radius of the drum, w 1z→2r the fraction of ZA phonons that transmit over the boundary towards the environment and c ZA is the ZA phonon propagation velocity. Both w 1z→2r and c ZA are tension-dependent parameters, which increase their value with increasing tension.
Based on Eq. 14, a linear relation between the parameters −C fast /C slow and 1/τ ZA is expected, assuming that  [30]. The second horizontal axis shows the temperature ratio calculated from Eq. 6.
τ LA+TA and P LA+TA are constant because they are relatively unsensitive to tension variations. A study of their correlations is therefore useful as a test for the hypothesis behind Eq. 14, since τ ZA and −C fast /C slow can be extracted independently from the measurement. Figure  5(a) shows a plot with −C fast /C slow on the horizontal axis and 1/τ slow on the vertical axis, for drums with a diameter of 6 µm. We find a significant linear correlation between a/τ slow and −C fast /C slow , which is in agreement with the models underlying Eqs. 14. Interestingly, in the Supplemental information [30], we find that this correlation is diameter dependent.
The ratio −C fast /C slow can be used to estimate the degree of thermal nonequilibrium in the system. It is assumed that the changes in modal temperatures of the in-plane phonons are equal ∆T LA+TA = ∆T LA = ∆T TA , based on the results obtained in Ref. 12. Combining Eq. 14 with the thermal expansion term C j = γ j ρc p,j ∆T j , where c p,j is the modal specific heat at constant pressure and ρ is the density, we obtain using Eq. 14: The temperature ratio at low frequencies is thus proportional to −C fast /C slow with a proportionality constant that can be evaluated from theory. Using γ LA = 1.06, γ TA = 0.40 and γ ZA = −4.17 [29], and the modal specific heats (c p,LA = 104 J/(kg·K), c p,TA = 225 J/(kg·K), c p,ZA = 358 J/(kg·K)) calculated at an environmental temperature of 293.15 K, we obtain: ∆T LA+TA /∆T ZA = −7.45C fast /C slow . Using this expression, a histogram of the temperature ratio is constructed as shown in Fig.  5(b). The average value of the ratio ∆T LA+TA /∆T ZA is of the order of 1 (see Fig. 5 and the Supplemental information [30]). This is surprising, because the observation that τ LA+TA τ ZA suggests that the ZA phonons have a very low thermal conductance and therefore, according to Eqs. 14 and 6, we should expect ∆T ZA ∆T LA+TA . This apparent contradiction between the observed thermal time constants and the temperature ratio is explained the selective electron-phonon coupling in graphene, which causes most of the heat supplied by the laser to end up in the LA and TA phonon bath, while the ZA phonons only receive a small fraction of this heat due to the weak coupling [12]. The small value of τ LA+TA /τ ZA in Eq. 14 is thus partially compensated by the large value of P LA+TA /P ZA , thereby causing the temperature of the in-plane and flexural acoustic phonon bath to be in the same order of magnitude.
In Fig. 5, although a few drums have almost the same value for in-plane and out-of-plane temperature, in many drums large variations in the temperature ratio are observed, with ∆T LA+TA /∆T ZA varying from 0.2 to 2.2. This provides evidence for the existence of a strong nonequilibrium thermal state. According to Eq. 5, τ ZA is tension-dependent, whilst τ LA+TA is not expected to be tension dependent. Consequently, according to Eq. 14, a linear correlation between τ slow and |C fast /C slow | as found in Fig. 5(a) shows that the large variations in the temperature ratio are dominated by device-to-device variations in the pre-tension via its effect on the thermal time constant τ ZA . Similar large variations in τ ZA have been observed in our previous work [19]. Some devices deviate from this linear correlation (see Supplemental Information [30]), this might suggest that other effects such as wrinkles and other imperfections are also playing a role in the variations in the temperature ratio.
The observed nonequilibrium effect has important consequences for the interpretation of thermal measurements on graphene, as it becomes difficult to determine the contribution of each phonon mode to the thermal conductivity. Moreover, the extracted thermal conductivity obtained from the classical heat equation can become geometry dependent. This affects any suspended graphene device that is locally heated, due to the inherent selective electron-phonon coupling and weak interaction between the phonon modes.
Other two-dimensional materials are expected to show similar effects as observed in this work if they exhibit weak mode interaction and a large negative Grüneisen parameter for the flexural acoustic phonons. This might hold for other monatomic two-dimensional materials at room temperature [28], but also transition metal dichalcogenides such as MoS 2 , MoSe 2 and WS 2 at low temperatures (< 100 K) [31].
To conclude, we have presented evidence that the motion of opto-thermally excited graphene resonators is the result of two counteracting contributions in the thermal expansion force. The amplitude of these contributions provides information on the ratio of the effective temperatures of the thermal baths of the in-plane and flexural acoustic phonons and based on a model it is shown that they are at different temperatures. These thermal nonequilibrium effects should be considered in the interpretation of the thermal conductivity measurements of 2D materials. Moreover, they are shown to lead to an unconventional time-dependent sign of thermal expansion forces in graphene.
The authors thank Applied Nanolayers B.V. for supply and transfer of the single-layer graphene. We furthermore thank C. Stampfer, J. Sonntag and J.E. Sader for fruitful discussions. This work is part of the research programme Integrated Graphene Pressure Sensors (IGPS) with project number 13307 which is financed by the Netherlands Organisation for Scientific Research (NWO). The research leading to these results also received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 785219 Graphene Flagship. Figure 5(a) shows the correlation between the parameters a/τ slow and −C fast /C slow . All diameters show a significant linear correlation between these parameters. However, the slopes of the straight lines that are fit to the data are diameter dependent. To investigate this diameter dependence further, a boxplot of −C fast /C slow and ∆T LA+TA /∆T ZA for each diameter is made as shown in Fig. 5(b), and a box plot of a/τ slow is made. Comparing the diameter dependence of −C fast /C slow to the diameter dependence of a/τ slow , the largest relative change in the mean values is observed in −C fast /C slow , suggesting this is the underlying mechanism behind the diameter dependent slope in Fig. 5(a).
Our model in Ref. [27] suggests no diameter-dependence of a/τ slow should occur, if the tension in the drums is not diameter dependent. No significant diameter-dependence of a/τ slow is discovered in Fig. 5(b), except for the 7-micron diameter drums where small values of τ slow appear to be missing. Drums with a small value of τ slow are expected to have a large tension. Therefore, this might be an indication of a selective breaking mechanism for these drums, where drums with large tensions are more likely to fail. Apart from the 7-micron drums, each diameter shows a considerable spread in a/τ slow , which may be attributed to device-to-device variations in the tension of the drums, that alter τ slow . This is also the underlying reason for spread along the straight lines in Fig. 5(a).
The diameter dependence of −C fast /C slow is unexpected, as our model [27] predicts that |C fast /C slow | increases as a function of diameter, while our data in Fig. 5(b) shows that it decreases. Several scenarios are investigated in Ref. [27] in order to explain this effect. The most probable explanation is that the diameter dependence is caused by ballistic effects in phonon transport.

S2: DERIVATION OF THE EXPRESSION FOR THE RATIO C fast /C slow
Here the expression for the ratio between the two contributions to the thermal expansion force C fast /C slow is derived. We assume for each force that ω τ i . From eq. (3) in the main section of the paper we have: Linearizing the system for small ∆T : Since the thermal expansion forces C slow and C fast are proportional to the respective phonon bath contributions to the thermal strain , it follows from Eqs. 7 and 8 that: C fast C slow = γ LA ρc p,LA ∆T LA + γ TA ρc p,TA ∆T TA γ ZA ρc p,ZA ∆T ZA Since ∆T j = R B,j P ac,j /A, where A is the area of the circumference and R B,j the thermal interface resistance, we can write: C fast C slow = γ LA ρc p,LA R B,LA P ac,LA + γ TA ρc p,TA R B,TA P ac,TA γ ZA ρc p,ZA R B,ZA P ac,ZA , Now it is convenient to convert the specific heat ρc p,j into the modal heat capacity C j , from [19]: For the thermal resistance: Using this and the relation τ j = R j C j , we arrive at: C fast C slow = γ LA τ LA P ac,LA + γ TA τ TA P ac,TA γ ZA τ ZA P ac,ZA , Since it is assumed the in-plane LA and TA phonons are at the same temperature, we can write this expression as: