Phonon Boltzmann equation non-local in space and time: the partial failure of the generalized Fourier law

Transport in bulk systems, with currents small enough to linearize, is given by simple linear response functions, such as the Fourier law for thermal conductivity, $\vec{J}_{\rm th}=-\kappa\vec{\nabla}T$, where the temperature is close to a background temperature $T_{0}$, deviating by a small correction $\vec{r}\cdot\vec{\nabla}T$, with $\vec{\nabla}T$ uniform in space. However, in small systems, the current and temperature need to be studied on distance scales less than the mean free paths of heat carriers. A sensible assumption is that a non-local Fourier law, $\vec{J}_{\rm F}(\vec{r})=-\int d\vec{r}^{\ \prime}\kappa(\vec{r},\vec{r}^{\ \prime})\vec{\nabla}T(\vec{r}^{\ \prime})$ should describe the linear response. However, for non-uniform heat driving, a correction has been noticed recently. An additional, non-Fourier term $\vec{J}_{\rm non-F}(\vec{r})$ appears. This was found recently by Hua and Lindsay Hua et al. (2019); Hua and Lindsay (2020). Hints of this can be found in earlier literature. Boltzmann theory shows how this happens. The derivation given below assumes heat carried by phonons in insulators. Similar results are expected for heat carried by electrons in metals, and also for electrical conduction. The derivation includes time-dependent driving of phonon heat.

The phonon, or Peierls, Boltzmann equation Peierls (1929); Ziman (1960) (PBE), linearized for small deviations from equilibrium, says that the phonon distribution $N_{Q}(\vec{r},t)$ evolves according to

	$\displaystyle\frac{dN_{Q}}{dt}$	$\displaystyle=$	$\displaystyle\left(\frac{dN_{Q}}{dt}\right)_{\rm drift}+\left(\frac{dN_{Q}}{dt}\right)_{\rm coll}+\left(\frac{dN_{Q}}{dt}\right)_{\rm ext},$
			$\displaystyle{\rm where}\ \left(\frac{dN_{Q}}{dt}\right)_{\rm drift}=-\vec{v}_{Q}\cdot\vec{\nabla}_{\vec{r}}N_{Q},$
			$\displaystyle{\rm and}\ \left(\frac{dN_{Q}}{dt}\right)_{\rm coll}=-\sum_{Q^{\prime}}R_{QQ^{\prime}}\Delta N_{Q^{\prime}}(\vec{r},t),$
			$\displaystyle{\rm and}\ \left(\frac{dN_{Q}}{dt}\right)_{\rm ext}=\frac{1}{C(T_{0})}\left(\frac{dn_{Q}}{dT}\right)_{T_{0}}P_{Q}(\vec{r},t).$
					(1)

The symbol $Q=(\vec{q},j)$ enumerates the $N=3nN_{\rm cell}$ eigenstates of the harmonic vibrational sustem. The scattering operator $R_{QQ^{\prime}}$ is linearized for small deviations $\Delta N_{Q}$ from the local equilibrium distribution, a Bose-Einstein distribution $n_{Q}$ at the local temperature $T(\vec{r},t)$; $P_{Q}(\vec{r},t)$ is the rate per unit volume of energy input into mode $Q$, and $T_{0}$ is the background temperature. The specific heat $C(T_{0})$ has units energy divided by temperature$\times$volume, and $P_{Q}$ has units power per volume. The external power input $P_{Q}$ causes mode $Q$ to increase in energy per unit time by an amount which would increase its effective temperature by $d\Delta T_{Q}(\vec{r},t)/dt=P_{Q}(\vec{r},t)/C$. For simple driving, with $P_{Q}$ independent of $Q$, the modal temperature increase $\Delta T_{Q}$ is the same for all modes $Q$.

Non-local Boltzmann theories require a definition of local temperature; the one that works in Boltzmann theory is that the local energy density $u(\vec{r},t)$ is $(1/V)\sum_{Q}\hbar\omega_{Q}n_{Q}(T(\vec{r},t))$. This means that there is no energy in the deviation $\Delta N_{Q}\equiv N_{Q}(\vec{r},t)-n_{Q}(T(\vec{r},t))$. The sum $(1/V)\sum_{Q}\hbar\omega_{Q}\Delta N_{Q}=0$ defines $T(\vec{r},t)$.

Each scattering event conserves phonon energy: $\sum_{Q}\hbar\omega_{Q}(d\Delta N_{Q}/dt)_{\rm scatt}=0$, whether linearized in $\Delta N_{Q}$ or not. In linear approximation, this is equivalent to the sum rule $\sum_{Q}\hbar\omega_{Q}R_{QQ^{\prime}}=0$. The operator $\hat{R}$ (where $R_{QQ^{\prime}}=\langle Q|\hat{R}|Q^{\prime}\rangle$) has a left null eigenvector $\langle\hbar\omega|\hat{R}=0$, where $\hbar\omega_{Q}=\langle\hbar\omega|Q\rangle$. The vector $|\cdot\rangle$ lies in the space of harmonic eigenstates. The states $|Q\rangle$ are a complete orthonormal basis. In the coordinate space enumerated by unit cell $\ell$ at $\vec{R}_{\ell}$, and atom number $n$, and Cartesian coordinate $i$,

$$\langle\ell ni|Q\rangle=\frac{1}{\sqrt{N_{\rm cells}}}e^{i\vec{q}\cdot\vec{R}_{\ell}}\epsilon_{Q}(n,i),$$

(2)

where the polarization vector $\hat{\epsilon}_{Q}$ is normalized, $\sum_{n,i}|\epsilon_{Q}(n,i)|^{2}=1$. But $\hat{R}$ is not symmetric. The corresponding null right eigenvector is $|n(n+1)\hbar\omega\rangle$. This is equivalent to the statement that if the distribution $N_{Q}=n_{Q}(T(\vec{r},t))+\Delta N_{Q}(\vec{r},t)$ consists of shifting the occupancy of every mode $Q$ away from local equilbrium $T(\vec{r},t)$ by $\Delta N_{Q}=(dn_{Q}/dT)\Delta T(\vec{r},t)$ with the same $\Delta T(\vec{r},t)$ for each mode, the shifted distribution is already locally thermalized and doesn’t experience net scattering.

A symmetrized scattering operator $\hat{\Omega}$ is Simons (1960); Krumhansl (1965)

	$\displaystyle\hat{\Omega}$	$\displaystyle=$	$\displaystyle[\hat{n}(\hat{n}+1)]^{-1/2}\hat{R}[\hat{n}(\hat{n}+1)]^{1/2}$
	$\displaystyle\Omega_{QQ^{\prime}}$	$\displaystyle=$	$\displaystyle[n_{Q}(n_{Q}+1)]^{-1/2}R_{QQ^{\prime}}[n_{Q^{\prime}}(n_{Q^{\prime}}+1)]^{1/2}.$		(3)

The operator $\hat{n}$ is defined as $\langle Q|n|Q^{\prime}\rangle=n_{Q}\delta_{QQ^{\prime}}$. The Bose factors $n_{Q}$ in Eq. 3 are taken at the background temperature $T_{0}$, not the local equilibrium temperature $T(\vec{r},t)$. Similarly the operator $\hat{n}$ is defined at $T=T_{0}$ unless explicitly written as $\hat{n}(\vec{r},t)$. The Boltzmann Eq. 1 is now

	$\displaystyle\frac{dN(\vec{r},t)\rangle}{dt}$	$\displaystyle=$	$\displaystyle-\hat{\vec{v}}\cdot\vec{\nabla}_{\vec{r}}|N(\vec{r},t)\rangle$		(4)
		$\displaystyle-$	$\displaystyle\sqrt{\hat{n}(\hat{n}+1)}\hat{\Omega}\frac{1}{\sqrt{\hat{n}(\hat{n}+1)}}|\Delta N(\vec{r},t)\rangle$
		$\displaystyle+$	$\displaystyle\hat{n}(\hat{n}+1)\frac{\hbar\hat{\omega}}{Ck_{B}T^{2}}|P(\vec{r},t)\rangle,$

where $\langle Q|P\rangle=P_{Q}$ The operators $\hat{\vec{v}}$ and $\hat{\omega}$ are diagonal in $Q$-representation, $\langle Q|\hat{\vec{v}}|Q^{\prime}\rangle=\vec{v}_{Q}\delta_{QQ^{\prime}}$ and $\langle Q|\hat{\omega}|Q^{\prime}\rangle=\omega_{Q}\delta_{QQ^{\prime}}$.

Now define from $\Delta N_{Q}$ a deviation function $\Phi_{Q}$,

	$\displaystyle\Delta N_{Q}$	$\displaystyle=$	$\displaystyle N_{Q}-n_{Q}(\vec{r},t)\equiv\sqrt{n_{Q}(n_{Q}+1)}\Phi_{Q}(\vec{r},t)$
	$\displaystyle|\Delta N\rangle$	$\displaystyle=$	$\displaystyle|N\rangle-|n(\vec{r},t)\rangle\equiv\sqrt{\hat{n}(\hat{n}+1)}|\Phi(\vec{r},t)\rangle$		(5)

Then Eq. 4 becomes

$$\frac{\partial|\Phi(\vec{r},t)\rangle}{\partial t}+\hat{\vec{v}}\cdot\vec{\nabla}_{\vec{r}}|\Phi(\vec{r},t)\rangle+\hat{\Omega}|\Phi(\vec{r},t)\rangle=|A(\vec{r},t)\rangle,$$

(6)

	$\displaystyle|A(\vec{r},t)\rangle$	$\displaystyle=$	$\displaystyle-[\hat{n}(\hat{n}+1)]^{-1/2}\left[\left(\frac{\partial}{\partial t}+\hat{\vec{v}}\cdot\vec{\nabla}_{\vec{r}}\right)\Delta T(\vec{r},t)\right]\left|\frac{dn}{dT}\right\rangle$		(7)
		$\displaystyle+$	$\displaystyle[\hat{n}(\hat{n}+1)]^{1/2}\frac{\hbar\hat{\omega}}{Ck_{B}T^{2}}|P(\vec{r},t)\rangle,$

The inhomogeneous part $|A\rangle$ has been separated. It drives the deviation $|\Phi\rangle$ from local equilibrium, $|\Phi\rangle=0$. Because of linearization, the equation is simplified by a Fourier transform,

	$\displaystyle|X(\vec{k},\omega)\rangle$	$\displaystyle=$	$\displaystyle\frac{1}{V}\int^{V}d\vec{r}\int_{-\infty}^{\infty}dt|X(\vec{r},t)\rangle\exp[-i(\vec{k}\cdot\vec{r}-\omega t)]$
	$\displaystyle|X(\vec{r},t)\rangle$	$\displaystyle=$	$\displaystyle\sum_{\vec{k}}\int_{-\infty}^{\infty}d\omega|X(\vec{k},\omega)\rangle\exp[i(\vec{k}\cdot\vec{r}-\omega t)]$		(8)

The different Fourier components are not coupled, and are treated one at a time, as if $|X(\vec{r},t)\rangle=|X(\vec{k},\omega)\rangle\exp[i(\vec{k}\cdot\vec{r}-\omega t)]$. The symbol $\omega$ is the external frequency, not to be confused with $\hat{\omega}$ which is the operator version of the phonon frequency $\omega_{Q}$. Equations 6 and 7 become

$$\left[\hat{\Omega}+i(\vec{k}\cdot\hat{\vec{v}}-\omega\hat{1})\right]|\Phi\rangle=|A(\vec{k},\omega)\rangle$$

(9)

$$|A\rangle=\frac{\sqrt{\hat{n}(\hat{n}+1)}}{k_{B}T^{2}}\left[-i(\vec{k}\cdot\hat{\vec{v}}-\omega\hat{1})\Delta T|\hbar\omega\rangle+\frac{\hbar\hat{\omega}}{C}|P\rangle\right]$$

(10)

Mode space has so far been described in $Q$-representation by harmonic eigenstates $|Q\rangle$. It is convenient to also use the “relaxon”-representation $|\alpha\rangle$ of eigenstates Guyer and Krumhansl (1966); Maris (1969); Cepellotti and Marzari (2016) of $\hat{\Omega}$.

$$\hat{\Omega}|\alpha\rangle=\gamma_{\alpha}|\alpha\rangle$$

(11)

where

$$\langle\alpha|\beta\rangle=\delta_{\alpha\beta}\ \ {\rm and}\ \ \sum_{\alpha}|\alpha\rangle\langle\alpha|=\hat{1}.$$

(12)

The eigenvalues $\gamma_{\alpha}$ are relaxation rates, $\gamma_{\alpha}=1/\tau_{\alpha}$. In this basis, Eq. 9 has the form

$$\sum_{\beta}\left[\gamma_{\alpha}\delta_{\alpha\beta}+i(\vec{k}\cdot\vec{v}_{\alpha\beta}-\omega\delta_{\alpha\beta})\right]\Phi_{\beta}=A_{\alpha},$$

(13)

where $\Phi_{\beta}=\langle\beta|\Phi\rangle$ and $A_{\alpha}=\langle\alpha|A\rangle$. There are $N$ modes $Q$ ($N=3N_{\rm cells}n_{\rm at}$ where $n_{\rm at}$ is the number of atoms in the unit cell), and Eq. 13 gives $N$ equations for the $N$ unknown components $\Phi_{\alpha}$ of the deviation function. But there are two fields ($|\Delta T\rangle$ and $|P\rangle$) driving the distribution out of equilibrium, of which one (typically $|\Delta T\rangle$) is unknown. An extra equation is needed. That equation is the definition of local temperature.

There is one null eigenvector, $|\alpha=0\rangle\equiv|0\rangle$ with eigenvalue $\gamma_{0}=0$. The vectors $\langle 0|$ and $|0\rangle$ deviate only by factors $[\hat{n}(\hat{n}+1)]^{\pm 1/2}$ from the null left $(\langle\hbar\omega|)$ and right $|dn/dT\rangle$ eigenvectors of $\hat{R}$. In $Q$-representation, the null eigenvector is

$$\langle Q|0\rangle=\sqrt{n_{Q}(n_{Q}+1)}\frac{\hbar\omega_{Q}}{\sqrt{CVk_{B}T^{2}}}$$

(14)

The factor $1/\sqrt{VCk_{B}T^{2}}$ normalizes the state, $\langle 0|0\rangle=1$. The definition of temperature takes the form

	$\displaystyle\frac{1}{V}\sum_{Q}$		$\displaystyle\hbar\omega_{Q}\Delta N_{Q}=0=\frac{1}{V}\sum_{Q}\hbar\omega_{Q}\sqrt{n_{Q}(n_{Q}+1)}\Phi_{Q}$		(15)
			$\displaystyle=\sqrt{\frac{Ck_{B}T^{2}}{V}}\langle 0|\Phi\rangle=\sqrt{\frac{Ck_{B}T^{2}}{V}}\Phi_{0}=0.$

This shows that there are actually only $N-1$ unknown parts of $|\Phi\rangle$, because the $\alpha=0$ component, $\langle 0|\Phi\rangle=\Phi_{0}$ must be zero by the definition of local temperature.

Now we can rewrite the formula for $|A\rangle$ using Eqs. 10 and 14,

$$|A\rangle=\sqrt{\frac{V}{Ck_{B}T^{2}}}\left[-i(\vec{k}\cdot\hat{\vec{v}}-\omega\hat{1})C\Delta T(\vec{k},\omega)+\hat{P}(\vec{k},\omega)\right]|0\rangle$$

(16)

where $\langle Q|\hat{P}|Q^{\prime}\rangle=P_{Q}\delta_{QQ^{\prime}}$. Now look at the $\alpha=0$ component of Eq. 13,

$$i\vec{k}\cdot\sum_{\beta}^{\neq 0}\vec{v}_{0\beta}\Phi_{\beta}=A_{0}$$

(17)

The left hand side of Eq. 17 is

$$\langle 0|i\vec{k}\cdot\vec{v}|\Phi\rangle=\sqrt{\frac{V}{Ck_{B}T^{2}}}i\vec{k}\cdot\vec{J}$$

(18)

where $\vec{J}(\vec{k},\omega)$ is the energy (or heat) current density, $(1/V)\sum_{Q}\hbar\omega_{Q}\vec{v}_{Q}[N_{Q}(\vec{k},\omega)-n_{Q}(\vec{k},\omega)]$. The right hand side of Eq. 17 is

$$\langle 0|A\rangle=\sqrt{\frac{V}{Ck_{B}T^{2}}}\left(i\omega C\Delta T+\bar{P}\right),$$

(19)

where $\bar{P}=\langle 0|\hat{P}|0\rangle$, or

$$\bar{P}\equiv\sum_{Q}\frac{C_{Q}}{C}P_{Q}\ \ {\rm and}\ \ C_{Q}=\frac{1}{V}\hbar\omega_{Q}\frac{dn_{Q}}{dT},$$

(20)

and $C=\sum_{Q}C_{Q}$ is the total specific heat. If $P_{Q}$ is independent of $Q$, then $P_{Q}=\bar{P}$ for all modes $Q$. Thus the $\alpha=0$ part of the Boltzmann equation expresses energy conservation. In $(\vec{r},t)$ representation,

$$\vec{\nabla}_{\vec{r}}\cdot\vec{J}=-\frac{\partial u}{\partial t}+\bar{P}.$$

(21)

where $\partial u/\partial t=C\partial\Delta T/\partial t\rightarrow-i\omega C\Delta T$.

We can now “reinterpret” Eq. 13 as a system of $N-1$ linear equations (the sum over $\beta$ does not include $\beta=0$ because $\Phi_{0}=0$). The $(N-1)\times(N-1)$ matrix $\Omega_{\alpha\beta}=\gamma_{\alpha}\delta_{\alpha\beta}$ is now positive-definite. We still need to invert the non-Hermitean matrix on the left side of Eq. 13 in order to solve for $|\Phi\rangle$.

The method of solution is given by Hua and Lindsay Hua and Lindsay (2020). A earlier version is in a paper by Cepellotti and Marzari Cepellotti and Marzari (2017). Rewrite Eq. 13 by rescaling the distribution function $\Phi_{\alpha}$ and driving term $A_{\alpha}$:

$$\Psi_{\beta}\equiv\gamma_{\beta}^{1/2}\Phi_{\beta}\ {\rm and}\ B_{\alpha}\equiv\gamma_{\alpha}^{-1/2}A_{\alpha}.$$

(22)

Equation 13 then becomes $[\hat{1}+i\hat{\Gamma}]|\Psi\rangle=|B\rangle$, or

			$\displaystyle\sum_{\beta}\left[\delta_{\alpha\beta}+i\Gamma_{\alpha\beta}\right]\Psi_{\beta}=B_{\alpha}$
	$\displaystyle\Gamma_{\alpha\beta}$	$\displaystyle=$	$\displaystyle\vec{k}\cdot\left(\gamma_{\alpha}^{-1/2}\vec{v}_{\alpha\beta}\gamma_{\beta}^{-1/2}\right)-\gamma_{\alpha}^{-1}\omega\delta_{\alpha\beta}$		(23)

The matrix $\hat{\Gamma}$ is real-symmetric, so it has real eigenvalues, $\lambda_{m}$:

$$\hat{\Gamma}(\vec{k},\omega)|m\rangle=\lambda_{m}(\vec{k},\omega)|m\rangle.$$

(24)

In this basis, Eq. 23 is $(1+i\lambda_{m})\Psi_{m}=B_{m}$, where $\Psi_{m}=\langle m|\Psi\rangle$, etc. Then the distribution function in the relaxon basis (the eigenbasis of $\hat{\Omega}$) is

	$\displaystyle\Psi_{\alpha}=\sum_{\beta}S_{\alpha\beta}B_{\beta},\ {\rm or}\ \ \Phi_{\alpha}=\sum_{\beta}\left(\gamma_{\alpha}^{-1/2}S_{\alpha\beta}\gamma_{\beta}^{-1/2}\right)A_{\beta},$
	$\displaystyle{\rm where}\ \ \langle\alpha|\hat{S}|\beta\rangle=S_{\alpha\beta}=\sum_{m}\langle\alpha|m\rangle\frac{1}{1+i\lambda_{m}(\vec{k},\omega)}\langle m|\beta\rangle.$
			(25)

This is the desired solution. In the spatially homogeneous ($\vec{k}=0$) and static ($\omega=0$) case, $\hat{\Gamma}=0$ and $\lambda_{m}=0$, so $S_{\alpha\beta}=\delta_{\alpha\beta}$ and the bulk solution $\Phi_{\alpha}=\gamma_{\alpha}^{-1}A_{\alpha}$ is recovered.

The operator $\hat{\Omega}$ is positive if we exclude the null space. The operator $\hat{S}=(\hat{1}+i\hat{\Gamma})^{-1}$ is defined in this “positive” or $p$-space. It is convenient to define another operator in the same $p$-space,

	$\displaystyle\hat{W}$	$\displaystyle\equiv$	$\displaystyle\hat{\Omega}^{-1/2}\hat{S}\hat{\Omega}^{-1/2}$
	$\displaystyle W_{\alpha\beta}$	$\displaystyle=$	$\displaystyle\gamma_{\alpha}^{-1/2}S_{\alpha\beta}\gamma_{\beta}^{-1/2}$		(26)

The solution of the Boltzmann equation is then

$$|\Phi\rangle=\hat{W}|A\rangle_{p}$$

(27)

where $|A\rangle_{p}$ is the part of the driving term that is orthogonal to $|0\rangle$ (and thus lies in $p$-space),

$$|A\rangle_{p}=\sqrt{\frac{V}{Ck_{B}T^{2}}}\left[-i\vec{k}\cdot\hat{\vec{v}}C\Delta T+(\hat{P}-\bar{P}\hat{1})\right]|0\rangle.$$

(28)

As required, the inner product $\langle 0|A\rangle_{p}$ vanishes, because it has a term proportional to $\langle 0|\hat{\vec{v}}|0\rangle$ which is zero because $\hat{\vec{v}}$ is an odd operator, and a term proportional to $\langle 0|(\hat{P}-\bar{P}\hat{1})|0\rangle$ that is zero because $\langle 0|\hat{P}|0\rangle\equiv\bar{P}$. From $|\Phi\rangle$ we get the heat current density $\vec{J}$,

$$\vec{J}=\sqrt{\frac{Ck_{B}T^{2}}{V}}\langle 0|\hat{\vec{v}}|\Phi\rangle=\vec{J}_{\rm F}+\vec{J}_{\rm non-F}.$$

(29)

The “generalized Fourier” component $\vec{J}_{\rm F}$ is

$$\vec{J}_{\rm F}(\vec{k},\omega)=-\kappa(\vec{k},\omega)\cdot(i\vec{k})\Delta T(\vec{k},\omega),$$

(30)

$$\kappa(\vec{k},\omega)=C\ \langle 0|\hat{\vec{v}}\ \hat{W}(\vec{k},\omega)\ \hat{\vec{v}}|0\rangle.$$

(31)

The “non-Fourier” term is

$$\vec{J}_{\rm non-F}(\vec{k},\omega)=\langle 0|\vec{\hat{v}}\hat{W}[\hat{P}(\vec{k},\omega)-\bar{P}(\vec{k},\omega)\hat{1}]|0\rangle.$$

(32)

This extra, non-Fourier, component of the current was found by Hua et al. Hua et al. (2019) in an RTA treatment. A version is in the paper by Mahan and Claro Mahan and Claro (1988). In their version, there was no external insertion $P$ except via boundary conditions. It was rederived by Hua and Lindsay Hua and Lindsay (2020) in a more complete treatment (not using RTA). Reference Hua et al., 2019 gives additional parts of the non-Fourier term that arise from boundary conditions, but solved only in RTA. Boundary terms do not appear here; the present derivation, like ref. Hua and Lindsay, 2020, assumes an infinite homogeneous sample.

In a d.c. situation ($\omega=0$), the insertion term $\hat{P}(\vec{r})$ should have a zero spatial average: $P(\vec{k}=0,\omega=0)=0$. Otherwise, the sample experiences net heating or cooling. The spatially homogeneous ($\vec{k}=0$) part of the “non-Fourier” current is zero. Also, when $\vec{k}=0$ and $\omega=0$, $\hat{W}(0,0)=\hat{\Omega}^{-1}$ and the bulk static thermal conductivity is recovered,

$$\kappa(\vec{k}=0,\omega=0)=C\ \langle 0|\hat{\vec{v}}\ \hat{\Omega}^{-1}\ \hat{\vec{v}}|0\rangle.$$

(33)

The local temperature distribution $T(\vec{r},t)$ can be found from energy conservation (Eq. 21),

$$\Delta T(\vec{k},\omega)=\frac{1}{i\omega C}\left[i\vec{k}\cdot\left(\vec{J}_{\rm F}(\vec{k},\omega)+\vec{J}_{\rm non-F}(\vec{k},\omega)\right)-\bar{P}(\vec{k},\omega)\right].$$

(34)

If $P_{Q}(\vec{k},\omega)$ is independent of $Q$ (i.e. $\hat{P}=\bar{P}$), the non-Fourier current is zero, and the local temperature variation is given by a simple nonlocal response function, the “thermal distributor” $\Theta(\vec{k},\omega)$,

$$\Delta T(\vec{k},\omega)=\Theta(\vec{k},\omega)\bar{P}(\vec{k},\omega)$$

(35)

where

$$\Theta(\vec{k},\omega)=\frac{1}{\vec{k}\cdot\kappa(\vec{k},\omega)\cdot\vec{k}-i\omega C}$$

(36)

This is Eq. 38 of ref. Hua and Lindsay, 2020. The function $\Theta$ was defined in ref. Allen and Perebeinos, 2018. The name “thermal distributor,” has been changed from the previous name “thermal susceptibility,” and a factor of C has been removed from the definition. When $\hat{P}\neq\bar{P}$, there is another term,

$$\Delta T(\vec{k},\omega)=\Theta(\vec{k},\omega)\bar{P}(\vec{k},\omega)-\frac{i\vec{k}\cdot\vec{J}_{\rm non-F}(\vec{k},\omega)}{\vec{k}\cdot\kappa(\vec{k},\omega)\cdot\vec{k}-i\omega C}$$

(37)

Griffin Griffin (1965) makes an interesting argument that may be related to thermal distributors. Guyer and Krumhansl Guyer and Krumhansl (1966) noticed that Griffin’s argument is related to the mean-field method relating non-interacting to interacting electrical susceptibility.

Various names are given to the operator $-\vec{v}_{Q}\cdot\vec{\nabla}_{\vec{r}}$: (a) “drift operator”, (b) “advection operator”, (c) “convection operator”, or (d) “diffusion operator”. Sometimes the combination $\partial/\partial t+\vec{v}_{Q}\cdot\vec{\nabla}_{\vec{r}}$ is called the “drifting operator” . The process these words are describing is that the occupancy $N_{Q}(\vec{r},t+\Delta t)$, in a completely ballistic system, is the same as $N_{Q}(\vec{r}-\vec{v}_{Q}\Delta t,t)$, so that

	$\displaystyle\left(\frac{dN_{Q}}{dt}\right)_{\rm drift}$	$\displaystyle=$	$\displaystyle\lim_{\Delta t\rightarrow 0}\left[\frac{N_{Q}(\vec{r},t+\Delta t)-N_{Q}(\vec{r},t)}{\Delta t}\right]_{\rm ballistic}$
		$\displaystyle=$	$\displaystyle-\vec{v}_{Q}\cdot\vec{\nabla}_{\vec{r}}N_{Q}.$

The term I use, “drift operator”, makes sense; “diffusion operator” does not.