Heat vortex in hydrodynamic phonon transport of two-dimensional materials

We study hydrodynamic phonon heat transport in two-dimensional (2D) materials. Starting from the Peierls-Boltzmann equation with the Callaway model approximation, we derive a 2D Guyer-Krumhansl-like equation describing hydrodynamic phonon transport, taking into account the quadratic dispersion of flexural phonons. In addition to Poiseuille flow, second sound propagation, the equation predicts heat current vortices and negative non-local thermal conductance in 2D materials, which are common in classical fluids but have not yet been considered in phonon transport. Our results also illustrate the universal transport behaviors of hydrodynamics, independent of the type of quasi-particles and their microscopic interactions.


The steady state distribution functions
The steady state distribution functions are obtained by considering the conserved quantities during the scattering process. For N-process, both crystal momentum and energy are conserved. If we take a 3-phonon scattering process as an example, we have ω k k k 1 1 1 + ω k k k 2 2 2 = ω k k k 3 3 3 , (S1) k k k 1 1 1 + k k k 2 2 2 = k k k 3 3 3 .
On the other hand, when reaching the steady state, the forward and backward scattering processes are balanced. Thus, the phonon distribution function satisfies f k k k 1 1 1 f k k k 2 2 2 (1 + f k k k 3 3 3 ) = (1 + f k k k 1 1 1 )(1 + f k k k 2 2 2 ) f k k k 3 3 3 . ).
Comparing Eqs. (S1-S4), we find that we can write the distribution function in terms of the conserved quantities ln(1 + f −1 k k k ) = β (hω −hk k k · u u u).
This results in a displaced Bose-Einstein distribution function Similarly, we can show that the U-process drives the phonons to the equilibrium Bose-Einstein distribution.

Derivation of the 2D G-K equation
Here, we give the derivation of 2D G-K equation in the main text. We use the recent developed multiscale expansion technique ? . The expansion is over both space and time as follows: where ε is a small parameter, defined as This means that we consider the situation where the scattering rates of N-process are much larger than that of R-process. This is the necessary condition for hydrodynamic phonon transport. The phonon distribution function f can be expanded similarly with the 0th, 1st and 2nd order distribution function f 0 , f 1 , f 2 . The macroscopic variables can be expressed by the sum of corresponding components. Substituting all of the equations above into the Peierls-Boltzmann equation, we obtain the explicit forms of f 0 and f 1 : (S11) and different order of approximations for the energy and heat-flux balance equations To derive the hydrodynamic equations, we need to write the macroscopic variables of different orders using f 0 and f 1 .

The zeroth-order result
The zeroth-order results are obtained by making the following approximation This approximation is valid if |hk k k · u u u| hω. It has been checked numerically ? .

Energy density
We consider the phonon energy density first. According to energy conservation, we know that E = E 0 . (S16) The second term in the above equation vanishes because it's an odd function of k k k. Thus we arrive at: Writing it as summation of the linear E L and quadratic E N contributions, we get with (S19) Note that we have included the factor 2 to account the degeneracy of the linear phonon modes. From the energy above, we obtain the specific heat capacity as We also define an effective average energy and specific heat capacity as whereĒ will be used in the expression for the heat flux.

Momentum density
Analogously, due to the momentum conservation of N-process, we have where the linear and quadratic terms take the following forms As discussed in the main text, to the first order in u u u, p p p N diverges logarithmically. While for small but finite u u u, p p p N is finite. Thus, we need to go beyond linear order in u u u when calculating p p p N .

Derivation of Eq. (S24)
We start from For simplicity, we set the drift velocity u u u to be parallel to the x axis, so and (p N ) y = 0.
Putting together, we get ) u u u.

Heat flux
The 0th order heat flux is calculated as Since at thermal equilibrium, the heat flux is zero. Only the second term contributes where q q q 0L , q q q 0N are the 0th order contributions of the linear and quadratic modes to the heat flux. Using Eq. (S21), we have q q q 0 =Ēu u u.

Thermal conductivity
The thermal conductivity tensor is Likewise, the second term in the above equation is an odd function of k k k and thus is zero. We then have It is noted that, given a constant τ R , the thermal conductivity contributed from a quadratic mode is two times larger than that contributed from a linear mode. Both of them are independent of the details of the dispersion.

The first-order result
Now, we calculate the 1st order terms using the expansion Based on the chain rule, we have To proceed, we have two options for the derivation of G-K equation. These two options are exactly the same for phonons with linear dispersion, i.e. the Debye model. But we must pay attention to the quadratic phonon mode here, it leads to different results. The underline physics is that, for phonons with linear dispersion, the momentum density and the heat flux are proportional to each other q q q L = v 2 g p p p L .
The group velocity v g is constant for phonons with linear dispersion. But for quadratic dispersion, this simple relationship does not hold, due to the wavevector dependence of v g , i.e. v g ∝ k.
Here, to consider heat transport, we start from Eqs. (S13-S14) and arrive at a lengthy expression for f 1 Note that, in the expression for f 1 we have one term ∝ q 1i . On the other hand, to get q 1i we need to know f 1 . This means the equations are coupled together and have to be solved self-consistently. Fortunately, the q 1i term in Eq. (S35) is odd in k and does not contribute to ← → κ κ κ 1 and the viscosity coefficients in the final G-K equation. We get for ← → κ κ κ 1

The heat G-K equation
Substituting the calculated 0th and 1st order results into neglecting the nonlinear effects resulting from the product of heat flux and temperature gradient in κ 1 and keeping only the lowest order term in the heat flux q q q ≈ q q q 0 , we get the G-K equation for heat flux The transport coefficients are κ 0 = αCv 2 g τ R , η = β v 2 g τ N , ζ = αv 2 g τ N , with α = C L /C, β = 9E L /8Ē, respectively. All the three transport coefficients can be divided into contributions from the linear and quadratic phonon modes, respectively. Denoting them with subscripts L and N, It can be checked that, our results reduce to that of Debye model if we ignore the quadratic phonon mode.

Steady state heat flow in a ribbon
We proceed to analyse the hydrodynamic heat transport in a graphene nano-ribbon. At steady state, we have ∂ q q q/∂t = ∂ E/∂t = 0 and Eq. (S42) reduce to η∇ 2 q q q − τ −1 R q q q = C L v 2 g ∇T.
Similar to Ref. ? , we introduce a stream function ψ with q q q = z z z × ∇ψ. A flow of heat current from a point source I(x) = Iδ (x) is injected into the ribbon and collected at the opposite side. So the normal heat flux is q y (x, y) y=0,w = I(x).