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.


Methods
We consider a prototype 2D material. It has one out-of-plane acoustic mode with a quadratic dispersion ak k 2 ω = (ZA mode), and two degenerate linear acoustic modes ω = v k k g (longitudinal and transverse). The magnitude of the linear group velocity is v g , and the magnitude of the wave vector is k = |k|. Here, in the spirit of Debye model, we ignore the possible anisotropic property and the difference between longitudinal and transverse branches. We will focus on the effects of ZA mode with quadratic dispersion on the hydrodynamic behaviors.
We first sketch the derivation of the 2D G-K equation 45,46 . Our starting point is the Peierls-Boltzmann equation under Callaway approximation 47,48 τ τ Here, s is the phonon index, τ N is the constant relaxation time for the N-process, while τ R is that for the resistive scattering process (R-process). It includes all scattering mechanisms that do not conserve crystal momentum, i. e., impurity scattering, electron scattering and other U-processes. We take the same τ R and τ N for all phonon branches, i.e., taking the wave vector and branch averaged values. This is the so-called gray approximation. For quantitative analysis of specific materials, one needs to go beyond the gray approximation and consider the wave vector and branch dependence of τ, which is beyond the scope of present study. The phonons may reach different steady state distributions due to different kinds of scattering processes, depending on the conserved quantities before and after the scattering. For the N-process, both energy and crystal momentum are conserved. It drives the system towards a displaced Bose-Einstein distribution function It has been shown numerically that, within some moderate temperature range (~100 K for graphene), the N-process is orders of magnitude faster than the R-process, meaning R N  τ τ 10,11 . When the system size is much larger than the normal-scattering mean free path l v g N τ ≈ . The relaxation from local to global equilibrium or steady state is governed by hydrodynamic equations describing the conserved quantities during the N-process, including energy and crystal momentum.
As we know, conservation laws play central roles in the derivation of hydrodynamic equations. Since phonons can be generated or destroyed during the scattering processes, the number of phonons is not conserved. We consider energy and crystal momentum conservation here. Both normal (N) and resistive (R) scattering processes conserve total energy, giving

Energy balance equation.
To get the equation governing the dynamics of the conserved quantity, we multiply by ω  sk , integrate over k and sum over the phonon index on both sides of Eq. (1). We then arrive at an equation describing energy conservation where is the energy density, and is the heat flux. The right hand side of Eq. (7) is zero because the scattering processes conserve energy.

Momentum balance equation. Multiplying k
 and performing the integration/summation, we can obtain an equation for the crystal momentum density p and its flux Φ ↔ from its conservation law during N-process R with the momentum density and momentum flux tensor Here, the presence of last term at the right side of the equation is due to the fact that normal scattering process conserves crystal momentum, but not the heat flux Here, q 0 is defined similarly by replacing f sk by f k N s eq , . We have defined the thermal conductivity tensor Multi-scale expansion. To obtain the hydrodynamic heat transport equation, we follow a multi-scale expansion technique 9 and extend it to the case of 2D material with one quadratic phonon dispersion. This quadratic phonon dispersion is important in 2D materials like graphene 11 . The expansion is over both space and time as follows: It is a perturbation expansion over a natural small parameter N R This means 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 www.nature.com/scientificreports www.nature.com/scientificreports/  Other quantities are defined similarly. According to energy and momentum conservation (Eqs. (4) and (6)), we know that E = E 0 , p = p 0 . If all the phonon modes follow linear dispersion, q is simplied proportional to p (Eq. (30)), thus q = q 0 . This is the case for the Debye model. However, for 2D materials, when we include the quadratic phonon mode, the simple relation does not hold any more, q ≠ q 0 .
All the defined quantities are calculated from the distribution function ( ) In most cases, we can approximate f N eq as Taking into account only the 0th order term f 0 , we can get Here, Z(x) is the Riemann Zeta function, the index L and N represent the linear and quadratic phonon contributions, respectively. The calculation of p 0N needs special care. Simple calculation using Eq. (24) leads to divergent result. We need to use the full form of f N eq instead of Eq. (24) and obtain We see that it depends on u non-linearly, and the linear coefficient diverges as → u 0. As a result, we have to consider a finite u when calculating p N . This makes the linear approximation in u and consequently the derivation of G-K equation using the momentum balance equation difficult. The other complication that the quadratic mode brings is the lack of simple proportionality between p N and q N . For linear dispersion, we have = q p v L g L 0 2 0 . As a result, when only considering phonons with linear dispersion, the same G-K equation can be derived starting from either the momentum balance or the heat-flux equation. The inclusion of the quadratic dispersion leads to two different G-K-like equations, describing momentum and heat flow, respectively. Here, we focus on the heat transport, and proceed from the heat-flux equation (13). We get the first order correction to the thermal conductivity The expression of ↔ Q and the details of the derivation can be found in the Supplementary Information (SI).
The 2D G-K equation. Substituting the above results into Eq. (13), and neglecting q 1 , we arrive at the G-K-like equation for heat transport Here, considering the isotropic case, the transport coefficients can be expressed as numbers. κ 0 is the magnitude of the zeroth order thermal conductivity, η and ζ are the first, second viscosity coefficients, 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, we have It can be checked that, our results reduce to that of Debye model if we ignore the quadratic phonon mode. We can also write κ α τ = Cv g R . Equation (33) with the above defined coefficients is the central result of this  , β = 1/5 for 3D 9 , respectively. We note by passing that, in deriving Eq. (33) we have made the approximation q q 0 ≈ . Inclusion of q 1 requires solution of higher order equations in the expansion over ε. Analytical treatment becomes difficult, if not impossible. Thus, we rely on fully numerical solution of the Callaway model to check its validity, as we have done in Fig. 3. We show that our G-K equation can reasonably re-produce the main features of the hydrodynamic heat flow that we focus in this work. Physically, this approximation should hold when 1 ε  or  τ τ N R , which is the parameter range we consider here. One additional support for our approximation is that the viscosity coefficients, which come from κ ↔ 1 , do not depend on q 1 (see Sec. 2.2 of SI).

Results
Second sound. The right side (RHS) of the G-K equation represents the effect of viscosity on the heat transport behavior. They come from the first order term in the expansion over ε. Before looking into these terms, we show here that the propagation of second sound can be analyzed without these terms. Replacing the RHS with zero, combining with Eq. (7), we arrive at where we have defined the second sound velocity This is the wave equation describing propagation of second sound with velocity v ss and damping coefficient Here, the presence of quadratic dispersion makes v ss temperature dependent, inherited from the different temperature dependence of C L and C N [see Fig. 2(a)]. We find that, in the presence of quadratic ZA modes, v ss increases with temperature and is much smaller than results from the Debye model in the relevant temperature ~100 K. The physical mechanism is the following. The quadratic dispersion of ZA modes have τ = − s, 0 5 χ = . the same as that in Fig. 1(c). This is compared to Fig. 1(c). The differences of heat current streamlines near the left and right boundaries are due to the periodic boundary condition used in the numerical calculation. www.nature.com/scientificreports www.nature.com/scientificreports/ frequency dependent phonon group velocity, i.e., ∝ v k g , and k-independent constant density of states. This is contrast to constant v g and linear-in-k density of states of phonons with linear dispersion (LA and TA modes). At low temperature, phonons with small k contribute dominantly to the propagation of second sound. Since ZA modes contribute dominantly, v ss is much smaller than the constant v g of linear phonons. When the temperature increases, more linear phonons contribute and v ss increases correspondingly.
Poiseuille flow. We now include the RHS of the G-K equation, and consider a nano-ribbon with length L ( x L 0 ≤ ≤ ) and width w ( y w 0 ≤ ≤ ) [ Fig. 1(a)]. A temperature difference is applied along the ribbon (x direction) [ Fig. 1(b)]. At steady state, ignoring τ q/ R , Eq. (33) reduces to a one-dimension form q y A / 2 2 ∂ ∂ = . This gives rise to a parabolic heat flux distribution perpendicular to the flow q y Ay y w ( ) , if we assume a non-slip boundary condition = = q qw (0) ( ) 0. By integration over y, the heat current is obtained The negative sign means heat flows opposite to the temperature gradient. The heat current scaling as w 3 is a signature of the Poiseuille flow. For diffusive phonon transport, the heat current scales linearly with the ribbon width I w ∝ , while for ballistic transport, the heat current can not go higher than linear scaling with the width. Thus, the cubic (super-linear) dependence of I on w can in principle be used as a signature of the Poiseuille flow. The Poiseuille flow in graphene ribbons has been studied numerically by solving the Boltzmann equation directly in refs. 10,11 . Similar behavior is also predicted for electronic transport in graphene nano-ribbons 12 . Length and width dependent thermal conductivity in suspended single layer graphene has been reported experimentally 49,50 . Thus, experimental confirmation of Poiseuille flow in graphene is already within reach.

Steady state heat flow in a ribbon. One important consequence of Eq. (33) is the formation of heat flow
vortices when there is a heat current source injecting into the 2D materials. As far as we know, this has not been considered before, on which we focus in this work. We consider steady state transport in a setup sketched in Fig. 1(c). Heat current source and drain are attached to a graphene nano-ribbon. The pattern of heat current flow at steady state can be obtained from the solution of simplified version of Eq. (33). At steady state, according to Eq.
The resulting equation has the form It shares the same form as the electronic case in ref. 12 . The 1st term on the left hand side (LHS) is the hydrodynamic term due to viscosity. When R 1 τ − is negligible, we get a pure hydrodynamic viscous flow. On the other hand, when η is negligible, we recover the normal diffusive heat transport governed by Fourier law. Actually, it is suggestive to define a dimensionless parameter to characterize the relative contribution of diffusive and hydrodynamic transport. χ basically characterizes the relative contribution of the first and the second term in Eq. (37). In the limit of χ → + ∞, the R-process is dominant, Eq. (37) reduces to the diffusive Fourier law. However, in the other limit 0 χ → , the viscosity term dominates.
Following ref. 12 , we have solved Eq. (37) analytically with the help of the streaming function (see Sec. 3 of the SI for details). As an example, we have plotted typical heat current flow patterns (lines) and the resulting temperature distribution (color) with 0 5 χ = . and 2 × 10 4 in Fig. 1(c,d), respectively. Here, a flow of heat current from a point source δ = I x I x ( ) ( ) is injected into the ribbon and collected at the opposite side. The non-slip boundary condition is used to solve Eq. (37). The heat current flow within the ribbon can be obtained from the solutions. For small χ or larger η [ Fig. 1(c)], hydrodynamic transport is dominant. The formation of vortices at both sides of the direct source-to-drain flow is a characteristic feature of the viscous flow. This feature is shared by different kinds of classical or quantum fluid. Similar behavior of electrons in graphene and other 2D materials has received intense research focus very recently [4][5][6]9,22,51 . As a results of vortices formation, there appears the separation of temperature gradient and heat flow. Even negative thermal resistance can be observed, where the heat current flows are from the low to the high temperature regime. This is an obvious violation of the Fourier law. For larger χ [ Fig. 1(d)], the system is in the diffusive Fourier transport regime. Heat current vortices and negative thermal resistance are absent. Thus, we can realize a transition from hydrodynamic to Fourier transport by changing the magnitude of χ.
Here, it is worth mentioning that, for much smaller nanoscale systems, the wave property of phonons and complicated elastic boundary scattering may also leads to the formation of vortices in the ballistic transport regime 34 . Despite the similarity, the physical mechanism and length scale are quite different from the viscous flow studied here. In the ballistic case, coherent phonon transport together with elastic boundary scattering is the physical mechanism to generate heat vortices. Here in the hydrodynamic case, vortex formation is due to frequent momentum-conserving N-processes together with the boundary conditions imposed here, i.e., heat current injection and collection at local positions. It is a result of frequent momentum exchange between different phonons, akin to the vortex formation in classical gas or liquid flow. While the ballistic phonon transport takes Scientific RepoRtS | (2020) 10:8272 | https://doi.org/10.1038/s41598-020-65221-8 www.nature.com/scientificreports www.nature.com/scientificreports/ place in nanoscale ribbons, the hydrodynamic transport takes place in the microscale (See below for the estimation of length scale in graphene).
As we mentioned above, in deriving Eq. (33), we have considered only the zeroth order term of the heat current q. To check the validity of this truncation and the results plotted in Fig. 1(c), we have performed fully numerical calculation by solving the Boltzmann equation under Callaway approximation using the discrete ordinate method 33,52 . More details of the numerical calculation can be found in refs. 33,52 . We consider the same ribbon at the same average temperature (T 100 0 = K) as in Fig. 1(c), with the same dimensionless parameter χ. The temperature difference in the numerical calculation is chosen such that the heat flux q flowing into the ribbon is the same as that in Fig. 1(c). In this way, we can compare directly the results obtained from the two methods. The numerical calculation serves as a benchmark for our truncation q 1 = 0. The results are plotted in Fig. 3(a). Comparing with Fig. 1(c), we can find that the main features of hydrodynamic heat flow are re-produced by the G-K heat equation derived here. Figure 3(b) shows the distribution of relative difference between q and q 0 both obtained from the numerical calculation. We have checked that, the maximum relative difference, defined as |(q − q 0 )/q|, locates at the upper and lower boundaries. Part of the difference comes from the singular boundary conditions, i.e., the heat current injected and collected are distributed as delta changes. The relative difference is rather small away from the two boundaries, especially around the vortices. As examples, we have plotted the distribution of heat flux along two line cuts (green dashed lines in Fig. 3) passing the center of the ribbon in y (Fig. 3(c)) and x (Fig. 3(d)) directions, respectively. The analytical and numerical results show good agreements in both cases. This further validates our study of the hydrodynamic heat flow using the G-K equation derived here.
Application to graphene. We now give an order-of-magnitude analysis using parameters of graphene. We obtain the phonon dispersion relation of graphene using density function theory based calculations (For the density functional theory calculation, we use the Vienna Ab-initio Simulation Package and the generalized gradient approximation for the exchange-correlation functional. The parameters are the same as ref. 53  ≈ . × π −  Jm −2 K −1 , respectively. To estimate the transport coefficients and the dimensionless factor χ, we use 10 N 10 τ = − s, τ = − 10 R 7 s 11 . We get 0 08 β ≈ . at = T 100 K, smaller than value obtained from the 2D Debye model 1/4 β = . In contrast to 2D or 3D Debye model, the different temperature dependences of C L , C N and E L , E N give rise to temperature dependent η (Fig. 2(b)). We get η ≈ .
0 002 m 2 /s at = T 100 K, which is orders of magnitude larger than that of water. We also get the dimensionless parameter w 5 10 9 2 χ ≈ × m −2 . Thus, phonon hydrodynamic transport can be realized in high quality graphene nano-ribbons of micrometer scale. The plot in Fig. 1(c) with 0 5 χ = . corresponds to sample size of ~10 μm.

conclusions
In summary, we have derived a 2D version of the G-K equation describing hydrodynamic phonon heat transport. We take into account the out of plane quadratic phonon dispersion of the ZA mode, normally present in 2D materials. Its effect on the hydrodynamic transport is analyzed. The derived equation serves as a starting point for investigating hydrodynamic phonon transport behavior in 2D materials. It shares a similar form as the Navier-Stokes equation that has been used to study electron hydrodynamic transport in graphene. Many interesting transport behaviors, including non-local negative resistance, higher-than-ballistic transport, predicted for electrons can be studied for phonons. Moreover, a large overlap of the parameter regime between electron and phonon hydrodynamic transport in graphene makes it promising to study the effect of their mutual interaction on the thermal transport behavior of the two kinds of fluids.