Optomechanical methodology for characterizing the thermal properties of 2D materials

Heat transport in two-dimensions is fundamentally different from that in three dimensions. As a consequence, the thermal properties of 2D materials are of great interest, both from scientific and application point of view. However, few techniques are available for accurate determination of these properties in ultrathin suspended membranes. Here, we present an optomechanical methodology for extracting the thermal expansion coefficient, specific heat and thermal conductivity of ultrathin membranes made of 2H-TaS2, FePS3, polycrystalline silicon, MoS2 and WSe2. The obtained thermal properties are in good agreement with values reported in the literature for the same materials. Our work provides an optomechanical method for determining thermal properties of ultrathin suspended membranes, that are difficult to measure otherwise. It can does provide a route towards improving our understanding of heat transport in the 2D limit and facilitates engineering of 2D structures with dedicated thermal performance.


I. INTRODUCTION
Soon after the discovery of monolayer graphene, it was found that 2D materials have unique thermal properties, which open opportunities for heat control at the nanoscale [1][2][3][4][5] .Due to their ultrasmall thickness, thermal properties of 2D materials are dominated by surface scattering of acoustic phonons, which is highly sensitive to strain 6 , grain size 7 and temperature 8 , as well as material imperfections such as defects and impurities 9 .To understand and optimize heat transport in 2D materials, precise thermal characterization methods are of great importance.
So far, a variety of experimental techniques have been developed to characterize thermal transport in 2D materials, of which the transient micro-bridge method 10,11 and the steadystate optothermal method based on Raman microscopy are most commonly used 12,13 .However, the construction of a micro-bridge is complicated and thermal contact resistances can affect measurement results, while for Raman measurements, the probed temperature resolution is usually relatively small, leading to large error bars.These limitations undermine the accuracy of probing heat transport in 2D materials, causing large variations in the thermal material parameters reported in literature.For example, literature values for the thermal conductivity vary from 2000 to 5000 Wm −1 K −1 for suspended monolayer graphene 14 .
In this work, we demonstrate an optomechanical noncontact method for measuring the thermal properties of nanomechanical resonators made of free-standing 2D materials.The presented methodology allows us to simultaneously extract the thermal expansion coefficient, the specific heat and the in-plane thermal conductivity of the material.It involves driving a suspended membrane using a power-modulated laser and measuring its time-dependent deflection with a second laser.Thus both the temperature-dependent mechanical fundamental resonance frequency of the membrane and charac-teristic thermal time constant at which the membrane cools down 15 are measured.A major advantage of the method is that no physical contact needs to be made to the membrane, such that its pristine properties are probed and no complex device fabrication is needed.Buckling effects are incorporated in the model to account for the induced compressive stress during temperature variations.Our results on 2H-TaS 2 , FePS 3 , polycrystalline silicon (Poly Si), MoS 2 and WSe 2 show good agreement with reported values in the literature.

II. FABRICATION AND METHODOLOGY
We fabricate 2D nanomechanical resonators by transferring 2D flakes over circular cavities with a depth of 285 nm and a radius R of 3 to 4 µm in a silicon (Si) substrate with a 285 nm thick silicon oxide (SiO 2 ) layer, as illustrated in Fig. 1a.The devices D1−D5 studied in this work are made of 2H-TaS 2 , FePS 3 , Poly Si, MoS 2 and WSe 2 , respectively.By using tapping mode Atomic Force Microscope (AFM), we determine the thickness, h, of each membrane (see Table I).All details about the device fabrication and thickness measurement can be found in SI section 1.To determine the Young's modulus E of each membrane, we use the AFM to indent the centre of suspended area with a force F while measuring the cantilever indentation δ 16 .The measured F versus δ , as depicted in Fig. 1b for device D1, is fitted with a model for point-force loading of a circular plate given by F = ( 16πD is the initial tension in the membrane, and ε 0 is the prestrain.We extract E = 108.45GPa and ε 0 = 6.75 × 10 −3 from the fit shown in Fig. 1b (drawn line), which are in good agreement with typical values found in literature for 2H-TaS 2 17 .The obtained values of E for devices D2−D5 are listed in Table I.
The setup for the optomechanical measurements 23,24 , is   [17][18][19][20][21][22] .shown in Fig. 1c.A power-modulated blue diode laser (λ = 405 nm) photothermally actuates the resonator, while a He-Ne laser (λ = 632 nm), of which the reflected laser power depends on the position of the membrane, is used to detect the motion of the resonator.The power-modulation of the blue laser is supplied by a Vector Network Analyzer (VNA), which also analyzes the photodiode signal containing the reflected laser power and converts that to the response amplitude, |z f |, of the resonator in the frequency domain.All measurements were done in vacuum at a pressure below 10 −5 mbar.As shown in Fig. 1d, |z f | shows a clear fundamental resonance peak, which we fit with a harmonic oscillator model, given , where f 0 is the fundamental resonance frequency, A res is the vibration amplitude at resonance and Q is the quality factor.For device D1, we obtain a f 0 = 9.53 MHz and Q = 160.In addition, we also find a maximum in the imaginary part of z f at kHz frequencies (see Fig. 1e), which we attribute to the thermal expansion of the membrane that is time-delayed with respect to the modulated blue laser power, because it takes a time τ for the temperature of the membrane to rise 25,26 .By solving the in-plane heat equation in the membrane, the thermal signal can be expressed as: where A th and τ are the thermal expansion amplitude and thermal time constant of the membrane, respectively.The red and blue laser powers are fixed at 0.9 and 0.13 mW respectively, to ensure linear vibration of the resonators with a negligible temperature raise of the membrane due to self-heating 25 .

A. Thermally-induced buckling phenomenon
When changing the temperature, the thermal expansion coefficient (TEC) α m of the membrane, which is higher than that of the silicon substrate α Si , changes the strain in the membrane by a quantity ∆ε.This results in a remarkable change in the dynamics of 2D nanomechanical resonators, which can be used for probing the thermal properties 27,28 .Therefore, we heat up the fabricated devices and investigate the dependence of resonance frequency f 0 on temperature T .As shown in Fig. 2a, we observe a decrease of f 0 with increasing T for device D1, which is in agreement with trends shown in literature 29 and can be attributed to a reduction in strain when the material thermally expands.However, the results obtained for devices D2 and D3 are substantially different: as shown in Figs.2b and 2c, we observe an initial decrease in f 0 with increasing T towards a minimum frequency (which we call the turning point), followed by a continuous increase.We attribute this to the thermally-induced buckling of the mechanical resonators as found in earlier studies [30][31][32] , which is caused by the loaded compression since α m > α Si .Here, as depicted in Fig. 2d, the thermal expansion of the membrane causes a compressive stress that triggers the membrane to buckle.We label the pre-buckling, the transition from pre-to postbuckling, and the post-buckling regions in Fig. 2e as cases I, II and III, respectively.
We use a Galerkin model for a clamped circular plate 33,34 , to find an approximate analytical expression of the fundamental resonance frequency f 0 under thermally-induced buckling 35 : where d is the diameter of the plate, U is the thermally changed in-plane displacement from boundary, ρ is the mass density, z is the central deflection of the plate, ν is the Poisson ratio, z f ree is the central deflection of the plate in the predeformed state when U = 0 (without loading), and β is a fitting factor determined by 35 β = 0.35ν + 0.42.Equation (2)  shows that f 0 depends on the in-plane displacement U(T ) from boundary and the central deflection z(T ) of the membrane.The relation between z(T ) and U(T ) can be found in SI section 1 (see Eq. (S1)).Following the literature 35 , z f ree can be extracted from the measured value of the fundamental resonance frequency at the turning point, f t , using the built Galerkin model.By substituting R = 4 µm, h = 23 nm, E = 108 GPa, ρ = 6860 kg/m 3 , ν = 0.35 and β = 0.54 into Eq.( 2) and Eq.(S1), we obtain f 2 0 versus U as shown in Fig. 2e.For case I, f 2 0 decreases as U increases; while for cases cases II and III, z increases as buckling happens (see Fig. S2b), leading to an increase of f 2 0 .The estimation in Fig. 2e can thus account for all measured results of f 0 versus T for devices D1 to D3.In the following, we describe how to extract the slope of thermal-changed displacement U versus temperature dU dT for cases I to III, which is related to the TEC α m of the membrane through 27 Case I. Pre-buckling regime For case I, the suspended membrane is nearly flat while the change of deflection z with increasing temperature T can be negligible.Therefore, assume dz dT = 0, the derivative of Eq. ( 2) can be simplified as (see details in SI section 2): where c t = 13.34Eπ 2 d 3 ρ(1−ν) .Therefore, in the pre-buckling regime, we can directly extract dU dT from the measured d f 2 0 dT using Eq.(4) (see the flow chart in Fig. 2f).Besides device D1, we also show that devices D4 and D5 are in case I according to their measured f 0 versus T (see SI section 4).
Case II.Transition from pre-to post-buckling For case II, Eq. ( 4) is not applicable anymore since z(T ) varies with temperature.We thus calculate the derivative of Eq. (2) (see SI section 2) and obtain: As depicted in the flow chart of Fig. 2f, we first extract z f ree = 20.6 nm for device D2 from the measured f 0 at the turning point using Eq. ( 2) and Eq.(S1), as well as z versus T (see Fig. S2).The obtained z f ree and z(T ) are then substituted into Eq.( 5) to extract dU dT .The result of U versus T for device D2 shows the expected transition of displacement from tensile (U > 0) to compressive (U < 0), as plotted in Fig. S2.  2) and Eq.(S1) under different z f ree .f 2 0 first decreases and then increases again as U varies from tensile to compressive, which is comparable to the measurement result for device D2. f The proposed procedure to determine thermal properties of 2D material membranes.

Case III. Post-buckling regime
For case III, Eq. ( 5) can be simplified as z 3 ≫ z f ree h 2 , which results in: Thus the result of dU dT for device D3 can be directly extracted from the measured f 0 versus T .The calculated curves in Fig. 2e also verify the linear relations given by Eq. ( 4) and Eq. ( 6).

B. Extracting in-plane thermal conductivity of 2D materials
The flow chart depicted in Fig. 2f also shows how optomechanical measurements as a function of temperature enable a precise pathway for studying the thermal properties of 2D resonators.We first extract the TEC α m of the membrane from the the results of dU dT for cases I to III, which are obtained from Eqs. (4) to (6), respectively, as discussed in the previous section.Then, we quantify the specific heat c v of the membrane from its thermodynamic relation with α m , which will be discussed in this section in more detail.Finally, from the solution of the 2D heat equation, we determine the in-plane thermal conductivity k of the membrane from the measured τ and the obtained c v .In the following, we go step by step through this procedure for device D1.
Let us start with extracting the TEC α m of the membrane.Since the in-plane displacement U originates from the boundary thermal expansion of the membrane, we can extract the TEC α m (T ) of the membrane from the obtained dU dT using Eq. ( 3), where the values of α Si (T ) are taken from literature 36 .The obtained α m versus T for device D1 is shown in Fig. 3a (left).
In the second step, since the specific heat at constant volume is approximately equal to that at constant pressure for solid, we can directly extract the specific heat, C v , of the membrane from the TEC α m using the thermodynamic relation 22 : where K = E 3(1−2ν) is the bulk modulus, V M = M/ρ is the molar volume, M is the atomic mass, and γ is the Grüneisen parameter of the membrane taken from literature.These parameters are listed in Table I for the used materials.Using the obtained α m , we extract C v versus T for device D1, as plotted in Fig. 3a (right).
In the last step, we focus on the heat transport in 2D membranes.As shown in Fig. 3b, we experimentally observe that τ is between 434.6 and 444.0 ns in the probed T range for device D1.Considering the heat transport in a circular membrane, we solve the heat equation in the membrane with an appropriate initial temperature distribution and well-defined boundary conditions (see SI section 3), and obtain the thermal time constant based on the thermal properties of the membrane: where τ rr and τ zz are the in-plane and out-of-plane thermal time constants of the membrane (see Fig. 3b, insert), respectively, c v = C v /M, µ 2 = 5 is the in-plane diffusive constant (see SI section 3), and k is the thermal conductivity of the membrane.Due to the low h/R ratio for 2D materials, we find that τ zz ≪ τ rr and thus the extracted τ from our measurement is equal to τ rr .By substituting the obtained C v and the measured τ into Eq.( 8), we extract k = 8.6 ± 0.3 Wm −1 K −1 for device D1, as plotted in Fig. 3c.
The obtained in-plane thermal conductivity k for all devices D1−D5 are listed in Table I, of which the raw data can be found in Fig. S6.For both 2H-TaS 2 (device D1) and FePS 3 (device D2), since relevant studies on their thermal properties are quite limited, we directly compare the obtained k with the values from the literature 18,39 and observe good agreements (see Fig. 4).For Poly Si (device D3), MoS 2 (device D4) and WSe 2 (device D5), we observe that k depends on the membrane thickness h.We attribute this to a smaller mean free path (MFP) of phonons in thin membranes compared to their bulk counterparts 8 .To account for this effect, we use the Fuchs-Sondheimer model 20,40 that evaluates thermal conductivity of 2D materials as a function of thickness: x dx, (9)  where k bulk and Λ bulk are the thermal conductivity and MFP of bulk, and x is a integration variable.The bulk thermal conductivities k bulk for Poly Si, MoS 2 and WSe 2 are 13.8 Wm −1 K −1 , 98.5 Wm −1 K −1 , and 35.3 Wm −1 K −1 , respectively 19,21,41 .We find that the given k versus h, including our results and literature values [18][19][20][21]39,[42][43][44] , is well described by Eq. ( 9) as indicated by the fitted solid lines in Fig. 4 using Λ bulk as fit parameter, obtaining 75 nm, 19 nm, and 19 nm for Poly Si, MoS 2 and WSe 2 , respectively. These fited values of Λ bulk are also in good agreement with previously reported phonon MFPs 19,45,46 , supporting the validity of employing Eq. ( 9) to predict the thickness dependent thermal conductivity of 2D materials.

IV. DISCUSSION
Compared to other methods for determining the thermal conductivity of 2D materials, the optomechanical approach has several advantages, as summarized in Table II.For the Raman microscopy method, since relatively large temperature changes are needed to resolve the resulting shift in Raman mode frequency, a very wide temperature range has to be measured to get an accurate slope χ T of the Raman peak shift with temperature.For example, χ T for MoS 2 is -0.013 cm −1 /K.
Considering a limited resolution 0.25 cm −1 for a Raman microscope, a temperature increase of at least 20 K is required to obtain meaningful results 47 .In our measurements, we require only a narrow T range to study the thermal transport (see Fig. S6).For the micro-bridge method, either thick crystals or stiff 2D materials like graphene are required to survive the complicated fabrication procedures including lithography and etching 1 .In contrast, for the presented contactless optomechanical method, one only needs to suspend membranes over cavities in a Si substrate, which is applicable for most 2D materials and can be done for any thickness.
Although we estimate the average MFP for bulk in Fig. 4, we note that the phonon MFP in 2D materials is highly related to the phonon dispersion relation, surface strain, crystal grain size, and temperature.These factors can be further studied using the presented optomechanical approach, which would help us to better understand the phonon scattering mechanisms in 2D materials.Moreover, our work suggests a new way to further investigate acoustic phonon transport in recently emerged 2D materials, such as phosphorene and MXenes with distinct thermal anisotropy 48,49 , as well as the magic-angle multilayer superconductor family 50 .It is also of interest to probe the dynamics of phonons across the interface in vdW heterostructures, so as to realize a coherent control of thermal transport across 2D interfaces 51,52 .

V. CONCLUSIONS
We demonstrated an optomechanical approach for probing the thermal transport in 2D nanomechanical resonators made of few-layer 2H-TaS 2 , FePS 3 , Poly Si, MoS 2 , and WSe 2 .We measured the resonance frequency and thermal time constant of the devices as a function of temperature, which are used to extract their thermal expansion coefficient, specific heat, as well as in-plane thermal conductivity.The obtained values of all these parameters (see Table I) are in good agreement with values reported in the literature.Compared to other methods for characterizing the thermal properties of 2D materials, the presented contactless optomechanical approach requires a smaller temperature range, allows for easy sample fabrication, and is applicable to any 2D material.This work not only advances the fundamental understanding of phonon transport in 2D materials, but potentially also enables studies into the use of strain engineering and heterostructures for controlling heat flow in 2D materials.A Si wafer with 285 nm dry SiO 2 is spin coated with positive e-beam resist and exposed by electron-beam lithography.Afterwards, the SiO 2 layer without protection is completely etched using CHF 3 and Ar plasma in a reactive ion etcher.The edges of cavities are examined to be well-defined by scanning electron microscopy (SEM) and AFM.After resist removal, 2D nanoflakes are exfoliated by Scotch tape, and then separately transferred onto the substrate at room temperature through a deterministic dry stamping technique.More details on the substrate fabrication and Scotach tape transfer method can be found in our previous work 53 .Using tapping mode atomic force microscopy (AFM), we measure the height difference between the membrane and the Si/SiO 2 substrate.As Fig. S1 shows, we find a membrane thickness h of 24.0 nm for device D3.
The fabricated devices is then fixed on a sample holder inside the vacuum chamber.A PID heater and a temperature sensor are connected with the sample holder, which allows to precisely monitor and control the temperature sweeping.A piezo-electric actuator below the sample holder is used to optimize the in-plane XY position of the sample to maintain blue and red lasers irradiating on the center of the circular sample.
Section 2: Mechanical model of 2D material membranes under cases I to III In the main text, we observe different types of dynamic response in the measured devices D1−D3, attributed to the thermallyinduced buckling in 2D nanomechanical resonators.Eq. ( 2) gives the expression of the resonance frequency f 0 (T ) in buckled resonators.Accodrding to our previous work 35 , Eq. ( 2) can be solved by the relation between central deflection z of the membrane and compressive displacement U from the boundary: where U is the thermally induced in-plane displacement of the plate, ρ is the mass density, z is the central deflection of the plate, ν is the Poisson ratio, z f ree is the central deflection of free plate without loading (when U = 0), and β is a fitting factor determined by ν = 0.35ν + 0.42.Therefore, using Eq.(S1) and Eq. ( 2), we can extract z f ree = 20.6 nm and z(T ) from the measured f 0 (T ) for FePS 3 device D2.As plotted in Fig. S2, we see as the boundary displacement U varies from tensile to compressive with increasing T , the central deflection z of the membrane gradually goes up.More details of the Galerkin buckling model can be found in our previous work 35 .
For pre-buckling regime (case I), assume the deflection z nearly keeps constant, the only time-dependent parameter in Eq. ( 2) is U, which allows us to obtain the derivative d f 2 0 dT = c t dU dT for case I.For transition regime from pre-to post-buckling (case II), we first calculate the derivative of Eq. (S1): Therefore, by substituting Eq. (S2) into the T -derivative of Eq. ( 2), we obtain Eq. ( 5).For post-buckling regime (case III), the thermally-induced buckling results in a large central deflection of the membrane, thus we assume z f ree h 2 ≪ z 3 and simplify Eq. ( 5) as: In this section, we explain how we derive thermal time constant τ with respect to the thermal properties of 2D membrane, as given in Eq. ( 8) in the main text.
Consider the situation that a modulated-laser irradiates at the center of the suspended 2D membrane, the Fourier heat conduction equation in cylindrical coordinate for this problem can be written as: where u(r, z,t) is time t-dependent temperature distribution in the membrane along with the radial r and perpendicular z directions, κ = k c p ρ is thermal diffusivity, dQ dt is the absorbed heat energy per unit time per unit volume at the center of the membrane, c p and k denote the specific heat and thermal conductivity of the membrane, respectively.

Transient state
We firstly discuss the transient state of heat transport in the membrane, corresponding to a laser pulse irradiates.As a result, Eq. (S4) is rewritten as: with the boundary conditions: and the initial condition:

Dependence on laser powers
Furthermore, we discuss the effect of laser powers on our optomechanical measurement, as depicted in Fig. S4.According to COMSOL simulation, we verify that the heating only plays a role on the amplitude of thermal signal, instead of its location which is related to τ (Fig. S4a).Next, we test MoS 2 device D4, and plot the extract τ as the function of red and blue laser powers, respectively (see Figs. S4b and S4c).As expected, the obtained τ nearly keep constant and are independent to both P r and P b .Therefore, we verify that the proposed optomechanical methodology does not need a laser calibration for determining the in-plane thermal conductivity of 2D materials.In Figs.S5a and S5b, we see that the measured f 0 decreases as T increases for devices D4 (MoS 2 ) and D5 (WSe 2 ).Therefore, both of them are in pre-buckling regime (case I), same as device D1 (2H-TaS 2 ).In addition, one could argue that the observed increase of f 0 with T for device D3 (Fig. 2c in the main text) is attributed to α m < α Si , instead of a post-buckling performance.As shown in Fig. S5c, we obtain the first decrease and then increase of f 0 with T increasing (case II) for another measured Poly Si device, which verify that α m > α Si for Poly Si and the observed increase of f 0 in device D3 corresponds to the mechanical response in post-buckling regime.
Figure .S6 shows the obtained specific heat c v , the measured τ and the obtained in-plane thermal conductivity as the function of T for the fabricated devices D2−D5, respectively.Corresponding results for 2H-TaS 2 device D1 has been given in Fig. 3  the main text.

FIG. 1 .
FIG. 1. Sample characterization and experimental setup.a Top, schematic diagram of 2D nanomechanical resonators, composed of 2D flake suspended on the etched SiO 2 /Si cavities; bottom, optical image of device D1 (2H-TaS 2 ).Scale bar is 5 µm.b AFM indentation results for device D1 (points), from which the Young's modulus E of the membrane is extracted by fitting the measured force F to the cantilever deflection δ (drawn line).c Laser interferometry setup used for the optomechanical measurement.The sample is mounted in a vacuum chamber (VC) with a pressure below 10 −5 mbar.The reflected red laser is detected by the photodiode (PD) and input to the vector network analyzer (VNA).PBS, polarized beam splitter; DM, Dirac mirror.d Resonant peak of device D1 measured at MHz regime (points), which is fitted with a harmonic model (drawn line) to extract the fundamental resonance frequency f 0 of device D1. e Thermal signal measured at kHz regime, including imaginary (red points) and real (blue points) parts.The imaginary part is fitted with Eq. (1) (drawn lines) to obtain the thermal time constant τ of device D1.

FIG. 2 .
FIG. 2. Optomechanical methodology for obtaining thermal properties.a-c Measured resonant peak measured as the function of temperature T for devices D1 to D3, corresponding to cases I to III, respectively.d Schematic diagram of the buckled device, where the central deflection z of the membrane increases as the boundary displacement U loaded.e Squared resonance frequency f 2 0 as the function of U in the membrane estimated by Eq. (2) and Eq.(S1) under different z f ree .f 2 0 first decreases and then increases again as U varies from tensile to compressive, which is comparable to the measurement result for device D2. f The proposed procedure to determine thermal properties of 2D material membranes.

FIG. 3 .
FIG. 3. Quantifying thermal characteristics of device D1. a Thermal expansion coefficient α m and specific heat C v of 2H-TaS 2 membrane versus temperature T .b Measured thermal time constant τ versus T .Insert, schematic diagram of heat transport in suspended 2D membrane.c In-plane thermal conductivity k versus T extract from Eq. (8) using the measured τ and the obtained C v .
FIG. S2.Analysis of thermally-induced buckling for device D2.(a) Points, measured fundamental resonance frequency f 0 versus temperature T ; drawn line, theoretical curve estimated by Eq. (2) and Eq.(S1) using the TEC of FePS 3 α m = 1.1 × 10 −5 K −1 .(b) Central deflection z of the membrane versus T .(c) Boundary displacement of the membrane U versus T .
FIG. S3.(a) 3D schematic diagram of 2D nanomechanical resonators in COMSOL simulation software.Insert, side view, where the boundary condition changes to the bottom of SiO 2 /Si substrate.(b) Average temperature of the membrane as the function of heating rate f .Points, simulation results; lines, fit with Eq. (1) in the main text to the simulation.

Section 4 :
Raw data of optomechanical measurements for devices D2 to D5 FIG. S5.Measured resonance frequency f 0 as the function temperature T for (a) MoS 2 device D4, (b) WSe 2 device D5, as well as (c) another Poly Si device with buckling transition.a b c

TABLE I .
Characteristics of devices D1−D5, including radius R, thickness h, mass density ρ, Young's modulus E, atomic mass M, Poisson ratio ν, Grüneisen parameter γ, as well as the obtained average TEC α m , specific heat C v and in-plane thermal conductivity k.The values of ρ, M, ν and γ are taken from literature

TABLE II .
Comparison of different thermal conductivity measurement methods, where the required temperature range is quantified by the studies of MoS 2 .