**1. Introduction**

The thermal conductivity of 2D materials is of grea<sup>t</sup> significance for both basic research [1–10] and practical application [11–13]. In basic research, Fourier's law has been successful in studying heat conduction in macroscopic systems. However, when down to micro or nanoscale, due to the existence of size effect [1], thermal rectification [2], and ballistic transport [3], this law is no longer usable. In addition, with the advent of new materials such as newly discovered borophene [4,5], MXene [6,7], and various heterostructures [8–10], it is also crucial to determine the thermal conductivity of these new materials. From a practical perspective, 2D materials are widely used in optoelectronic devices [11], biological monitoring [12], and energy storage [13] due to their excellent optical, electrical, and mechanical properties. It is necessary to explore the thermal conductivities of these 2D materials to optimize heat dissipation in optoelectronic devices.

Various theoretical calculation methods such as molecular dynamics simulation [14–17], phonon Boltzmann transport equation [18–23], and atomistic Green's functions [24–26] have been developed to study the underlying physical mechanism of heat transfer in 2D materials. Yet, due to the ignorance of surface defects, the accuracy of these methods is limited. Experimental methods, such as the suspended micro-bridge method [27–32], 3ω method [33–36], timedomain thermoreflectance method [37–46], and Raman method [47–53], have been developed to study the thermal conductivity of 2D materials. Bao et al. [54] introduced heat transfer research methods in micro-nano structures from the perspective of theoretical calculation. Experimental-based thermal characterization techniques for low-dimensional materials were also reviewed [55,56]. Considering all these different techniques, however, there is still a lack of comprehensive review on the thermal conductivity measurement methods of 2D materials. In

**Citation:** Dai, H.; Wang, R. Methods for Measuring Thermal Conductivity of Two-Dimensional Materials: A Review. *Nanomaterials* **2022**, *12*, 589. https://doi.org/10.3390/nano12040589

Academic Editor: Gyaneshwar P. Srivastava

Received: 31 December 2021 Accepted: 28 January 2022 Published: 9 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

153

this paper, both theoretical and experimental methods for studying the thermal conductivity of 2D materials are reviewed. In addition, the factors affecting the thermal conductivity of 2D materials are also discussed.

#### **2. Theoretical Methods**

The theoretical calculation is an effective way to deeply understand the potential mechanism of phonon transport in 2D materials. Currently, the molecular dynamics (MD) simulation, phonon Boltzmann transport equation (PBTE), and atomistic Green's functions (AGF) were the 3 mainstream theoretical methods.

#### *2.1. MD Simulation*

In MD simulation, the motion of each particle in the dynamic process was described based on Newton's second law, and the position, velocity, and force of each atom were calculated at each step. The interatomic forces were derived from the potential function, and the commonly used empirical potential functions were Lennard–Jones (LJ) potential for interlayer van der Waals (vdW) interaction, Stillinger–Weber (SW) [15] potential for atomic interaction, and REBO potential [16] for the covalent bonding of the carbon atoms in diamond and graphite. Two common methods were used to calculate the thermal conductivity of 2D materials: the equilibrium MD (EMD) method based on the Green– Kubo formalism and the nonequilibrium MD (NEMD) method based on Fourier's law.

In EMD method, the thermal conductivity is expressed as the integration of the heat current autocorrelation function (HCACF) with respect to a given correlation time *t*,

$$\kappa\_{\text{\tiny\mu\text{V}}}(t) = \frac{1}{\kappa\_B T^2 V} \int\_0^t \langle J\_{\mu}(0) J\_{\text{\tiny\text{\tiny\mu\text{V}}}}(t') \rangle dt' \tag{1}$$

where κ*B* is Boltzmann's constant, *T* is the absolute temperature of the system, *V* is the volume and *J*μ is the μth component of the full heat current vector *J*. The heat current at a given time depends on the positions and velocities of the particles in the system. The key point of EMD based thermal conductivity calculation is to calculate the time integral with the upper limit of infinity in Equation (1) and ensure its convergence. In addition, the size effect of thermal conductivity is also difficult to be studied in EMD, which can be solved in NEMD. The NEMD technique can be employed to characterize the in-plane thermal conductivity of a sample with finite length L by driving the system out of equilibrium. When steady state is achieved after sufficient time, the heat current (flux) Q and temperature gradient ∇ *T* are obtained to calculate thermal conductivity κ(*L*) according to Fourier's law:

$$\kappa(L) = -\frac{\mathcal{Q}}{\nabla T} \tag{2}$$

For the MD method, one advantage is that the simulation is based on a real physical model in space, which makes it convenient to study the effects of physical parameters, such as strain, defect, doping, etc. However, the accuracy of MD simulation is highly dependent on the potential empirical function used, which is usually developed by fitting the existing material properties. Moreover, as the MD method is modeled in real space, the calculation range is limited due to the simulation time, which makes the calculation of thermal conductivity not very accurate. In addition, the phonon scattering rate in MD is related to the Maxwell Boltzmann distribution while ignoring the quantum effect below the Debye temperature. Thus, erroneous results were obtained for the calculated thermal conductivity below the Debye temperature [17].

#### *2.2. PBTE Method*

Recently, PBTE, combined with first principles, was used much more frequently to explore the thermal conductivity of non-metallic materials, such as 2D selinene [18], phosphorene [19], borophane [20], and transition metal dichalcogenide (TMDC) MX2 [21]. Under the effect of temperature gradient (∇*T*), the phonon distribution *fλ* deviates from the Bose–Einstein distribution in equilibrium *f* 0*λ*, which can be obtained by solving PBTE:

$$\left.\frac{\partial f\_{\lambda}}{\partial t}\right|\_{\text{diff}} + \left.\frac{\partial f\_{\lambda}}{\partial t}\right|\_{\text{scatt}} = 0\tag{3}$$

where the diffusion term ( *∂ fλ ∂t* diff) is caused by the temperature gradient ∇*T* and given by:

$$\frac{\partial f\_{\lambda}}{\partial t}\Big|\_{\text{diff}} = -v\_{\lambda} \nabla T \frac{\partial f\_{\lambda}^{0}}{\partial T} \tag{4}$$

where *vλ* is the group velocity of phonon mode *λ*. The scattering term ( *∂ fλ ∂t* scatt) in Equation (5) is determined by the scattering process in the system. Under the relaxation time approximation (RTA), the scattering term can be written as:

$$\left.\frac{\partial f\_{\lambda}}{\partial t}\right|\_{\text{scatt}} = \left.\frac{f\_{\lambda} - f\_{\lambda}^{0}}{\tau\_{\lambda}}\right|\tag{5}$$

where τ*λ* is the relaxation time. Considering anharmonic phonon–phonon interactions τ*λ* can be obtained by perturbation theory. The heat flow *J*α in α direction can be written as:

$$J^{\alpha} = \sum\_{\lambda} \int \hbar \omega\_{\lambda} f\_{\lambda} v\_{\lambda}^{\alpha} \frac{dk}{2\pi^{3}} \tag{6}$$

where *k* denotes the phonon wave vector, ω*λ* is the frequency of phonon mode. According to Fourier's law *J*α = −∑<sup>β</sup> <sup>κ</sup>αβ(∇*T*)<sup>β</sup>, the lattice thermal conductivity tensor καβ under the RTA can be written as:

$$\kappa^{\alpha\beta} = \frac{1}{k\_B T^2 N V} \sum\_{\lambda} (\hbar \omega\_{\lambda})^2 f\_{\lambda}^0 (1 + f\_{\lambda}^0) v\_{\lambda}^{\alpha} v\_{\lambda}^{\beta} \,\mathrm{\tau}\_{\lambda} \tag{7}$$

where *N* is the total number of phonon wave vectors included in the summation, *V* is the volume of the unit cell, <sup>ν</sup><sup>α</sup>*λ* and νβ*λ* are the group velocity of phonon mode *λ* with Cartesian coordinates indexed by α and β, respectively. In the actual simulation, the parameters <sup>ν</sup><sup>α</sup>*λ* and νβ*λ* were obtained from the interatomic force constants (IFCs), which can be extracted from DFT packages such as VASP. Some open-source software packages such as ShengBTE [22] were available to predict the lattice thermal conductivity of solid materials with the input files of these IFCs.

In PBTE, the calculation accuracy depends on the accuracy of the scattering mechanism in the 2D material. Anharmonicity causes inelastic scattering of phonons. Meanwhile, many factors such as isotopes, holes, and interfaces may disturb the lattice vibration. At present, PBTE lacks the description of some scattering mechanisms, such as holes [23].

#### *2.3. AGF Method*

The AGF method, which is based on a dynamical equation and the quantum mechanical phonon energy distribution, is an effective tool to simulate ballistic transport in nanoscale devices. As shown in Figure 1, the quantum thermal transport system can be divided into 3 parts: central scattering region (abbreviated as C), left and right lead (abbreviated as L, R). Under the harmonic approximation, the phonon waves in the system can be described as:

$$(\omega^2 I - H)\Phi(\omega) = \mathbf{0} \tag{8}$$

where *ω* is the angular frequency of lattice vibration, *I* is the identity matrix, *H* is the harmonic matrix, and **<sup>Φ</sup>**(ω) is the eigenvector of *H*. The response of the system under small disturbance can be obtained by Green's function:

$$(\omega^2 I - H)\mathbf{G} = \mathbf{0} \tag{9}$$

**Figure 1.** Schematic diagram of heat transport model for low dimensional system. (The L, C and R in the figure represents three parts of the system: central scattering region, left and right lead, and Q indicates the direction of heat flow).

The atomic interactions in each region are described by constructing a harmonic matrix for AGF calculation [24]. The phonon transmission function <sup>Ξ</sup>(ω) is calculated by:

$$\Xi(\omega) \;= \; \text{Tr} [\Gamma\_L G\_\mathbb{C}^\prime \Gamma\_R G\_\mathbb{C}^\prime \; ^\*] \tag{10}$$

where *GrC* and *GrC*∗ are the Green's function of the central region and its complex conjugate [25], Γ*L* and Γ*R* are phonon escape rates from left contact and right contact, and *Tr* represents the trace of the matrix.

According to the Landauer formula and the phonon transmission function, the thermal conductivity κ of the system can be calculated as:

$$\kappa = \frac{\hbar}{2\pi} \int\_0^\infty \frac{\partial f(\omega, T)}{\partial T} \Xi(\omega) \omega d\omega \tag{11}$$

where *f*(<sup>ω</sup>, *T*) is the Bose–Einstein distribution and - is the Planck's constant. AGF studies the thermal conductivity based on the harmonic approximation condition and does not consider the anharmonic interaction, that is, phonon–phonon scattering. Therefore, AGF, which studies the structure dominated by elastic scattering, is mainly used in nanostructures dominated by harmonic scatterings, such as defects, interfaces, lattice mismatch [26].

#### **3. Experimental Methods**

The theoretical simulation methods, including MD, PBTE, and AGF, have become effective tools for calculating the thermal conductivity of 2D materials. However, it is challenging to ensure the accuracy when considering the impurities, defects, and rough surface of real samples. That is, it is of grea<sup>t</sup> significance to developing experimental methods to improve the measuring accuracy. The measuring accuracy can also be further improved by combining the experimental methods with the theoretical calculation. At present, experimental measurement methods mainly include electro-thermal and optothermal methods.

#### *3.1. Electro-Thermal Techniques*

Electro-thermal techniques, which include the suspended micro-bridge method, 3ω method, electron beam self-heating, T-bridge, four-probe transport measurements techniques, characterize the thermal conductivities of materials based on the temperature dependence of thermal resistance. The suspended micro-bridge method and 3ω method are two typical techniques and are introduced in detail.

3.1.1. Suspended Micro-Bridge Method

The suspended micro-bridge method was first used by Majumdar et al. [27,28] in 2001 to measure the thermal conductivity of a single multi-walled nanotube. Since then, the suspended micro-bridge method has been applied to measure the thermal conductivities of graphene [29], hexagonal boron nitride [30], MoS2 [31] and other 2D materials. As shown in Figure 2, the suspended device is composed of 2 adjacent silicon nitride (SiNx) membranes suspended by 5 SiNx beams. The platinum resistance thermometer coil designed on each membrane is connected to the substrate through a platinum (Pt) leads on the long SiNx beam. A mixed current of DC (microampere level) and AC (nano ampere level) is introduced to the heating membrane.

**Figure 2.** Schematic diagram of microbridge measurement. Reprinted with permission from Ref. [32]. Copyright 2017, John Wiley and Sons.

The DC current is used to generate Joule heat (*Qtot*) on one side, and the Pt resistance is measured by AC to characterize the temperature change (Δ*Th*, Δ*Ts*) of heating and sensing membrane caused by *Qtot*. Heat is transferred between the heating membrane and the sensing membrane only through the sample. Since the *Qtot* is transferred only from the heating membrane to the substrate with an environment temperature *Ta* and sample, we can express the heat flux distribution on the whole device and sample as follows:

$$Q\_{tot} = Q\_1 + Q\_2 \tag{12}$$

$$Q\_1 = G\_b \times \Delta T\_h \tag{13}$$

$$Q\_2 = G\_5(\Delta T\_h - \Delta T\_s) \ = G\_b \times \Delta T\_s \tag{14}$$

where *Qtot* is the total heat on the heating membrane, *Q*1 and *Q*2 are the heat transferred from the heating membrane to the substrate and the sample, *G*s and *Gb* are the conductance of the sample and SiNx beams, respectively. The thermal conductivity (κ) of the sample can be written as:

κ

$$\mathbf{I} = \mathbf{G}\_s \frac{\mathbf{L}}{\mathbf{S}} \tag{15}$$

where *L* and *S* are the length and sectional area of the sample, respectively. In practice, there is thermal resistance ( *R*c) at the interface between the sample and the heating/sensing membrane. The measured total thermal resistance is *R* = *RG* + 2*RC*, where *RG* is the actual thermal resistance of the sample. Some methods have been designed to reduce the effect of *RC*. One is to calculate the temperature rise between the sample and membrane through numerical simulation [30]. In addition, it can be considered to add high thermal conductivity materials to the membrane to reduce *RC* and improve the uniformity of membrane temperature [29].

#### 3.1.2. 3 ω Method

The 3 ω method is based on the frequency-domain feedback characteristic that the temperature of the heating resistor varies with the frequency of the applied AC electrical current. As shown in Figure 3a, a metal electrode such as Pt with a certain shape and thickness (the yellow part) was prepared on the surface of the thin film sample (the blue part) by photolithography and thermal evaporation, which was used both as a heater and a thermometer. Thin-film samples are usually deposited on the substrate (the bottom gray part) by chemical vapor deposition (CVD) and high-temperature oxidation. When an AC power supply with a frequency of 1 ω is connected to the metal electrode, the internal resistance of the metal electrode changes approximately at a frequency of 2 ω due to the linear relationship with the temperature change. Finally, the voltage signal with 3 ω frequency variation can be extracted by the lock-in amplifier (shown in Figure 3b).

**Figure 3.** Schematic diagram of (**a**) 3 ω method. (**b**) experimental circuit. Reprinted with permission from Ref. [33]. Copyright 2008, AIP Publishing.

In this technique, 2 structures were prepared: substrate and film-substrate structure. Metal electrodes were deposited on the 2 structures to measure the corresponding temperature changes ( Δ *Ts*, Δ *Ts* + *f*). Then, temperature change caused by the film can be written as Δ *Tf* = Δ *Ts* + *f* − Δ *Ts*. The thermal conductivity (<sup>κ</sup>*f*) of a thin film is determined using Equation (16)

$$\kappa\_f = \frac{Pt}{\Delta T\_f \cdot S} \tag{16}$$

where *P* and *t* are heating power and film thickness, respectively.

For the 3 ω method, thermal contact resistance measurements between graphene and SiO2 based on a differential 3 ω technique were made [34]. However, as the fabrication of a

metal electrode with high quality and the signal extraction of the phase-locked amplifier were required, it was difficult to measure the thermal conductivity of 2D material with atomic level thickness. Zhang et al. [35] reported the thermal conductivity measurement of 100 nm thickness silicon nitride (SiN) and 64 nm thickness amorphous boron nitride (BN) based on the 3ω method. As the thermal conductivity was obtained by the frequencydependent temperature oscillation, the 3ω method was free of the effect from contact thermal resistance between sample and substrate. Then, due to the small surface area of metal electrode, the effect of heat radiation was also limited [36]. The 3ω method also has some drawbacks that further limit its application for supported samples. The thermal conductivity of the substrate should be much higher than that of the film deposited on it to ensure a high sensitivity. A lower surface roughness of the sample is needed to prevent damage to the thin metal wires.

#### *3.2. Opto-Thermal Techniques*

Compared with electrothermal method, opto-thermal techniques, which can realize non-contact measurement with simple sample preparation, have been widely used in thermal conductivity characterization of 2D materials. Two representative methods, timedomain thermoreflectance (TDTR) method, and Raman-based methods, are introduced in detail in this section.

#### 3.2.1. Time-Domain Thermoreflectance (TDTR) Method

The TDTR method is based on the change of surface reflectance caused by temperature change. Figure 4a shows a typical setup. The emitted laser is divided into pump light and probe light through a polarizing beam splitter (PBS). The pump light modulated by the electro-optic modulator is used to heat the sample surface. The sample surface is usually covered with a metal film in order to ensure that the pump laser is absorbed at the surface. The detection beam is delayed relative to the pump light by the mechanical delay stage and received by the photodiode detector. The converted electrical signal is extracted by the lock-in amplifier with two outputs: in-phase (*V*in) signal and out-of-phase (*V*out) signal, which represent the phase of the reflected beam related to the temperature response and can be written as *R* = −*V*in/*V*out. By continuously changing the delay time, the curve of R versus time shown in Figure 4b can be obtained. Combined with the heat transfer model established by Cahil et al. [37] in 2004, the thermal conductivity can be extracted. Schmidt et al. [38] further applied the model in anisotropic thermal conduction of highly ordered pyrolytic graphite (HOPG).

**Figure 4.** (**a**) Schematic of a typical TDTR setup. (**b**) The ratio between in-phase and out-of-phase signals, −*V*in/*V*out as a function of delay time is compared with the thermal modeling to extract the thermal conductivity. Reprinted with permission from Ref. [39]. Copyright 2021, AIP Publishing.

In TDTR, due to the deposition of metal films on the sample surface, the intrinsic thermal conductivity of the sample cannot be measured accurately. The main limitation for TDTR is the requirement for a highly smooth surface to minimize diffuse reflection and the complex experimental device. In the later development, many improvements were made based on TDTR, such as FDTR [40] TDTR based on time-resolved magneto-optical Kerr effect (TR-MOKE) [41]. In FDTR, the relationship between thermal reflection signal and modulation frequency rather than the delay time was established, in which continuous laser can be used thus it is simpler and cheaper. TDTR can also be combined with TR-MOKE to probe the sample's surface temperature, which depends on the temperature dependence of the polarization rather than the intensity of the reflected beam. This temperature measurement method allows us to use thinner ferromagnetic metal film as a transducer, reduces lateral heat flow in the metal film, and improves the accuracy of measurement results. Currently, for thermal conductivity measurement, the TDTR method is mostly used in thin films, such as transition metal dichalcogenides MX2 (M = Mo, W and X = S, Se) [42], h-BN [43], BP [44], which requires a relatively large thickness (>100 nm). For 2D materials, this technique can realize the characterization of the interfacial heat transfer [45]. Additionally, FDTR, which is an improved TDTR method, is used in measuring the thermal conductivity of 2D materials. Rahman et al. [46] implemented frequency domain magnetooptical Kerr effect (FD-MOKE) to measure the thermal conductivity of various 2D materials, such as graphene, monolayer MoS2, and four-layer h-BN.

#### 3.2.2. Optothermal Raman Methods

Compared with the complex measurement device of TDTR, Raman-based methods are simpler and have been widely used in the thermal conductivity measurement of 2D materials. By constructing different heat transfer states in the time and space domain, various Raman-based measurement methods were developed. Among them, the optothermal Raman method based on steady-state heating is the most commonly used.

In this method, the sample can be heated optically or electrically. Taking optical heating as an example, Figure 5a shows the MoS2 sample is suspended on a Si2N4 substrate and heated by a focused laser light. The heat can only diffuse around the sample and eventually to the substrate. As shown in Figure 5c, the Raman shift (Δω) of MoS2 is linearly related to the local temperature change (Δ*T*) of the sample upon laser heating, and can be written as: Δω = χ*T*Δ*T*, where χ*T* is the first-order temperature coefficient. Varying laser power will also produce different thermal effects, which means that there is a similar linear relationship between Raman shift and laser power (Figure 5d). For the sample suspended on a hole with radius *R*, the temperature at *r* from the center of the hole can be calculated from the heat conduction equation as follows:

$$\kappa\_{\rm sus} \frac{1}{r} \frac{d}{dr} \left[ r \frac{dT(r)}{dr} \right] + q(r) = 0, \text{ for } r \le R \tag{17}$$

$$\kappa\_{\rm sup} \frac{1}{r} \frac{d}{dr} \left[ r \frac{dT(r)}{dr} \right] + \frac{\mathcal{G}}{t} [T(r) - T(a)] = 0, \text{ for } r \ge R \tag{18}$$

where κsus, <sup>κ</sup>sup, *q*(*r*),*g*, *t* and *<sup>T</sup>*(*a*) are the thermal conductivity of suspended and supported structure, volume optical heating, interface thermal conductivity, thickness of sample and environmental temperature, respectively. Figure 5b shows the calculated temperature distribution of the sample. The weighted average temperature rise in the laser spot can be written as:

$$T\_{calculated} \approx \frac{\int\_0^R T(r)q(r)r dr}{\int\_0^R q(r)r dr} \tag{19}$$

**Figure 5.** Illustration of optothermal Raman methods. (**a**) Schematic of the thermal conductivity measurement showing suspended MoS2 flakes and excitation laser light. (**b**) Simulation of laser heating temperature rise. (**c**) Raman peak frequency shift as a function of temperature. (**d**) Experimental data for Raman shift as a function of laser power, which determines the local temperature rise in response to the dissipated power. Reprinted with permission from Ref. [47]. Copyright 2014, American Chemical Society.

By matching the calculated temperature rise (*Tcalculated* − *Ta*) with the temperature rise measured by Raman spectroscopy (Δ*Tmeasured*), the thermal conductivity κ can be extracted.

One drawback of optothermal Raman method is the measurement of absolute laser absorption power. Under laser heating, part of the energy is absorbed by the sample, while the rest of the energy is reflected by the sample or transmitted to the substrate. Currently, it is very difficult to determine the laser absorption coefficient accurately. In addition, a temperature calibration process, which is time-consuming and can introduce large errors, is also needed.

#### 3.2.3. Time-Resolved Raman Methods

Time-domain differential Raman (TD-Raman) [48], which uses a square wave modulated laser with variable duty cycle, can be applied to measure the thermal conductivity of 2D materials. As shown in Figure 6a, the modulated laser is used for sample heating and Raman excitation, which consists of a variable excitation period *te* and a fixed thermal relaxation period *tr*. Here, the thermal relaxation time *tr* needs to be long enough for the sample to cool completely before the next pulse period. Figure 6b shows the corresponding temporally accumulative Raman spectra of one laser pulse cycle in 3 cases. It can be seen that longer excitation time *te* leads to higher temperature rise, and the corresponding Raman spectra also change. From cases 1 to 3, the intensity of the Raman peak increases gradually, and the softening phenomenon of Raman peak position is also observed. By analyzing the changes of Raman signals mentioned above, the average temperature rise (Δ*T*) of the sample in the heating zone can be determined by Raman spectroscopy. Moreover, the accumulative Raman emission for one excitation cycle (0~*te*) is written as:

$$E\_{\omega}(\omega, t\_{\varepsilon}) = I\_0 \int\_0^{t\_{\varepsilon}} \left(1 - A\Delta T^\*\right) \exp\left[-\frac{4ln2\cdot(\omega - \omega\_0 + B\Delta T^\*)^2}{\left(\Gamma\_0 + C\Delta T^\*\right)^2}\right] dt\tag{20}$$

where *I*0, ω0, Γ0 are the corresponding Raman properties at the beginning of laser heating, A, B, C are the changing rate of Raman intensity, Raman shift, and linewidth against the normalized temperature Δ*T*<sup>∗</sup>. Combining with the transient heat transfer model, the thermal conductivity of the sample can be determined by fitting the variation of Raman peak with time.

**Figure 6.** Concept of TD−Raman. (**a**) The change of temperature evolution(ΔT), and instant changes of Raman peak intensity (I), peak shift (ω) and linewidth (Γ). (**b**) The corresponding temporally accumulative Raman spectra of one laser pulse cycle in Case 1, 2, and 3. Reprinted with permission from Ref. [48] © The Optical Society. Copyright 2015, The Optical Society.

In TD-Raman, thermal conductivity of 2D materials is measured through the Raman characterization of transient heat transfer. However, in practice, when the heating time is too short, a long time of Raman signal acquisition is needed, which makes it hard for fast thermal transport characterization produces more environmental interference and affects the measurement accuracy. To solve this problem, Wang's lab further developed frequency-resolved (FR) Raman technique Frequency-resolved Raman for transient thermal probing and thermal diffusivity measurement [49]. As shown in Figure 7, an amplitudemodulated square-wave with different frequencies is employed to heat the sample and excite Raman signals. When the sample is irradiated by the high-frequency laser pulse, the temperature of the sample is almost constant in the whole process, which is defined, as "quasi-steady state," and the temperature rise is regarded as <sup>Δ</sup>*Tqs*. On the contrary, when low-frequency laser pulse irradiates the sample because the sample has enough time to rise to a stable state in the excitation time, the temperature of the sample is approximately regarded as a constant in the excitation time, which is defined as "steady state," and the temperature rise is regarded as Δ*Ts*. Here we have Δ*Tqs* = Δ*T*s/2, which shows that the temperature decreases with increasing frequency. The thermal conductivity can be extracted based on the transient heat transfer model and the collected Raman signal. TD-Raman has been applied to the measurement of the anisotropic thermal conductivity of

black phosphorus [50]. Compared with TD-Raman, the Raman signal collection of FR-Raman is more efficient, but the sensitivity is lower because the time between pulses is not enough to completely cool the sample.

**Figure 7.** Concept of FR-Raman. (**a**) Time profiles of laser pulse, temperature evolution (*T*), Raman peak intensity (*I*), peak shift (*ω*) and linewidth (Γ). (**b**) Temperature variation (*TQS*) at quasi-steady state. (**c**) Temperature variation (*TS*) at very low frequency. Reprinted with permission from [49] © The Optical Society. Copyright 2016, The Optical Society.

#### 3.2.4. Energy Transport State Resolved Raman (ET-Raman)

Besides the time-domain modulation, the energy transport states can also be modulated in the spatial domain. Based on this, a technique named ET-Raman was developed to measure the in-plane thermal conductivities of supported or suspended 2D materials. For supported 2D samples, both a CW laser and a picosecond laser are used. As shown in Figure 8, 5 energy transport states were constructed both in time and spatial domains [51]. Three physical processes occur with laser heating. The first is hot carrier generation, diffusion in space, and electron–hole recombination. This process introduces heat transfer and energy redistribution, which is determined by the hot carrier diffusivity (*D*). The subsequent process is the heat conduction by phonons, which receives energy from the hot carriers or electron–hole recombination, which mainly happens in the in-plane and depends on the thermal conductivity (*κ*). The third is the heat conduction from sample to substrate, and this process is dominated by the local thermal resistance (*R*).

By using different laser power (*P*), a parameter named Raman shift power coefficient (RSC) can be obtained and expressed as: χ = *∂*ω/*∂P*, where ω is Raman peak shift. Moreover, χ is determined by κ, *D*, *R*, laser absorption coefficient and temperature coefficient of Raman shift. According to the 5 heating states in Figure 8, 3 normalized RSC were obtained: Θ*n* = <sup>χ</sup>*cw*,*<sup>n</sup>*/(<sup>χ</sup>*ps*<sup>1</sup> − <sup>χ</sup>*ps*<sup>2</sup>), *n* = 1, 2, 3. The error caused by laser absorption, Raman temperature coefficient were eliminated. Meanwhile, the heat accumulation effect was removed by the difference of the heating between the 2 objectives (50<sup>×</sup>, 100×) under picosecond pulse laser. Then, a 3D numerical model was employed to determine κ, *D* and *R*. Figure 9 shows the evolution of the distribution of Ω(<sup>κ</sup>, *D*, *R*). Yuan et al. measured that the in-plane thermal conductivity spans from 31.0 to 76.2 W/(m·K) of 2D few layers MoS2 samples (thickness ranging from 2.4 nm to 37.8 nm) supported on a glass substrate by 5 state picosecond ET-Raman method.

**Figure 8.** The Schematic diagram for mechanism of five-state energy transport state-resolved Raman (ET-Raman) technique. (**a**) The generation, diffusion, and recombination of the hot carrier in MoS2 upon laser irradiating. (**b**,**<sup>c</sup>**) Transient state heating using picosecond laser heating under 50× and 100× objective lenses. (**d**–**f**) Steady state heating using CW laser with 20×, 50×, and 100× objective lenses. Reproduced from Ref. [51] with permission from the Royal Society of Chemistry. Copyright 2018, Royal Society of Chemistry.

**Figure 9.** The evolution of distribution of <sup>Ω</sup>(<sup>κ</sup>, *D*, *<sup>R</sup>*). (**a**) <sup>Ω</sup>(<sup>κ</sup>, *D*, *R*) ≤ 0.65; (**b**) <sup>Ω</sup>(<sup>κ</sup>, *D*, *R*) ≥ 0.80; (**c**) <sup>Ω</sup>(<sup>κ</sup>, *D*, *R*) ≥ 0.95; (**d**) <sup>Ω</sup>(<sup>κ</sup>, *D*, *R*) = 1.0. Reproduced from Ref. [51] with permission from the Royal Society of Chemistry.

Due to the short pulse interval, the picosecond laser, which would generate heat accumulation in the suspended structure, is replaced by a nanosecond laser. As shown in Figure 10, Zobeiri et al. measured the κ and *D* of suspended WS2 by constructing 3 heating states with a continuous laser and a nanosecond pulse laser [52]. The influence of κ and *D* can be distinguished by changing the size of the heating area with a different objective lens. Similarly, 3 RSCs were defined: ψ*CW*, ψ*ns*<sup>20</sup> and ψ*ns*100. Two normalized RSCs were further defined: Θ20 = ψ*ns*20/ψ*CW* and Θ100 = ψ*ns*100/ψ*CW*. Theoretical values Θ under different κ and *D* values were obtained by temperature rise simulation under 3 states, and the matching κ and *D* were obtained by comparing with the experimental values. The thermal conductivity of suspended WS2 was observed to increase from 15.1 to 38.8 W/(m·K) as the sample thickness increased from 13 nm to 107 nm with nanosecond ET-Raman technique.

**Figure 10.** (**<sup>a</sup>**,**b**) Schematic diagram of suspended WS2 illuminated by continuous and nanosecond lasers. (**<sup>c</sup>**,**d**) Energy transport states are constructed by continuous and nanosecond lasers in the temporal and spatial domain. (**<sup>e</sup>**–**g**) Thermal diffusion length, laser radius, and carrier diffusion length under three states. Reprinted with permission from Ref. [52]. Copyright 2019, Elsevier.

Since the Raman signal comes from optical phonons (OPs), but the heat transfer in the sample is related to acoustic phonons (APs). Considering the measurement error caused by ignoring the temperature difference between the 2 phonons, Wang et al. developed 6 heating states nanosecond ET-Raman technique by changing the objective lens: 3 steady states and 3 transient states, and realized the measurement of intrinsic *κ* of MoS2 and MoSe2 nanofilms and phonon coupling factors [53]. As shown in Figure 11a, the Raman spectrum reflects the temperature rise (Δ*Tm*) of OPs, which is the sum of the temperature difference (Δ*TOA*) between OPs and APs and the temperature rise (Δ*TAP*) of APs. Figure 11b shows that the Δ*T*OA decreases to zero faster than Δ*T*AP, which means that phonon coupling between OPs and APs is negligible when the laser spot is very large. Figure 11c shows the Δ*T*m with different laser spot radius, which can be written as:

$$
\Delta T\_{\rm m} = \Delta T\_{\rm OA} + \Delta T\_{\rm AP} \propto \left. \text{Ar}\_{0}^{-2} + f(\kappa) \cdot r\_{0}^{-\eta} \left( \text{n} < 2 \right) \right. \tag{21}
$$

where *r*0 is the laser spot radius, and *f*(κ) is a function of thermal conductivity κ. Based on this more accurate temperature rise fitting process, the intrinsic thermal conductivity of the sample is approximately extracted.

In order to compare the different experimental methods for measuring the thermal conductivity of 2D materials much more conveniently, these methods are summarized in Table 1. The thermal conductivity values of 2D materials with a similar thickness measured by different experimental methods were quite different, which may be attributed to the differences in sample quality and different measurement methods.

**Figure 11.** (**a**) Energy transport process among different energy carriers in suspended 2D materials under laser irradiation. (**b**) The temperature difference between OP and AP against laser spot size. (**c**) Acquisition of thermal conductivity of 2D materials and coupling coefficient between OP and AP. Reprinted with permission from Ref. [53]. Copyright 2020, John Wiley and Sons.


**Table 1.** Application and comparison of various experimental methods.

1 It is measured by the variation of TDTR: FDTR.

#### **4. Analysis of Factors Affecting Thermal Conductivity of 2D Materials**

The thermal conductivity of 2D materials is affected by many factors, such as length, thickness, temperature, substrate, strain, and so on. These factors will affect the process of phonon transmission and scattering and further affect the thermal conductivity.

#### *4.1. Size Effect*

Unlike bulk materials, the thermal conductivity of nm-thick 2D materials usually exhibits an abnormal size dependence. In Zhang's work [59], the in-plane thermal conductivity of h-BCN monolayer calculated by NEMD increases with sample length increasing from 10 nm to 250 nm.

The factor for the size-dependent thermal conductivity originates in phonon scatterings at the sample boundaries. When the phonon mean free path ( *λ*) is larger than the length (l) of the system, heat transfer is ballistic. Certain phonon modes can transmit from the heat-source to the heat-sink without scattering. When *l* > *λ*, phonon scattering is suppressed. Therefore the calculated κ results change with length *l* on small scales. In order to extract the thermal conductivity κ∞ in an infinitely long system, Schelling et al. [60] proposed an extrapolation formula:

$$\frac{1}{\kappa(l)} = \frac{1}{\kappa\_{\infty}} (1 + \frac{\lambda}{l}) \tag{22}$$

#### *4.2. Thickness Effect*

The thermal conductivity of 2D materials is also thickness-dependent. In the work of Smith et al. [61], the thermal conductivity of BP was observed to increase with the thickness increasing from 10 to 1000 nm. However, when the thickness is reduced to less than 10 nm, the thickness dependence of thermal conductivity may show the opposite trend. Yuan et al. [51] reported that the thermal conductivity of 1 to 10 layers decreases with the increase of layers. All these are related to different phonon scattering modes. In monolayer materials, the thermal conductivity is mostly affected by the boundary scattering. Moreover, Umklapp scattering is quenched. Nevertheless, Umklapp scattering has a more significant effect in thicker materials with long phonon mean free path, and the boundary scattering effect is weak which leads to low thermal conductivity.

#### *4.3. Temperature Effect*

Temperature, which can directly affect the thermal performances and cause adverse effects on the structural stability of 2D materials, is also an important factor affecting κ. Hong et al. [62] studied κ of phosphorene/graphene under different temperatures using NEMD. The result showed that the κ of phosphorene and graphene decreased with the increase of temperature, which was as expected for phonon-dominated crystalline materials. As the system temperature increases, more high-frequency phonons are activated, which accelerates heat conduction. Meanwhile, high temperature also promotes Umklapp scattering, which suppresses phonon transmission. The strong scattering effect plays a leading role in the process of heat transfer and eventually leads to the decrease of thermal conductivity. The results show that the maximum reduction of thermal conductivity κ of phosphorene and graphene from 100 K to 400 K is, respectively, 64%, 58%, 11%, and 13%. The calculated thermal conductivity is inversely proportional to the temperature, indicating that the Umklapp scattering is dominant in the temperature range.

#### *4.4. Other Influence Factors*

As the main heat carrier in 2D materials, the propagation of phonons can be adjusted by many other factors, such as lattice deformation caused by strain, substrate coupling, isotope-engineering, which leads to the change of thermal conductivity. Zhang et al. reported that a small strain has a positive effect on the heat conduction of monolayer h-BCN [59]. With further stretching, the thermal conductivity of h-BCN monolayer begins

to decrease. Chen et al. reported that the thermal conductivity of SLG supported on amorphous SiO2 substrate decreased by 40% compared with suspended SLG structure. Through spectral energy density (SED) analysis, it was found that substrate coupling inhibits the thermal transmission of ZA phonons, resulting in a significant reduction in thermal conductivity [63]. In addition, isotope engineering can also affect the thermal conductivity of 2D materials. Through photothermal Raman measurement, Li et al. reported that the in-plane thermal conductivity of isotopic pure 100MoS2 monolayer was 50% higher than that synthesized from naturally abundant isotope mixtures. They attribute this to the former having fewer defects, which reduces phonon-defect scattering [64].

#### **5. Conclusions and Outlook**

In this paper, we systematically introduce the theoretical and experimental methods for the thermal conductivity measurement of 2D materials. The basic principles, advantages, and disadvantages are discussed in detail. Furthermore, some factors (size, temperature, thickness, strain, substrate, and isotope-engineering) that affect the thermal conductivity of 2D materials are also introduced. Based on the thorough analysis, there are many works to conduct to further develop the theoretical and experimental methods.

For the theoretical methods, the accuracy can be further improved by taking more parameters of actual materials into consideration. For example, the growth of 2D materials obtained by chemical vapor deposition (CVD) is controlled by macro physical conditions and parameters, such as partial pressures of each gas in the CVD environment, substrate, defects, furnace configuration, temperature conditions, and gas-phase reactions. One idea is to use growth kinetics and parameter settings to establish the growth model for describing the growth mechanism of the material, which makes the established model consistent with the actual growth model. Netto et al. have reported the continuous growth process of CVD diamond films using timedependent Monte Carlo algorithm with the chemical reaction mechanism [65]. Recently, due to the high calculation requirements for theoretical methods, machine learning has been employed to accelerate the estimation of material thermal conductivity while ensuring the accuracy of measurement results. Mortazavi et al. [66] employed machine-learning interatomic potentials (MLIPs) trained over short ab initio molecular dynamics (AIMD) trajectories instead of density functional theory (DFT) calculation to evaluate anharmonic interatomic force constants, examining the thermal conductivity conveniently, efficiently, and accurately.

For experimental methods, there is also a lot of work to conduct. For suspended micro bridge devices, the contact thermal resistance is an important factor affecting the accuracy of measurement results, which is quantified in the subsequent development of electron beam self-heating method [67]. Furthermore, the suspended device can combine with TDTR for an ultrafast heat pump and probe. In this way, the influence of contact thermal resistance can be eliminated. Moreover, the non-diffusion heat transfer in 2D materials can be characterized. For the 3 ω method, in order to realize its application in measuring the thermal conductivity of 2D materials, the fabrication of a metal electrode with high quality and the signal extraction of a phase-locked amplifier should be considered. Compared with other methods, the Raman method is more widely used in the measurement of 2D material thermal conductivity. However, there is still room for further improvement. First, higher spectral resolution means more accurate temperature measurement. For Raman spectrometer, the higher the grating line density, the higher the corresponding spectral resolution. The grating line density of the commonly used Raman spectrometer ranges from 300 g/mm to 1800 g/mm. If higher density gratings, such as 2400 g/mm and 3600 g/mm, are used, the temperature measurement accuracy will be improved accordingly. Besides, the spatial modulation of laser spot can be considered more. Nowadays, the modulation of laser spot size is realized by using objective lenses with different magnification. In addition, the shape of the laser spot can also be modulated to measure the thermal conductivity of anisotropic 2D materials. Moreover, aside from temporal and spatial modulation, the excitation energy can also be modulated by lasers with different wavelengths. There are few reports on the thermal conductivity of 2D materials using larger wavelength lasers such

as 660 nm laser in the visible light band due to its long exposure time and low excitation efficiency. However, at the same time, the long-wavelength laser also has some advantages, such as reducing fluorescence interference and not easily damaging the sample. Therefore, its application in 2D heat transfer measurement is expected.

**Author Contributions:** H.D. conceived and wrote the review; R.W. conceived and revised the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China (52005367 for Ridong Wang) and National Key R&D Program of China (2020YFC2004600 for Ridong Wang).

**Institutional Review Board Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.
