Non-local thermal transport modeling using the thermal distributor

Advances in technology and experimental studies beyond micro-scale dimensions of materials require new insights into theoretical models that had been developed initially based on the continuum transport theories. The Peierls formulation of thermal transport in solids (the Peierls-Boltzmann equation, PBE [1]) is based on the quasiparticle picture of phonons. The temperature gradient, $\nabla T$, enters the PBE as a driving force. At macroscopic scales and in steady-state, the PBE leads to the Fourier’s law, $\vec{J}=-\kappa(T_{0})\nabla T$, where $\kappa(T_{0})$ is the thermal conductivity at the background temperature, $T_{0}$, and $\vec{J}$ is the heat current density. Small variations of the temperature gradient are ignored. This is called the diffusive regime.

Early experiments on thermal transport in submicron devices [2, 3, 4, 5, 6] showed that the temperature gradient is not constant, but varies on length scales shorter than or comparable to the mean free path (MFP) of phonons. Often analysis of experiment suggests a version of Fourier’s law using an “effective thermal conductivity.” The experiments indicate a non-diffusive regime, with a non-local relation between heat current and temperature gradient. Fourier’s law requires generalization, and Boltzmann theory does this well  [4, 7, 8, 9, 10, 11]. The Boltzmann equation describes the evolution of the phonon distribution function $N_{Q}(r,t)$. When phonons are driven away from equilibrium by local power insertion, it is necessary to add a new term $(dN_{Q}/dt)_{\rm ext}$ describing the power $P_{Q}(r,t)$ added to phonon mode $Q$. The local temperature $T(r,t)=T_{0}+\Delta T(r,t)$ is an ultimate goal, but is not needed to find the non-equilibrium distribution. The inserted power $P_{Q}(r,t)$ is enough to determine $N_{Q}$ and the corresponding heat current $\vec{J}(r,t)$. The local temperature deviation $\Delta T(r,t)$ is an important measure of the behavior of the system, but Boltzmann theory does not contain a definition of $\Delta T(r,t)$; it is necessary to choose a definition. The correct definition is that $C\Delta T(r,t)=\Delta E(r,t)$, where $\Delta E(r,t)$ is the deviation of the total non-equilibrium phonon energy of the system when it is driven away from the equilibrium state at the background temperature $T_{0}$, and $C$ is the specific heat. Unfortunately, when the Boltzmann scattering operator is approximated by its relaxation time approximation (or RTA), an alternate and less physical definition is necessary to restore the energy conservation that is broken by RTA. In this paper, we find $\Delta T(r,t)$ by solving the linearized PBE (or LPBE) using the full scattering operator and the correct definition, and compare it with the RTA version.

Solving the PBE requires a matrix inversion, which is often avoided by using the relaxation time approximation,

$$\left(\frac{\partial N_{Q}}{\partial t}\right)_{\rm scatt}^{\rm RTA}=-\frac{N_{Q}(\vec{r},t)-n_{Q}(T(\vec{r},t))}{\tau_{Q}}.$$

(1)

Here $Q=(\vec{q},s)$ labels phonon modes: $\vec{q}$ is the wavevector, and $s$ is the branch index. The Bose-Einstein distribution $n_{Q}(T(\vec{r},t))$ is evaluated at the local temperature $T(\vec{r},t)$. The phonon relaxation rate $1/\tau_{Q}$ is evaluated using the Fermi golden rule for anharmonic three- phonon scatterings. It is also the diagonal part of the linearized scattering operator, $\hat{S}^{0}$,

$$\left(\frac{\partial N_{Q}}{\partial t}\right)_{\rm scatt}^{\rm LPBE}=-\sum_{Q^{\prime}}S_{QQ^{\prime}}^{0}(N_{Q^{\prime}}-n_{Q^{\prime}}).$$

(2)

In this version labeled with superscript 0, $1/\tau_{Q}=S_{QQ}^{0}$. The correctly linearized operator $\hat{S}^{0}$ is non-Hermitian. For numerical inversion, it is preferable instead to define [12],

	$\displaystyle N_{Q}(\vec{r},t)$	$\displaystyle\equiv$	$\displaystyle n_{Q}(T(\vec{r},t))+$		(3)
		$\displaystyle+$	$\displaystyle n_{Q}^{0}(T_{0})(n_{Q}^{0}(T_{0})+1)\phi_{Q}(\vec{r},t)$

$$\left(\frac{\partial N_{Q}}{\partial t}\right)_{\rm scatt}^{\rm LPBE}=-\sum_{Q^{\prime}}S_{QQ^{\prime}}\phi_{Q^{\prime}}.$$

(4)

where the $n_{Q}^{0}(T_{0})$ is the Bose-Einstein distribution at equilibrium temperature $T_{0}$. Then the operator $\hat{S}$ is Hermitian, and the diagonal element is

$$S_{QQ}=n_{Q}^{0}(T_{0})[n_{Q}^{0}(T_{0})+1]/\tau_{Q}.$$

(5)

The spatially homogeneous PBE driven by a constant $\nabla T$ has been solved by inversion of this Hermitian operator $\hat{S}$ [13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. For spatially inhomogeneous situations, the LPBE (in Fourier space $(\vec{k},\eta)$ rather than coordinate space $(\vec{r},t)$) requires much more difficult inversion of the non-Hermitian operator $S_{QQ^{\prime}}+i(\vec{k}\cdot\vec{v}_{Q}-\eta)\delta_{QQ^{\prime}}$ where $\vec{v}_{Q}$ is the velocity of the phonon mode $Q$. The difficult inversion is avoided by using the RTA approximation $S_{QQ^{\prime}}\rightarrow\delta_{QQ^{\prime}}n_{Q}(n_{Q}+1)/\tau_{Q}$ [23, 24, 25, 26]. Recently inversions with the correct scattering operator for inhomogeneous transport have been done [10, 27].

In [8], the authors introduced a new concept called thermal susceptibility, inspired by the definition of the electrical susceptibility. Thermal susceptibility relates the temperature deviation at $(\vec{r},t)$ to the heat insertion at $(\vec{r}^{\ \prime},t^{\prime})$. This study aims to investigate the capability of the thermal distributor function $\Theta$, which is a redefined version of the thermal susceptibility function:

$$\Theta(\vec{r}-\vec{r}^{\ \prime},t-t^{\prime})\equiv\frac{\delta T(\vec{r},t)}{\delta P(\vec{r}^{\ \prime},t^{\prime})}$$

(6)

for the analysis of non-local thermal transport. Thus the temperature deviation is obtained,

ΔT(→r,t) = 1V∫d→r^  ′ ∫_-∞^t dt^′Θ(→r-→r^  ′,t-t^′) P(→r^  ′,t^′),

(7)

where $V$ is the sample volume. In reciprocal space, $\Delta T(\vec{k},\eta)=\Theta(\vec{k},\eta)P(\vec{k},\eta)$. We apply our analysis to graphene, depicted in Fig. 1.

Because graphene is a two-dimensional crystal, the vectors $\vec{r}$ and $\vec{k}$ are two-dimensional. Because the heat source and sink are parallel to the $y$ axis, the relevant wavevector is $\vec{k}=(k_{x},0)$. Heat current density $\vec{J}_{\rm 2D}$ has units W/m; the more familiar unit in 3D is W/m². Thermal conductivity in 2D has units W/K; in order to compare with 3D, it is conventional to choose the somewhat arbitrary “thickness” for graphene to be $h=3.4$ Å. This paper uses 3D units for current density $J=J_{\rm 2D}/h$, input power $P=P_{\rm 2D}/h$, energy density $U=U_{\rm 2D}/h$, and specific heat $C=C_{\rm 2D}/h$. Then $\kappa$ has conventional units W/mK, and $\Theta$ has units Km³/W.

The measured thermal conductivity of graphene (in 3D units) is reported to lie in the range of 2600 to 5300 W/mK at room temperature [28, 29]. The theoretical thermal conductivity of pristine infinite-size graphene at room temperature is reported in the range of 2800-4300 W/mK [13, 19, 30]. In devices smaller than mean free paths $\Lambda_{Q}$ of important phonons, the measured heat current divided by an approximate measurement of temperature gradient gives an “effective thermal conductivity” ($\kappa_{\rm eff}$) of smaller value. For phonons with bulk $\Lambda_{Q}=|v_{Qx}|\tau_{Q}$ ($\vec{v}_{Q}$ is phonon group velocity) greater than device size $L$, the contribution to $\kappa_{\rm eff}$ is reduced from $C_{Q}|v_{Q,x}|\Lambda_{Q}$ to $C_{Q}|v_{Qx}|L$, where $C_{Q}$ is the contribution of mode $Q$ to the specific heat. We will describe this effect using the thermal distributor function [23].

Under the assumptions of well-defined quasiparticles, the PBE in a crystalline solid is:

$$\frac{dN_{Q}}{dt}=\left(\frac{\partial N_{Q}}{\partial t}\right)_{\rm drift}+\left(\frac{\partial N_{Q}}{\partial t}\right)_{\rm scat}+\left(\frac{\partial N_{Q}}{\partial t}\right)_{\rm ext},$$

(8)

The first term on the right-hand side of Eq. 8 is the change of $N_{Q}$ caused by phonon drift in the distribution gradient:

$$\left(\frac{\partial N_{Q}}{\partial t}\right)_{\rm drift}=-\vec{v_{Q}}.\vec{\nabla}_{\vec{r}}N_{Q}$$

(9)

where $\vec{\nabla}_{\vec{r}}$ is the spatial gradient. The second term in Eq. 8 contains all scattering processes in the crystal. The term from anharmonic three-phonon scatterings ($Q\rightarrow Q^{\prime}+Q^{\prime\prime},Q+Q^{\prime}\rightarrow Q^{\prime\prime}$) can be found from the Fermi golden rule [31]:

$$\begin{split}\left(\frac{\partial N_{Q}}{\partial t}\right)_{\rm scat}&=\frac{\pi\hbar}{16N_{q}}\sum_{Q^{\prime}Q^{\prime\prime}}\left|V_{QQ^{\prime}Q^{\prime\prime}}\right|^{2}\\ &\frac{1}{2}\bigg{\{}N_{Q}(N_{Q^{\prime}}+1)(N_{Q^{\prime\prime}}+1)-(N_{Q}+1)N_{Q^{\prime}}N_{Q^{\prime\prime}}\bigg{\}}\delta(\omega_{Q}-\omega_{Q^{\prime}}-\omega_{Q^{\prime\prime}})\\ &+\bigg{\{}N_{Q}N_{Q^{\prime}}(N_{Q^{\prime\prime}}+1)-(N_{Q}+1)(N_{Q^{\prime}}+1)N_{Q^{\prime\prime}}\bigg{\}}\delta(\omega_{Q}+\omega_{Q^{\prime}}-\omega_{Q^{\prime\prime}}),\end{split}$$

(10)

where $\omega_{Q}$ is the phonon frequency and $N_{q}$ is a number of wavevectors in the Brillouin zone ¹¹1In the simulations, we use $N_{q}=120\times 120$ $q$-points to sample the Brillouin zone of graphene and for delta functions in Eq. (10) we use Gaussian broadening function $\delta(x)=\exp{-(x/\sigma)^{2}}/\sqrt{\pi}\sigma$ with $\hbar\sigma=0.9$ meV.. $V_{QQ^{\prime}Q^{\prime\prime}}$ is the matrix element of the three-phonon process, given by,

$$\begin{split}V_{QQ^{\prime}Q^{\prime\prime}}=\sum_{MLnml}\sum_{\alpha\beta\gamma}\frac{\epsilon_{Q}^{n\alpha}\epsilon_{Q^{\prime}}^{m\beta}\epsilon_{Q^{\prime\prime}}^{l\gamma}}{\sqrt{M_{\rm c}^{3}\omega_{Q}\omega_{Q^{\prime}}\omega_{Q^{\prime\prime}}}}\Psi_{\alpha\beta\gamma}(0n,Mm,Ll)e^{i\bm{q^{\prime}.R_{M}}}e^{i\bm{q^{\prime\prime}.R_{L}}}\delta(q+q^{\prime}+q^{\prime\prime},G),\end{split}$$

(11)

where $\Psi_{\alpha\beta\gamma}(0n,Mm,Ll)$ is the third derivative of the crystal potential by displacements of atoms in positions $(n,m,l)$ inside the unit cells $(0,M,L)$. The supercell with index $0$ is the central unit cell. $\epsilon_{Q}^{n\alpha}$ is the $\alpha^{\rm th}$ Cartesian component of the polarization vector of mode $Q$ at atom $n$, and $M_{\rm c}$ is the mass of a carbon atom. The Kronecker delta ensures the conservation of the lattice momentum, where $G$ is a reciprocal lattice vector. The local equilibrium phonon population $n_{Q}=n_{Q}(T(\vec{r},t))$ depends implicitly on the space and time through its explicit dependence on the temperature $T(\vec{r},t)$. For small deviations from equilibrium, expand the phonon population $N_{Q}$ as in Eq. 3 to first order in $\phi_{Q}$. The anharmonic scattering matrix $\hat{S}$ (Eq. 4) then has diagonal and off-diagonal elements,

	$\displaystyle S_{QQ}$	$\displaystyle=$	$\displaystyle 1/\tau_{Q}=2\pi\hbar\sum_{Q^{\prime}Q^{\prime\prime}}\left|V_{QQ^{\prime}Q^{\prime\prime}}\right|^{2}$
			$\displaystyle\bigg{\{}\frac{(n_{Q}^{0}+1)n_{Q^{\prime}}^{0}n_{Q^{\prime\prime}}^{0}}{2}\delta(\omega_{Q}-\omega_{Q^{\prime}}-\omega_{Q^{\prime\prime}})+n_{Q}^{0}n_{Q^{\prime}}^{0}(n_{Q^{\prime\prime}}^{0}+1)\delta(\omega_{Q}+\omega_{Q^{\prime}}-\omega_{Q^{\prime\prime}})\bigg{\}}$
	$\displaystyle S_{QQ^{\prime}}$	$\displaystyle=$	$\displaystyle 2\pi\hbar\sum_{Q^{\prime\prime}}\left|V_{QQ^{\prime}Q^{\prime\prime}}\right|^{2}\bigg{\{}n_{Q}^{0}n_{Q^{\prime}}^{0}(n_{Q^{\prime\prime}}^{0}+1)\delta(\omega_{Q}+\omega_{Q^{\prime}}-\omega_{Q^{\prime\prime}})$		(12)
			$\displaystyle-(n_{Q}^{0}+1)n_{Q^{\prime}}^{0}n_{Q^{\prime\prime}}^{0}\delta(\omega_{Q}-\omega_{Q^{\prime}}-\omega_{Q^{\prime\prime}})-n_{Q}^{0}(n_{Q^{\prime}}^{0}+1)n_{Q^{\prime\prime}}^{0}\delta(\omega_{Q}-\omega_{Q^{\prime}}+\omega_{Q^{\prime\prime}})\bigg{\}}$

This version of the scattering matrix is real-symmetric ($S_{QQ^{\prime}}=S_{Q^{\prime}Q}$) i.e. Hermitian. Each collision conserves phonon energy, which is assured by

$$S_{QQ}\omega_{Q}+\sum_{Q^{\prime},Q^{\prime}\neq Q}S_{QQ^{\prime}}\omega_{Q^{\prime}}=0,\ {\rm or}\ \hat{S}|\omega\rangle=0.$$

(13)

The mode frequency $\omega_{Q}=\langle Q|\omega\rangle$ is an eigenvector, in fact, the only “null eigenvector”, of the linearized Hermitian scattering operator $S_{QQ^{\prime}}=\langle Q|\hat{S}|Q^{\prime}\rangle$.

The last term in Eq. (8) models external heat sources and sinks. The form is usually [33, 34]

$$\left(\frac{\partial N_{Q}}{\partial t}\right)_{\rm ext}=\frac{P_{Q}(\vec{r},t)}{C}\frac{dn_{Q}}{dT},$$

(14)

The heat source/sink $P$, its geometry $(\vec{r},t)$, and its spectral distribution $Q$ determine whether the heat transport is quasiballistic or diffusive. We use the simplest version where $P_{Q}=P$ is independent of $Q$ [35, 36, 8, 34],

$$\left(\frac{\partial N_{Q}}{\partial t}\right)_{\rm ext}=\frac{P}{C}\frac{dn_{Q}}{dT},$$

(15)

where $P$ is the heat power added per unit volume of the system. Each mode gets the same boost $\Delta T$ from $P(\vec{r},t)$. Detailed knowledge of source and sink would cause a $Q$-dependence of $P$ [23, 37], but missing this knowledge, $P_{Q}=P$ is a sensible guess.

In vector-space notation, the LPBE, Eq. 8, is

$$\frac{\partial}{\partial t}[|n\rangle+|n^{0}(n^{0}+1)\phi\rangle]=-\vec{v_{Q}}\vec{\nabla}_{\vec{r}}[|n\rangle+|n^{0}(n^{0}+1)\phi\rangle]-\hat{S}|\phi\rangle+\frac{P(\vec{r},t)}{C}|\frac{dn}{dT}\rangle.$$

(16)

In this notation, the kets (like $|n\rangle$) are vectors in the space of phonon modes, with components $\langle Q|n\rangle=n_{Q}$. The solution $\phi_{Q}(\vec{r},t)$ is found from its Fourier $(\vec{k},\eta)$ representation,

$$\phi_{Q}(\vec{r},t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\eta e^{-i\eta t}\sum_{\vec{k}}e^{i\vec{k}\cdot\vec{r}}\phi(\vec{k},\eta).$$

(17)

The temperature deviation $\Delta T(\vec{r},t)$ and the power input $P(\vec{r},t)$ are also transformed to Fourier space. To simplify the algebra, define a vector $|X\rangle$ and an operator $\hat{W}$.

$$|X\rangle\equiv|n^{0}(n^{0}+1)\hbar\omega\rangle$$

(18)

$$\hat{W}\equiv\hat{S}+i(\vec{k}\cdot\hat{\vec{v}}-\eta\hat{1})\hat{n}^{0}(\hat{n}^{0}+\hat{1})$$

(19)

where $\hat{\vec{v}}$ and $\hat{n}^{0}$ are diagonal in $Q$ space (i.e. $\langle Q|\hat{n}^{0}|Q^{\prime}\rangle=n_{Q}^{0}\delta(Q,Q^{\prime})=\langle Q|n^{0}\rangle\delta(Q,Q^{\prime})$). The LPBE in Fourier space is

$$\begin{split}\hat{W}|\phi\rangle=\frac{1}{k_{B}T^{2}}&\Big{[}-i(\vec{k}\cdot\hat{\vec{v}}-\eta\hat{1})\Delta T(\vec{k},\eta)+\\ &\frac{P(\vec{k},\eta)}{C}\hat{1}\Big{]}|X\rangle\end{split}$$

(20)

Quasiparticle heat transport in solids can be roughly categorized into three regimes: ballistic, quasi-ballistic, and diffusive. In the diffusive regime, where Fourier’s law applies, the channel length of the heat conductors is much longer than the MFPs; phonons experience multiple scattering events, so that a local equilibrium is established, with a temperature gradient that is constant everywhere except in the small region close to heat sources and sinks. Thermal transport is ballistic when the channel length is comparable to or less than phonon MFPs. In this case, thermal energy is mainly dissipated near the heat sources. The temperature gradient becomes thermally inhomogeneous, and heat current has a non-local relation to temperature. When the channel length is similar to the averaged phonon MFP, some phonons travel ballistically and others diffusely; heat transfer is “quasi-ballistic”. The wide range of phonon MFPs in graphene  [38] makes it hard to differentiate transport regimes. Moreover, the heat transport regime in a given device is temperature-dependent since the phonon MFPs depend on temperature.

The thermal distributor $\Theta(\vec{r})$, Eq. 6, simplifies the description of heat transport in different regimes. The spatial variation of temperature is given by

$$T(\vec{r})=\sum_{\vec{k}}T(\vec{k})e^{i\vec{k}\cdot\vec{r}}=\sum_{\vec{k}}\Theta(\vec{k})P(\vec{k})e^{i\vec{k}\cdot\vec{r}}.$$

(30)

The Fourier transform $\Theta(\vec{k})$ is related to the non-local generalization of the bulk thermal conductivity $\kappa_{\rm bulk}={\rm lim}_{\vec{k}\rightarrow 0}\kappa(\vec{k})$ by Eq. 26. First, we calculate LPBE and RTA versions of $\Theta(\vec{k})$ of graphene from Eq. (27) and Eq. (29), using the modified Tersoff potential, with the parameters given in ref. 39, to model the crystal potential. The scattering rates are shown in Fig. 2 at $T=300$ K.

Note that there are numerical challenges in calculating the scattering rates of phonons using the Gaussian broadening [40]. Following Ref. [13, 31], we apply linear interpolation of the ZA scattering rates with phonon energy, as shown by the blue circles for phonon energies below 0.3 THz in Fig. 2b. This linear scaling is explained by noting that the bending energy is given by $E_{b}=\int d^{2}r\kappa_{b}/(2R^{2})$, where $\kappa_{b}=2.1$ eV [41] is the bending stiffness, and $R$ is the radius of curvature given by $1/R=d^{2}z/dx^{2}$. Applying the Bloch theorem for a discrete atomic model with a lattice constant $a$, one can show that bending energy for mode $q$ is given by $E_{b}=8\kappa_{b}A_{c}z_{q}^{2}\sin^{4}(qa/2)/a^{4}\approx A_{c}\kappa_{b}z_{q}^{2}q^{4}/2$, where $A_{c}$ is the area per atom. By comparing it to the harmonic oscillator potential energy $m\omega_{q}^{2}z_{q}^{2}/2$, one can obtain the expected result for flexural phonon frequency: $\omega_{q}=q^{2}\sqrt{\kappa_{b}A_{c}/m}$. The third order anharmonicity can be introduced by coupling the ZA mode to an LA mode by modifying the bending stiffness $\kappa_{b}=\kappa_{b0}-\alpha_{b}(x_{i+1}-x_{i-1})$. After applying the Bloch theorem, the third-order anharmonic potential for mode $q$ in the small $q$-limit becomes $H_{3}=\sum_{q_{1}}i\alpha_{b}A_{c}z_{q}z_{q_{1}}x_{-q-q_{1}}q^{2}q_{1}^{2}(q+q_{1})a$. Using second-quantized amplitudes for the phonon displacements $z_{q}$ and $x_{q}$, one can show that $V_{qq_{1}}\sim q^{2}q_{1}^{2}(q+q_{1})(\omega_{q}^{ZA}\omega_{q_{1}}^{ZA}\omega_{q+q_{1}}^{LA})^{-1/2}$. The $q$-ZA phonon scattering with $q_{1}$-ZA phonon into an $(q+q_{1})$-LA phonon has a rate of $1/\tau^{ZA}_{q}\sim\sum_{q_{1}}|V_{q,q_{1}}|^{2}n^{0,ZA}_{q_{1}}n^{0,LA}_{q+q_{1}}\delta(\omega^{ZA}_{q}+\omega^{ZA}_{q_{1}}-\omega^{LA}_{q+q_{1}})$. The $\delta$-function ensures that $q_{1}\sim v_{s}\sqrt{m/(\kappa_{b}A_{c})}$, where $v_{s}$ is the sound velocity. Therefore, the interpolation function $1/\tau^{ZA}_{q}\sim q^{2}\sim\omega^{ZA}_{q}$ used in Fig. 2b can be justified.

Fig. 3 shows the calculated $\Theta(k)$ for graphene using LPBE and RTA approximations. The width of the sample is much broader than the MFP; this allows a one-dimensional treatment. The spatial variation $\vec{r}=(x,0)$ is only along $\hat{x}$, parallel to the channel, so the spatial Fourier variable is $\vec{k}=(k,0)$. The spatial resolution of $T(x)$ at small distances $x$ requires values of $\Theta(k)$ at correspondingly large $k\sim 2\pi/x$ (see Eq. 30). The range of $k$ in Fig. 3 corresponds to distances from a few nanometers to a meter-long channel length. A lower limit of the length scale is imposed by the validity of the quasiparticle picture of phonons used in the PBE formalism  [26]. Note that the thermal distributor diverges for both LPBE and RTA solutions when $k\to 0$. According to Eq. (26), for a finite thermal conductivity in a diffusive regime, $\Theta(k)$ must diverge in $k\to 0$ limit as $\Theta(k)\sim\frac{1}{k^{2}}$. Curve-fitting shows that the calculated $\Theta(k)$ for both versions agrees with $\frac{1}{k^{2}}$ very accurately in the small $k$ limit, namely:

$$\begin{split}&\Theta_{LBPE}(k)=\frac{2.4\times 10^{-4}\ {\rm Km/W}}{k^{2}}\\ &\Theta_{RTA}(k)=\frac{1.51\times 10^{-3}\ {\rm Km/W}}{k^{2}}\end{split}$$

(31)

Using Eq. (26), this corresponds to bulk thermal conductivities 4145 W/mK and 662 W/mK for LPBE and RTA, respectively. The large discrepancy between LPBE and RTA values of thermal conductivities of graphene has been noticed previously [13].

The thermal conductivities evaluated according to Eq. (25) for LPBE and Eqs. 26, 29 for RTA are shown in Fig. 4. The sharp fall-off of the thermal conductivity in Fig. 4 indicates the ballistic-to-diffusive crossover; it happens at larger $k$, and, therefore, smaller characteristic lengths in the RTA treatment than in the correct LPBE treatment.

Now we discuss thermal conduction in the geometry of Fig. 1 using the PBE to evaluate the thermal distributor. We can calculate the temperature profile for any pattern of heat input and removal from this response function, provided the graphene sample has one-dimensional periodicity. Energy conservation requires $dJ/dx=-P(x)$. Steady-state ($\eta=0$) requires equal external heat addition and removal. $L_{\rm ch}$ is the length of the channel between the two heat reservoirs, source and sink; each of length $L_{\rm s}$. For simplicity, our heat input has odd symmetry ($P(-x)=-P(x)$) around $x=0$, so $J$ is even in $x$. The period of the supercell is $L_{d}=2(L_{\rm ch}+L_{\rm s})$, which determines the shortest non-zero wavevector i.e. $k_{min}=2\pi/L_{d}$. A large $L_{d}$ value allows a fine Fourier mesh to describe nanoscale physics, such as the ballistic-to-diffusive crossover. Our reservoirs are ideal thermal baths, with zero interfacial thermal resistance between the channel and reservoirs. Note that experimental thermal conductivity measurements are often done using periodic structures such as periodic metallic gratings or transient thermal gratings  [42, 43]. Our periodic geometry of Fig. 1 works for both periodic structures and single-channel devices. In the latter case, it is necessary to make $L_{\rm s}$ larger than phonon MFPs.

Figure 1 shows how heat insertion and removal $P(x)$ is distributed uniformly (with magnitude $P_{0}$) over lengths $L_{\rm s}$ on either side of the sample (or channel) of length $L_{\rm ch}$. Energy conservation $dJ/dx=P$ then gives heat current density $J=P_{0}L_{\rm s}/2$ in the channels. Figure 5 shows the resulting $\Delta T$ profiles in three devices with channel lengths spanning from ballistic to diffusive regimes. In the ballistic device, Fig. 5a, both $\Theta_{RTA}$ and $\Theta_{LBTE}$ predict similar values for temperature profiles in the channel and in the source/sink reservoirs. This behavior can be understood from the two thermal distributors being similar in magnitude in the large $k$-limit, as shown in Fig. 3. The temperature profile of Fig. 5a clearly shows the source regions are hotter than the channel, which is a characteristic of non-diffusive thermal transport [34, 33, 44]. The long MFP phonons (ZA phonons) dissipate the thermal energy in the source regions while flying in the channel ballistically. As a result, the temperature is higher in the heat region. This observation manifests the nonlocal effect.

For the larger channel lengths in Fig. 5b and 5c, RTA predictions for the temperature are significantly higher than LPBE predictions in agreement with Fig. 3 which shows that at smaller $k$ (corresponding to larger distances) $\Theta_{RTA}$ is greater than $\Theta_{LPBE}(k)$. In the quasi-ballistic regime of Fig. 5b, there is significant non-linearity near the heat source/sink. Although in this regime both long and short MFP phonons contribute to thermal transport, the share of thermal energy carried by short MFP phonons is almost negligible.

Non-local heat transport is often described by an “effective thermal conductivity” $\kappa_{\rm eff}$. The definition varies depending on the experiment. For example, experimental studies such as time-domain thermoreflectance [45] measure the transient temperature response to a heat pulse. This can be used to extract $\kappa_{\rm eff}$ [42]. Many theoretical studies also have applied a similar procedure.

For our steady-state thermal transport computations, there are three sensible definitions of $\kappa_{\rm eff}$, shown in table 1. The first of these defines $\kappa_{\rm eff}$ by applying Fourier’s law at the middle of the channel where the curvature of $\Delta T(x)$ is zero. This resembles Fourier’s law in thermally homogeneous structures. We note that in the fully diffusive regime, when both $L_{\rm s}$ and $L_{\rm ch}$ are larger than the MFPs, the $\kappa_{\rm eff,mp}$ and $\kappa_{\rm eff,ch}$ values coincide, and the temperature drop in the channel is $\Delta T_{\rm ch}=J_{0}L_{\rm ch}/\kappa_{\rm eff,mp}$. In this limit, the temperature drop under the contact is equal to $J_{0}L_{\rm s}/(4\kappa_{\rm eff,mp})$, because current $J(x)$ varies linearly with distance from the middle of the contact. Therefore, the maximum temperature variation across the unit cell equals $\Delta T_{\rm max}=\Delta T_{\rm ch}+J_{0}L_{\rm s}/(2\kappa_{\rm eff,mp})$ and effective thermal conductivity $\kappa_{\rm eff,min}=J_{0}(L_{\rm ch}+L_{\rm s})/\Delta T_{\rm max}=\kappa_{\rm eff,mp}(L_{\rm ch}+L_{\rm s})/(L_{\rm ch}+L_{\rm s}/2)$. As can be seen from table 1, the device geometry in Fig. 5c approaches the diffusive limit. However, one should note the fully diffusive limit relationship between $\kappa_{\rm eff,min}$ and $\kappa_{\rm eff,mp}$ is not realized in the devices we considered.

The $\kappa_{\rm eff,mp}$ definition is plotted in Fig. 6 and 7, for both LPBE and RTA versions of the thermal distributor, as a function of the source length and channel length, respectively. When the chosen length $L_{\rm ch}$ in Fig. 6 (or $L_{\rm s}$ in Fig. 7) has a value less than the diffusive length scale ($\sim 1-3\>\mu$m for LPBE and $\sim 100-200\>$ nm for RTA), $\kappa_{\rm eff,mp}$ decreases with decreasing periodicity due to suppression of the ballistic phonon contribution to the heat transport. At large values of $L_{\rm ch}$, thermal conductivity saturates and approaches the diffusive limit when both $L_{\rm ch}$ and $L_{\rm s}$ are large.

Similarly, Fig. 7 shows that $\kappa_{\rm eff,mp}$ saturates with increasing $L_{\rm s}$ and reaches the diffusive limit value at large $L_{\rm ch}$. For small $L_{\rm ch}$ values, the $\kappa_{\rm eff,mp}$ dependence on $L_{\rm s}$ shows ballistic to diffusive crossover as in Fig. 6. Those dependencies happen because of superposing interaction of ballistic phonons in the ballistic channels and recovery of diffusive-like thermal transport. This phenomenon shows the deterministic role of the geometry of the heat source and has already been observed experimentally  [43, 46].

To quantify our results in Fig. 6 and 7, we fit our $\kappa_{\rm eff,mp}$ to the phenomenological ballistic-to-diffusive crossover equations used to describe electrical transport [47]:

	$\displaystyle\kappa_{\rm eff,a}$	$\displaystyle=$	$\displaystyle\kappa_{0}+(\kappa_{\rm dif}-\kappa_{0})\frac{L_{\rm ch}}{L_{\rm ch}+\lambda}\$		(32)
	$\displaystyle\kappa_{\rm eff,b}$	$\displaystyle=$	$\displaystyle\kappa_{0}+(\kappa_{\rm dif}-\kappa_{0})\frac{L_{\rm s}}{L_{\rm s}+\lambda}$		(33)

where $\kappa_{\rm dif}$ is thermal conductivity in the diffusive regime. From the best fits, we find a characteristic length scale $\lambda\sim 1$ $\mu$m ($\sim 90$ nm) for the LPBE (RTA) solution using fits of $\kappa_{\rm eff,a}$ versus $L_{\rm ch}$, and $\lambda\sim 3$ $\mu$m ($\sim 200$ nm) for LPBE (RTA) solution using fits of $\kappa_{\rm eff,b}$ versus $L_{\rm s}$. The value of $\kappa_{0}$ approaches zero only when both $L_{\rm s}$ and $L_{\rm ch}$ are small, as discussed above. Those estimates for the characteristic crossover length scales are consistent with the $k$-dependences of thermal conductivity in Fig. 4. Using a characteristic value of $k_{0}$, when thermal conductivity drops by a factor of two in Fig. 4, we find $2\pi/k_{0}\sim 5$ $\mu$m for LPBE and $2\pi/k_{0}\sim 300$ nm for RTA. Finally, we can estimate an average phonon MFP using the standard expression for thermal conductivity in 2D: $\kappa=Cv\lambda_{\rm ph}/2$, where the heat capacity is $C\times h=4.5\times 10^{-4}{\rm J/Km}^{2}$ (calculated at $T=300$ K) and $v$ is an averaged phonon velocity. Using the velocity of the ZA parabolic band at room temperature $v\approx 10$ km/s, we obtain $\lambda_{\rm ph}=625$ nm for LPBE and $\lambda_{\rm ph}=100$ nm for RTA, consistent with the above estimates for the diffusive to ballistic crossover length scales.

Steady microscale addition and removal of heat energy, together with a scattering of heat carriers, leads to carrier distributions $N$ centered around local equilibrium $n$ with a local temperature $T(x)$. The thermal distributor $\Theta(k)$ contains all relevant information about thermal transport in all the regimes. By computing the thermal distributor $\Theta(k)$ of graphene from the PBE with correct inversion using the full scattering operator, we investigated thermal transport at submicron length scales. At long length scales, heat transport occurs with constant temperature gradients. Non-local effects (where $\nabla T$ varies with $x$) are seen at length scales of the order of or less than mean free paths of phonons. Details of the geometrical structure of the device and heat sources and sinks cause the inhomogeneous temperature profile $T(x)$. The long-wavelength phonons are forcefully scattered in the source and sink regions while flying the channel without appreciable scatterings with other phonons. This causes local $\Delta T(x)$ to be higher near the source/sink. This is often ascribed to the suppression of heat transport by phonons with long mean free paths. The RTA contains these effects but, especially in graphene, overestimates the temperature inhomogeneity needed to drive a heat current.

We acknowledge support from the Vice President for Research and Economic Development (VPRED), SUNY Research Seed Grant Program, and the Center for Computational Research at the University at Buffalo [48].