## Abstract

Heat flow along two-dimensional strips as a function of the Knudsen number is examined in two different versions of heat-transport equations with non-local terms, with or without heat slip flow. In both of them, a parabolic heat profile corresponding to the Poiseuille phonon flow may appear in some domains of temperature, or of the Knudsen number, in the transition from the Fourier regime to the ballistic one. The influence of the slip heat flow on such a transition is discussed.

## 1. Introduction

Material systems with characteristic sizes of the order of nanometres exhibit interesting optical, magnetic, electrical and/or photoelectrical properties. In recent years, several research groups have focused their attention on the contribution of collective hydrodynamic-like phonon behaviour on two-dimensional nanosystems [1–4]. Silicon nanolayers, boron nitride nanolayers [5–7] and graphene [8–10] have especially attracted attention. In particular, it has been recently argued that, in suspended graphene [4] and other two-dimensional systems [3], Poiseuille phonon flow could be observed at higher temperatures and in wider temperature ranges than in three-dimensional systems, because in two-dimensional systems the normal (momentum conserving) phonon–phonon collisions are two orders of magnitude higher than those in three-dimensional systems. Thus, analysis of this regime in two-dimensional situations has become a topic of current interest.

Phonons, which are the main heat carriers in non-metallic materials, can undergo not only a diffusive regime but also a hydrodynamic regime (Poiseuille-like) and a ballistic one. In the first case, heat transfer can be described by means of Fourier’s law. In the latter cases, instead, Fourier’s law breaks down [4,11–17], and more general constitutive equations for the local heat flux are needed [18,19].

Different theories and/or approaches can be found in the literature to face with the breakdown of Fourier’s law at the nanoscale and to describe thermal transport in nanosystems [1,11,20–25]. Here, we focus our attention on phonon hydrodynamics [4,13,26], since it connects mesoscopic and microscopic approaches.

It is worth noting that the term *phonon hydrodynamics* may be interpreted in two different ways:

(1) as a particular regime of phonon flow, wherein normal phonon collisions and collective effects are dominant, thus leading to a Poiseuille-like flow;

(2) as a more general thermodynamic model of heat transfer, wherein non-local effects play a relevant role [13,26–33] in such a way that, for very small systems, the form of the constitutive equation for the local heat flux is analogous to that for the velocity of a viscous fluid in classical hydrodynamics.

The usual starting point of phonon hydrodynamics (in the second meaning) is the so-called Guyer–Krumhansl (GK) equation [27,28], which is characterized by the inclusion of second-order non-local effects in the evolution equation for the heat flux. One version of that equation is
*τ*_{R} is the relaxation time related to resistive (i.e. momentum non-conserving) collisions between phonons, *T* is the temperature, λ is the usual Ziman limit for the bulk thermal conductivity [34], ℓ is an average phonon mean-free path, and *a* (*T*) means a dimensionless scalar function [35].

In the hydrodynamic regime (i.e. when in equation (1.1) the contribution of the heat flux may be neglected with respect to its spatial derivatives [26,29,30]), the phonons exhibit a macroscopic drift motion [4], which is the main difference between that regime and the more well-known diffusive and ballistic regimes. Due to that drift motion, in the phonon–hydrodynamic regime, the heat-flux profile is not uniform everywhere and may exhibit a parabolic profile analogous to that of Poiseuille hydrodynamic flow.

In steady states, in [26,29,30], equation (1.1) has been used to describe the size-dependency of the effective thermal conductivity (ETC) in nanosystems of different sizes and shapes. In doing so, it was assumed that the local heat flux is the sum of two different contributions: the bulk contribution **q**_{b}, arising from the solution of equation (1.1) with vanishing boundary conditions, and a wall contribution **q**_{w}, which may be related to **q**_{b} by a constitutive equation of the form
*ξ* means the normal direction to the wall cross section (pointing towards the flow), and *γ* is the curve accounting for the outer surface of the transversal section of the system. Moreover, *C* is a positive constant related to the properties of the walls, namely to the combination of diffusive and specular phonon–wall scattering [26,29,36–38] as
*p* being the relative number of phonon–wall specular collisions as compared with the total number of specular and diffusive collisions. This is well known for gases and for electron collisions [34] and it can also be adopted for the phonons [35]. Since *p* may only vary in the range [0;1], the parameter *C* should be larger than 2. Small values of *C* mean that diffusive phonon–wall collisions are predominant over specular ones. Conversely, large values of *C* mean that the specular phonon–wall scattering dominates over the diffusive one. The relation between the parameters *a*, *C* and ℓ in the boundary condition (1.2) has been investigated in [35]. In [26,29,30], equations (1.1) and (1.2) were used to provide a phenomenological way to explicitly describe the effects of phonon–wall collisions [5,39,40] through a formalism which is parallel to that of rarefied fluid dynamics [36–38].

In principle, the wall contribution **q**_{w} is restricted to a thin region near the walls, the so-called Knudsen layer, whose thickness is of the order of the mean-free path of the heat carriers. Far from the walls, **q**_{w} is vanishingly small [26,29,30]. The use of equation (1.2) allows us to account for the boundary scattering, which is the main cause of the non-uniform heat-flux profile in the hydrodynamic regime. We note that, from the theoretical point of view, the approach used in [26,29,30] prescribes that only the bulk contribution **q**_{b} can be considered as a state-space variable.

An alternative way to account for the boundary scattering by means of a slip heat flux is to assume that the local heat flux **q**, instead of **q**_{b}, is the solution of equation (1.1) together with the following Robin-type boundary condition:

In this way, one can assume that **q** (and not its partial contribution **q**_{b}) belongs to the state spaces, as is usual in extended thermodynamics [13–15]. In this paper, we principally investigate the alternative method described above.

As a brief summary, the paper is organized as follows.

In §2, we use equation (1.1), complemented by equation (1.3), to show how phonon hydrodynamics may be used to investigate heat transfer in two-dimensional narrow strips in a wide temperature range (or Knudsen-number range). Therein, we also derive the ETC and compare the theoretical expression with that obtained in [26,29].

In §3, we derive the theoretical behaviour both of the local heat flux and of the ETC by means of an approach closer to the original proposal by Guyer & Krumhansl [27,28], which directly introduces the phonon–wall interactions in the differential equation (1.1) and assumes **q**=**0** on the walls. In other words, in this section, we assume that in equation (1.1), λ is a size-dependent resistive thermal conductivity, and the quantity *a*ℓ in the non-local term is related to normal phonon–phonon collisions, the dynamical role of which has been recently stressed in a microscopic kinetic-collective model of heat transport in [1,2].

In §4, the main results are finally summarized.

## 2. Phonon hydrodynamics in two-dimensional narrow strips

In the present section, we apply equations (1.1) and (1.3) to study the thermal transport in a two-dimensional strip (or a very thin layer). In particular, we assume that the thickness of the layer is much smaller than its other two characteristic sizes. This is the case, for example, for graphene, which is a flat monolayer of carbon atoms tightly packed into a two-dimensional honeycomb lattice [8]. For the sake of illustration, in figure 1, a sketch of the system we are analysing in this paper can be seen. Therein, the dashed lines indicate a generic cross section *yz* which, in the case of two-dimensional nanosystems, reduces to a line along the *y*-axis.

In particular, we assume steady states in such a way that, from the local balance of internal energy per unit volume *u* in the absence of heat source, namely
*u*, one has ∇⋅**q**=0 and, consequently, equation (1.1) reduces to

If we assume that **q** may only flow along the *x*-axis, and that the temperature gradient is constant along the layer, the general solution of equation (2.2) is
*T*/*L*=−∇*T*, and *A* and *B* are two arbitrary integration constants. Moreover, in this case, the boundary condition (1.3) becomes
*C*_{1}≡*a*(*T*)*C*, and *w* is the width of the strip such that −*w*/2≤*y*≤*w*/2. For symmetry reasons, if we use the further boundary condition
*yz* plane), then the coupling of equations (2.3)–(2.5) yields
*w* is the Knudsen number.

### (a) Results for the local heat-flux profile

In figure 2, we plot the behaviour of the non-dimensional heat flux *q*_{*}=*q*/(λΔ*T*/*L*) for different values of the Knudsen number (namely, Kn=0.01, corresponding to the line with the right-pointing triangle markers; Kn=0.1, corresponding to the line with the square markers; Kn=0.5, corresponding to the line with the circle markers; Kn=1, corresponding to the line with the upward-pointing triangle markers; and Kn=10, corresponding to the line with plus-sign markers), and for different values of the non-dimensional product *C*_{1}≡*a* (*T*) *C* (i.e. *C*_{1}=0.05, *C*_{1}=0.5, *C*_{1}=2 and *C*_{1}=10). The results in figure 2 were obtained from equation (2.6).

On the vertical axis of each subfigure, the dimensionless quantity *y*_{*}≡*y*/*w* spans the range [−0.5;0.5]. In other words, the total width of the strip is taken as the reference length, to which the other lengths are referred. This means that changes in Kn may be due to changes in *w* at a constant ℓ, or to changes in ℓ at a constant *w*. Since the phonon mean-free path ℓ depends on temperature (in general, the smaller the *T*, the larger the ℓ), we may therefore equivalently claim that *q*_{*} in figure 2 is plotted as a function of the average temperature of the system, for a given *w*. Thus, the approach based on equations (1.1)–(1.3) also allows us to investigate the behaviour of the heat flux in a wide temperature range.

Depending on Kn, the different characteristic behaviours of the local heat flux can be seen in figure 2. Their main features are discussed below.

#### (i) Low values of Kn (Fourier diffusive regime)

When the Knudsen number is small, the regime of heat transfer is diffusive. As can be seen from figure 2, in such a case (see, for example, the curves corresponding to Kn=0.01) the heat-flux profile has a shape which is almost uniform across the strip, whatever the values of *C*_{1} are. In particular, if

The uniform profile, which is that predicted by the classical Fourier’s law, characterizes the diffusive regime. In the diffusive regime, the phonons undergo multiple scattering inside the strip, since their mean-free path is much smaller than the characteristic size, but they do not have a macroscopic drift motion [4]. The heat-flux profile is uniform because the umklapp phonon–phonon collisions, as well as the phonon-impurity and the phonon-defect collisions (which are felt everywhere in the bulk of the system), are prevalent over the phonon–boundary scattering (which, instead, mainly occurs near the boundary of the system, that is, in the so-called Knudsen layer [26,29]).

We explicitly note that for low values of Kn, from figure 2, it is possible to see a slight discrepancy between the results predicted by equation (2.6) and those arising from the strict use of a Fourier-like equation. The latter equation, in fact, predicts a completely flat profile, whereas we have, instead, a flat profile except for a very narrow zone near the walls, where there is a steep increase in the heat-flux profile from a low value to the flat value. This mathematical discrepancy is not a truly physical discrepancy, because in strict terms Fourier’s law is only valid for spatial scales larger than the phonon mean-free path. Thus, since for small values of Kn, the phonon mean-free path is much smaller than the width of the strip, the discrepancy is not only quantitatively negligible in reference to the total heat flux, but it is also conceptually admissible because it is found out of the range of validity of Fourier’s law. Anyway, it is an interesting question which should be explored in greater detail in computer simulations. In fact, such a steep increase in the heat flow near the walls has been modelled on other grounds than those described here (i.e. without incorporating the Laplacian term in equation (1.1)) by assuming that the thermal resistance is higher than that predicted by Fourier’s law near the walls, by a factor of *D* being the diameter of the channel [41,42]. The subtle details of this simplified, but efficient, model and our model should be considered in a more in-depth way.

#### (ii) Moderate values of Kn (Poiseuille phonon flow)

For moderate values of the Knudsen number, the heat-transfer regime is in transition from the diffusive regime to the ballistic one. In particular, for increasing values of Kn (i.e. whenever ℓ becomes comparable to, or larger than, the width *w* of the strip), the phonon bulk scattering gradually reduces, and the phonon–wall scattering increases. Moreover, the phonons show a macroscopic drift motion. In this case, the boundary–phonon scattering becomes prevalent among the other scattering mechanisms. The heat flux is smaller near the walls because the boundary–phonon scattering is mainly felt in their neighbourhood. This explains the parabolic profile of the local heat flux we can see, for example, in figure 2 when Kn=0.1, Kn=0.5 and Kn=1. Since in this case, the higher the Kn, the smaller the maximum value of *q*, we may infer that phonon–boundary collisions are the main mechanism of momentum loss and resistance to the heat flow, especially for small values of *C*_{1} (i.e. for diffusive phonon–wall collisions predominating over specular ones). Instead, for higher values of *C*_{1}, the parabolic profile is not so evident, as the slip flow makes them closer to a flat profile.

#### (iii) Large values of Kn (ballistic regime)

When Kn becomes large (for example, when Kn≈10 in the case of the study described in this paper), the regime of heat transfer becomes ballistic. In this regime, the phonons suffer only scant internal (i.e. in the bulk of the system) scattering. As can be seen from figure 2, in this case the heat flux again attains a uniform profile, similar to that observed in a diffusive regime. In particular, if

It is worth noting that, in contrast to the flat heat-flux profile, we find that, for very low values of the Knudsen number (i.e. when equation (2.7) holds), the flat heat-flux profile arising from equation (2.8) depends both on the width of the strip, for a given ℓ (i.e. for a given temperature) through the Knudsen number, and on *C*_{1}, as it is logical to expect.

From the physical point of view, equation (2.8) can be explained in the following way. When Kn reaches high enough values, the Knudsen layer pervades the main part of the cross section. In this case, the frequency of phonon–boundary interactions is so high that phonon–boundary scattering is the main event, and the consequent reduction in the local heat-flux amplitude due to the momentum loss becomes again homogeneous in the whole transversal section. It is worth observing that, in [26,29,30], the authors draw the same conclusion, but with an alternative way of accounting for the boundary conditions (namely, by using equation (1.2) instead of equation (1.3)). Let us also observe that the higher the *C*_{1}, the smaller the values of Kn at which this second limit behaviour occurs.

As interpreted from the point of view of the width *w*, the three regimes mentioned above are as follows: diffusive (Fourier transport) is found for *w*>ℓ_{R}, with ℓ_{R} being the mean-free path of resistive collisions; the Poiseuille phonon flow for ℓ_{N}<*w*<ℓ_{R}, with ℓ_{N} being the mean-free path of normal collisions; and the ballistic flow for *w*<ℓ_{N},ℓ_{R}.

### (b) Results for the effective thermal conductivity

The ETC of narrow strips is much reduced with respect to that of wide strips of the same material [43]. Since the thermal conductivity of solid thin layers is important for the design of field-effect transistors in electronic circuits, coated lenses in laser systems and microfabricated superconducting radiation detectors [44,45], then theoretical model predictions it may be useful in practical applications.

Starting from the results obtained in §2a, here we derive a theoretical model for the ETC in thin layers, defined as
*A* is the area of its cross section and *Q*_{tot} is the total heat per unit of time flowing in the system. Indeed, in the two-dimensional case, instead of a cross section *A*, only the width *w* has to be used in the denominator of equation (2.9).

Combining equations (2.6) and (2.9), one is led to the following expression for the ETC in narrow strips (or very thin layers):

The theoretical prediction in equation (2.10) differs from the expression
^{−1} for increasing Knudsen number in the ballistic regime [5,26,29], i.e. when Kn≫1. Similarly, for very low values of Kn (i.e. in the Fourier diffusive regime), both equations (2.10) and (2.11) turn out to be λ_{eff}=λ_{eff,diff}≡λ.

By means of equation (2.10), in figure 3, we plot the behaviour of the non-dimensional ETC λ_{*}=λ_{eff}/λ in a two-dimensional strip as a function of the Knudsen number Kn=ℓ/*w* for different values of the non-dimensional product *C*_{1}≡*a* (*T*) *C* (i.e. *C*_{1}=0.05, *C*_{1}=0.5, *C*_{1}=2 and *C*_{1}=10) accounting for the diffusive and specular phonon–wall interactions.

As can be seen from figure 3, in the ballistic regime (i.e. when Kn≫1), the ETC is much reduced with respect to that characterizing the diffusive regime (i.e. when Kn≪1).

## 3. Comparison with the Guyer–Krumhansl formalism

Though we previously referred to equation (1.1) as the GK equation, it has in fact important differences with respect to the original proposal of the authors [27,28], especially regarding the way in which the phonon–wall collisions are included in the model. In the approach used in §2 (as well as in that of [26,29]), in the differential equation (1.1), the usual bulk thermal conductivity obtained as the usual Ziman limit [34] is used, that is, *c*_{v} is the specific heat at constant volume per unit volume, *τ*_{R0} is the characteristic time of resistive phonon collisions in the bulk, given by
*τ*_{u} is the relaxation time of umklapp phonon–phonon collisions, *τ*_{i} is the relaxation time of phonon–impurity collisions and *τ*_{d} is the relaxation time of phonon–defect collisions. The phonon–wall interactions have then been accounted for by suitable boundary conditions [24,26,29], which in the present case are given by equation (1.3).

Indeed, in [28], the authors considered a boundary relaxation-time *τ*_{b} that they combined with the usual relaxation time *τ*_{R0} due to bulk resistive mechanisms by means of the Matthiessen rule as
*b*(*T*) is a suitable non-dimensional temperature function depending on the form of the cross section of the system and on the relative abundance of specular and diffusive phonon–wall collisions.

Once the combined resistive-boundary (phonon–wall) collision time had been obtained, the thermal conductivity λ_{GK} (depending on the size of the system through *τ*_{b}) was calculated, and used in the first term on the right-hand side of equation (1.1). In particular,

It is also important to note that, in equation (1.1), ℓ means the average phonon mean-free path obtained from the Ziman expression for the bulk thermal conductivity, i.e. *a*^{2}ℓ^{2}, *τ*_{N} is the characteristic collision time of normal (momentum conserving) phonon–phonon collisions. In terms of *a*^{2}ℓ^{2}, the coefficient *a*^{2} could be obtained by equating *a*^{2}ℓ^{2}=ℓ^{2}_{GK}, namely *a*^{2} (*T*)=*τ*_{N} (*T*)/5*τ*_{R} (*T*).

A further observation is that, in the GK approach, the non-slip boundary condition for **q** is considered, namely **q**|_{γ}=0, and the effects of the boundary collisions are considered in the form (3.3) of the ETC. The term in ∇^{2}**q** is supposed to describe only an additional contribution of the walls restricted to the case when phonons display a collective behaviour when normal collisions are dominating with respect to resistive ones [1,2].

This alternative approach, in steady states, yields the following behaviour for the local heat flux across any cross section of the strip:

When

For small Kn, this is the Fourier result in equation (2.7), whereas for high Kn it is analogous to equation (2.8) if *b* (*T*) is identified as *b* (*T*)=[*a* (*T*) *C*]^{−1} (for a different geometry, the relation between *b* (*T*) and *a* (*T*)*C* would also have a form-dependent numerical factor). When

Finally, by the definition of ETC in equation (2.9) in this case, we have

In figure 4, we compare the results arising from our phonon–hydrodynamic approach with those arising from the GK formalism. The comparison involves both the heat-flux profile (figure 4*a*), and the ETC (figure 4*b*). To compare the heat-flux profile, we used equations (2.6) and (3.4) in the case of Kn=0.2. To compare the ETC we used, instead, equations (2.10) and (3.7). For the sake of computation, we assumed *a*=1, *C*=5 and *b*=*C*^{−1}≡0.2.

As can be seen from figure 4, both for *q*_{*}, and for λ_{*} the phonon–hydrodynamic approach predicts larger values than the GK formalism.

In closing the present section, it is important to stress that both the theoretical heat-flux profile in equation (2.6) and the theoretical ETC in equation (2.10) are strictly related to the non-dimensional parameters *a* (*T*) and *C*. Although in [26,30,35], for example, more in-depth investigations of their physical interpretations can be found, the very important role these parameters play requires further study. To this end, it would be useful to infer information about them from the comparison of our theoretical model with other existing models, experimental data or simulation results.

## 4. Conclusion

In this paper, we have examined the heat-flux profile across a narrow two-dimensional strip as a function of the Knudsen number Kn≡ℓ/*w* using two different forms of a generalized heat-transport equation containing second-order non-local terms. From the model based on equations (1.1)–(1.3), we have found that the heat-flux profile in a generic transversal section is given by equation (2.6). In particular, that profile has a flat shape given by equation (2.7) for Kn≪1 (i.e. in the Fourier regime), a parabolic one at intermediate Kn (between 0.1 and 1, for instance, depending on the values of the slip coefficient *C*), and again a practically flat profile as in equation (2.8) for high Kn (i.e. in the ballistic regime). Thus, the Poiseuille flow with a parabolic profile is found only in a given temperature range, or a given size range for the strip width.

In the original GK model based on equations (1.1) and (3.3), the parabolic heat profile is also found in some restricted temperature (or size) range. Conceptually, both formalisms of non-local heat transport agree on this specific range for hydrodynamic phonon flow.

The most interesting difference between both models rests on the boundary conditions used for **q** on the walls, and on the way of accounting for the phonon–wall collisions. In the model based on equations (1.1) and (1.3), the non-vanishing boundary conditions imply a slip heat flow along the walls, whereas in the model based on equations (1.1) and (3.3) non-slip conditions are imposed on the heat flux. Furthermore, in the approach based on equations (1.1) and (3.3), the size-dependent thermal conductivity (3.3) describes the phonon–wall collisions, whereas in the description based on equations (1.1) and (1.3) the role of the shape of the system is more explicit (through the use of equation (1.3), and of the relative role of specular or diffusive collisions). In this case, the presence of a slip flow also connects in a smooth way the hydrodynamic–phonon regime to the ballistic regime.

In summary, a possible way of checking how phonon–wall effects should be more faithfully included in a generalized heat-transport equation with non-local terms would be to check whether the heat flux exhibits some slip effects along the walls, especially in a strip of suspended graphene, where the hydrodynamic regime is expected to show at a relatively wide range of temperatures around 100 *K*, in contrast to bulk materials, where it is shown at low temperatures (1–2 *K*) in a narrow temperature range [4]. If there is such a slip flow, the dependence of the ETC on the width *w* of the layer would not be on *w*^{2}, as in Poiseuille flow analogy, but is slightly different because the profile would be slightly flatter. This would not mean, however, that collective hydrodynamic effects are not present, but that their manifestation is somewhat obscured by the effects of the slip flux.

Finally, it is worth noting that here we have considered narrow but long strips, in such a way that only boundary effects on the lateral walls have been taken into account. In the case of short strips, instead, the effects on the entrance and exit boundaries have to also be taken into account. The physical mechanism for these boundary effects is different from that on the lateral walls [46,47], as is manifested, for instance, in the anisotropy in-plane and cross-plane thermal conductivity in narrow films [46]. Therefore, in the future it would be interesting to consider both these aspects, and the so-called *entrance effects* which are found in the hydrodynamic analyses of short channels.

## Data accessibility

The non-dimensional heat flux *q*_{*}=*q*/(λΔ*T*/*L*) and the non-dimensional ETC λ_{*}=λ_{eff}/λ, which we plotted in figures 2–4, allowed us to avoid the use of experimental data. All the results of the present article are fully reproducible by means of the data we provided here.

## Author' contributions

All authors contributed equally.

## Competing interests

All authors have no competing interests.

## Funding

A.S. acknowledges the financial support of the Italian Gruppo Nazionale per la Fisica Matematica (GNFM-INdAM). I.C. acknowledges financial support from the University of Basilicata and from the University of Salento for her stay in the Autonomous University of Barcelona as part of her doctoral studies. D.J. acknowledges the financial support of the Spanish Ministry of Economy and Competitiveness under grant no. FIS 2012-33099, and the Ministry of Science and Innovation under the Consolider Project Nano Therm (grant no. CSD-2010-00044).

- Received June 6, 2015.
- Accepted August 25, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.