## Abstract

The travel of heat in insulators is commonly pictured as a flow of phonons scattered along their individual trajectory. In rare circumstances, momentum-conserving collision events dominate, and thermal transport becomes hydrodynamic. One of these cases, dubbed the Poiseuille flow of phonons, can occur in a temperature window just below the peak temperature of thermal conductivity. We report on a study of heat flow in bulk black phosphorus between 0.1 and 80 K. We find a thermal conductivity showing a faster than cubic temperature dependence between 5 and 12 K. Consequently, the effective phonon mean free path shows a nonmonotonic temperature dependence at the onset of the ballistic regime, with a size-dependent Knudsen minimum. These are hallmarks of Poiseuille flow previously observed in a handful of solids. Comparing the phonon dispersion in black phosphorus and silicon, we show that the phase space for normal scattering events in black phosphorus is much larger. Our results imply that the most important requirement for the emergence of Poiseuille flow is the facility of momentum exchange between acoustic phonon branches. Proximity to a structural transition can be beneficial for the emergence of this behavior in clean systems, even when they do not exceed silicon in purity.

## INTRODUCTION

The finite thermal resistivity of an insulating solid is a manifestation of the anharmonicity of the underlying lattice. The most common approach to calculate the thermal conductivity of a given solid is to assume the harmonic approximation and introduce a finite lifetime for phonons that captures the effects beyond the harmonic approximation. In the Boltzmann-Peierls picture, heat-carrying phonons are scattered by other phonons or by crystal imperfections and boundaries. The rate of collisions relaxing the momentum of the traveling phonon set the magnitude of thermal resistivity. Only a subset of phonon-phonon scattering events, called a Umklapp (or U) process, degrades the heat current. The initial and final momenta in a U process differ by a multiple of a reciprocal lattice vector. On the other hand, a normal (or N) phonon-phonon collision cannot lead to thermal resistance in the absence of Umklapp scattering. In an infinite sample, if all collisions were normal, then one expects the flow of phonons to be undamped in the absence of any external field applied to sustain it (*1*).

Black phosphorus (BP) has attracted much recent attention as a cleavable solid with a promising exfoliating potential toward two-dimensional phosphorene (*2*, *3*). It has an orthorhombic crystal structure with puckered honeycomb layers in the *ac* planes and van der Waals coupling between the layers (Fig. 1, A and B). Unlike graphene, it has a tunable direct bandgap ranging from 0.3 eV in bulk to 2 eV in a monolayer (*4*). These features make BP a promising material for applications. In addition, the presence of a significant in-plane anisotropy may induce spatially anisotropic transport (*5*) absent in other graphene-like materials. Although electronic conduction in BP has been extensively investigated (*6*–*9*), few studies have been devoted to its thermal transport (*10*–*15*). They are restricted to relatively high temperatures and did not explore the low-temperature range, the focus of the present study.

Here, we report on a study of thermal conductivity along the *a* and *c* axes of BP single crystals and establish the magnitude and temperature dependence of the thermal conductivity down to 0.1 K. We find a moderately anisotropic thermal conductivity mainly reflecting the anisotropy of the sound velocity. Unexpectedly, we find that below its peak temperature, the thermal conductivity evolves faster than the ballistic *T*^{3}. This feature, combined with size dependence of the Knudsen minimum in the effective phonon mean free path, provides compelling evidence for Poiseuille phonon flow, a feature previously restricted to a handful of solids (*16*). We argue that this arises because of the peculiar phonon dispersions enhancing the phase space for normal scattering events in BP. This observation has important consequences for the ongoing research in finding the signatures of hydrodynamic phonon flow in crystals.

## RESULTS

Our samples are insulators, as seen in Fig. 1C, which shows the temperature dependence of the electrical resistivity ρ along the *a* and *c* axes. For both directions, ρ shows a peak around 250 K, reflecting the change in the number of thermally excited carriers across more than one gap. The magnitude of the resistivity, its temperature dependence, and its anisotropy are in agreement with an early study (*6*). Below 50 K, ρ shows an activated behavior, as seen in the upper inset of Fig. 1C. The distance in the energy between the top of the valence band and the acceptor levels set the activation energy of hole-like carriers (*6*). Using the Arrhenius equation ρ = exp(*E*_{g}/2*k*_{B}*T*) between 40 and 20 K, one finds an activation energy of *E*_{g} ~ 15 meV, in agreement with the previous result (*6*). Below 20 K, the temperature dependence of ρ becomes appreciably weaker, indicating that the system enters the variable range hopping (VRH) regime where the electrical transport is governed by the carriers trapped in local defects hopping from one trap to another. One can describe the resistivity by the expression ρ ∝ exp[(*T*/*T*_{0})^{−1/(d+1)}] in a system with dimension *d*. As seen in Fig. 1 (D to F), it is hard to distinguish the hopping space dimensionality from our results. The resistivity anisotropy ratio, ρ_{a}/ρ_{c} (lower inset of Fig. 1C), shows little temperature dependence in the activated regime. The ratio, almost constant (~3.5) in the activated regime and reflecting the intrinsic anisotropy in the carrier mobility (*6*), sharply increases in the VRH regime and attains ~10.

Because of the extremely low electrical conductivity, phonons almost entirely transport heat, and we are going to focus on them. Figure 2A shows the temperature dependence of the thermal conductivity, κ, along the *a* and *c* axes in a BP single crystal. Four different groups recently measured the thermal conductivity of bulk BP (*12*–*15*). Sun *et al*. (*13*) found a large and anisotropic thermal conductivity, much larger than what was found in an early polycrystal study by Slack (*10*) and a single-crystal study by other groups (*12*, *15*). As one can see in the figure, our data for the samples of ~100-μm thickness (solid circles) is in good agreement with the data reported by Sun *et al*. (*13*), including the anisotropy essentially due to the sound velocity (see the Supplementary Materials for details). Both sets of data point to an intrinsic bulk conductivity much larger than what was found in a polycrystal (*10*). The large thickness dependence of the thermal conductivity in BP found by our study (see below), in agreement with what was recently reported by Smith *et al*. (*14*), provides a clue to the striking discrepancy between the single-crystal and polycrystal data. If, even at temperatures exceeding 100 K, the phonon thermal conductivity is affected by sample dimensions, then the attenuation of thermal conductivity in the presence of grain boundaries is not surprising.

In the absence of anisotropy, the lattice thermal conductivity can be represented by the following simple expression(1)where *C* is the phonon specific heat, 〈*v*〉 is the mean phonon velocity, and *l*_{ph} is the phonon mean free path. Note that the ^{1}/_{3} prefactor, a consequence of averaging over the whole solid angle in isotropic medium, can be different in the presence of anisotropy. At temperatures exceeding the Debye temperature, where the Dulong-Petit law applies, the specific heat saturates to a constant value. In this regime, the temperature dependence of the phonon mean free path, which becomes shorter with increasing temperature, governs the temperature dependence of the thermal conductivity. At the other extreme, that is, at low temperature, the mean free path saturates to its maximum magnitude, of the order of the sample dimensions, and the *T*^{3} behavior of the phonon specific heat sets the temperature dependence of thermal conductivity. Between these two extreme regimes, dubbed the kinetic (at high temperature) and the ballistic (at low temperature), the thermal conductivity passes through a peak. As demonstrated by extensive studies on LiF (*17*), introducing isotopic impurities damps κ_{max}, the peak value of the thermal conductivity. On the other hand, decreasing the size of the sample shifts the position of κ_{max} with little effect on the magnitude of thermal conductivity in temperatures exceeding the peak.

Our key finding is the observation of a κ evolving faster than *T*^{3} just below the peak temperature in BP. This can be seen in Fig. 2B. As a consequence, κ divided by *T*^{3} shows a maximum or a shoulder around 10 K (see the inset of Fig. 2A). Upon further cooling, the heat transport eventually becomes ballistic as evidenced by the *T*^{3} behavior of κ (saturation of κ/*T*^{3}) observed in the low-temperature data of the thickest sample (open circles in Fig. 2A). As a consequence, the effective phonon mean free path, *l*_{ph}, extracted from Eq. 1, presents a local maximum (around 10 K) and a local minimum (around 4 K; solid circles in Fig. 3, C and D). Here, we calculated *l*_{ph} using 〈*v*_{a}〉 = 0.536 × 10^{4} m/s, 〈*v*_{c}〉 = 0.354 × 10^{4} m/s [these are calculated velocities in agreement with previous calculations (*18*) and neutron scattering measurements (*19*)], and the measured specific heat of the crystal (see the Supplementary Materials for details). In milliKelvin temperatures, *l*_{ph} becomes comparable to the sample thickness (see the inset of Fig. 3C). Decades ago, the phonon Poiseuille flow (*16*) was detected in Bi (*20*) by observing these two features, namely, a faster than *T*^{3} evolution of κ (Fig. 2C) and a peak in the temperature dependence of *l*_{ph} just below the maximum thermal conductivity (Fig. 3E). Similar features have also been found in crystalline ^{4}He (*21*), ^{3}He (*22*), and H (*23*).

To exclude any experimental artifact, we measured the thermal conductivity of P-doped Si using the same experimental setup. As seen from Fig. 2A, the temperature dependence of κ for P-doped Si conforms to what was found previously (*24*). κ shows a cubic temperature dependence between 2 and 5 K. (Fig. 2D). To show this point more explicitly, we plotted *l*_{ph} that was calculated from Eq. 1 as a function of temperature in the inset of Fig. 3D. The physical parameters used in the calculation are θ_{D} = 674 K and 〈*v*〉 = 6700 m/s (*25*). As expected, *l*_{ph} does not show a local maximum, and at low temperature, its magnitude is in reasonable agreement with the sample dimensions (0.3 mm × 1.4 mm × 4.0 mm).

Only Umklapp and impurity scattering events degrade the momentum of the traveling phonons. Normal scattering events conserve crystal momentum and do not contribute to the thermal resistance. When normal scattering is sufficiently strong and the momentum-degrading scattering events are significantly rare, normal scattering can even enhance the heat flow (*1*, *16*). The essential condition for this is given by the inequality(2)

Here, *d* is a characteristic sample dimension, and *l*^{N} and *l*^{R} are normal and resistive mean free paths, respectively. In such a situation, the phonons lose their momentum only by diffuse boundary scattering, and they flow freely while continuously exchanging momentum in the center of substances, analogous to the Poiseuille flow in fluids (Fig. 2E). As a result, moving like a Brownian particle, a phonon traverses the path length of the order of *l*_{ph} ~ *d*^{2}/*l*^{N} between successive collisions with the boundaries, and in the ideal case, *l*_{ph} can reach a value larger than the characteristic sample dimension *d*. If *l*^{N} increases with decreasing temperature as *T*^{−5}, then since κ ~ *C*〈*v*〉*d*^{2}/*l*^{N}, one can expect thermal conductivity to follow as *T*^{8} (*1*). Experimentally, in all systems in which Poiseuille flow has been reported (namely, Bi, solid He, and H), what has been observed is a faster than *T*^{3} thermal evolution of κ and a nonmonotonic *l*_{ph}(*T*) (*20*–*23*), not a superlinear size dependence. This is also the case in the present study on BP.

In the narrow temperature window of the Poiseuille flow, the dominance of normal scattering creates an unusual situation where warming the system enhances the mean free path because it enforces momentum exchange among phonons. In other words, the temperature dependence of thermal conductivity is set by the temperature dependence of the phonon viscosity and not exclusively by the change in the population of thermally excited phonons or the frequency of momentum-relaxing events, as it is the case in the ballistic and diffusive regimes. At sufficiently low temperatures, the phonons start to behave just like a Knudsen flow of gas (*26*). This transition from the Poiseuille flow to the Knudsen flow, *l*_{ph} is accompanied by a mean-free-path minimum, dubbed Knudsen minimum, which occurs when *d*/*l*^{N} ~ 1. This is where the diffuse boundary scattering rate is effectively increased because of normal scattering (Fig. 4A).

By changing the thickness of the samples along the *b* axis, we also examined a size dependence of the thermal conductivity and confirmed that the effect is more prominent in larger samples (Fig. 3, A and B). As expected, in thicker BP samples, the maximum *l*_{ph} is longer. Moreover, as seen in the case of Bi (Fig. 3E) (*20*), both the local maximum and the local minimum in *l*_{ph}(*T*) persist (Fig. 3, C and D). In addition, note that despite the anisotropic thermal conductivity, *l*_{ph} attains values comparable to the sample dimensions at low temperatures, irrespective of the orientation of the heat current, meaning that *l*_{ph} is quasi-isotropic along the three principle axes and is set by the average sample size. We derived this from the low anisotropy in the sound velocity between the *ac*-plane and out-of-plane orientations despite the layered structure of BP (see fig. S4 and table S2). The peak thermal conductivity in larger BP crystals attains a magnitude as large as 800 W/Km, much higher than the values found in thinner samples (Fig. 3, A and B) (*14*) but still smaller than the maximum thermal conductivity of large (~cm) crystals of Si (3000 W/Km) (*24*) or Bi (4000 W/Km) (*20*).

In the original Gurzhi (*1*) picture, the quadratic dependence of *l*_{ph} on the sample dimension, *d*, would lead to a superlinear size dependence of the thermal conductivity in the Poiseuille regime. As one can see in Fig. 5B, this is absent in BP, as it was in the case of Bi (*20*). Now, the Poiseuille flow of phonons is expected to emerge concomitantly with the second sound, a wavelike propagation of entropy (*16*, *27*). Both require a peculiar hierarchy of scattering events (*16*, *28*), where the normal scattering rate is much larger than the boundary scattering and the latter much larger than resistive (Umklapp and impurity) scattering. Its experimental observation in Bi in the same temperature range (*29*) confirms the interpretation of static thermal conductivity data.

The absence of superlinear size dependence can be traced back to the expression for phonon kinematic viscosity, proposed by Gurzhi (*1*), *ν* = *v*_{T}*l*^{N}, which implies that the phonon fluid is non-Newtonian. Here, *v*_{T} is the phonon thermal velocity whose magnitude is comparable to the velocity of the second sound, so that *v*_{T} ~ 〈*v*〉/3 = 2000 to 3000 m/s and *l*^{N} is estimated to be of the order of 10 μm at the Knudsen minimum. This yields *ν* ~ 0.02 to 0.03 m^{2}/s, several orders of magnitude larger than the kinematic viscosity of water. This comparison is to be contrasted with the case of electrons in PdCoO_{2} (*30*) whose dynamic viscosity was found to be comparable to the dynamic viscosity of water. Since the phonon viscosity depends on the normal scattering rate and local velocity, it is not a constant number at a given temperature. The velocity profile of such a non-Newtonian fluid is much flatter than parabolic (*31*). This makes the size dependence much less trivial than in the parabolic and Newtonian case. A recent theoretical study (*32*) has shown that the thickness dependence can be either sublinear or superlinear according to the relative weight of boundary and normal scattering rates.

There is, however, another signature of the Poiseuille flow in the thickness dependence, which our experiment detected. In the Poiseuille regime, boundary scattering is rare. Therefore, with increasing thickness, phonons can travel a longer distance between successive collisions with the boundaries, giving rise to an enhancement of the thermal conductivity and the maximum in *l*_{ph} with the thickness. On the other hand, at the onset of the ballistic regime, the fraction of phonons that contribute to the diffuse boundary scattering due to numerous normal scattering increases with the thickness. As a consequence, the Knudsen minimum shifts to lower temperatures in the thicker samples (Fig. 4A). One can see both these features in a plot of *l*_{ph} normalized by its minimum value , clearly along the *c* axis (Fig. 4C) but not clearly along the *a* axis (Fig. 4B), as it was in the case of Bi (Fig. 4D) (*20*). We conclude that the Poiseuille flow is most prominent along the *c* axis. This suggests that the relative weights of normal and resistive scattering rates depend on orientation and one needs to go beyond a Boltzmann picture with a scalar scattering time (*33*).

The Ziman regime is another hydrodynamic regime (*16*). Here, normal scattering still dominates, but resistive scattering becomes larger than boundary scattering. In this regime, expected to occur just above the peak temperature, one expects κ to decrease exponentially with increasing temperature before showing a *T*^{−1} temperature dependence at higher temperatures. Our measurements do not detect a regime where the thermal conductivity follows an exponential behavior in BP. Instead, all samples showed a robust *T*^{−1} behavior in a wide temperature window extending down to 40 K, an order of magnitude lower than the Debye temperature (Fig. 5A). Moreover, the slope of the *T*^{−1} temperature dependence was larger in the thicker samples (see the inset of Fig. 5A).

Changing the thickness from 6(15) to 20(30) μm of the *a*(*c*)-axis sample is to multiply the number of parallel heat-conducting phosphorene layers by a factor of 3(2). Surprisingly, while no change in phonon dispersion is expected to occur in such a context, the enhancement in thermal conductivity persists up to 80 K, the highest temperature of our range of study (Fig. 5B). Such a large thickness dependence [first reported in submicronic samples by Smith *et al*. (*14*)] is unusual, as shown by a comparison with the much more modest effect observed in Ge (*34*). In Bi, a comparable enhancement is seen in its Poiseuille regime (Fig. 5B), but it rapidly dies away with heating. A plausible explanation is that in BP, some heat carriers remain ballistic at high temperature. Previous theoretical treatments of the thermal conductivity in graphene and graphite have argued for the presence of normal processes up to room temperature leading to collective phononic modes with a mean free path on the order of hundreds of micrometers (*28*, *35*, *36*). Our observation may also indicate the presence of substantial normal processes above the peak temperature of thermal conductivity in BP.

## DISCUSSION

Our theoretical calculations confirm that phonon-phonon scattering in BP is large at low energies. Figure 6 shows the phonon dispersions of BP and Si, which agree well with previous calculations (*37*, *38*). The phonon dispersions of Si show highly dispersive acoustic branches that lead to high thermal conductivity. The acoustic branches of BP are less dispersive because of a relatively weaker bonding than the one in Si, which has strong sp^{3} bonds. This leads to a higher phonon density of states (PDOS) in the low-frequency region, which allows for a larger phase space for the momentum-conserving normal three-phonon scattering events in BP.

Figure 6 (G and H) shows the calculated PDOS of BP and Si. The PDOS of BP is significantly larger than that of Si in the low-frequency region below 100 cm^{−1} as a consequence of the smaller dispersion of the acoustic branches. The phonon scattering processes involving the linearly dispersive acoustic modes conserve momentum and are not resistive to heat flow. The large low-frequency PDOS in BP provides much more phase space for these nondissipative scattering events. This can favor the emergence of Poiseuille flow by establishing the required hierarchy of time scales as discussed above.

In summary, by detecting that phonons present a faster than *T*^{3} thermal conductivity in a finite temperature window, we have identified a hallmark of Poiseuille flow of phonons in BP, which was previously seen only in four other solids (*39*). Theoretical calculations indicate that the peculiarities of phonon dispersion lead to a significant enhancement of the phonon-phonon normal scattering in BP compared to an archetypal insulator such as Si. We note that BP is close to a structural instability. Like Bi and other column V elements, its crystal structure results from a slight distortion of the cubic structure (*40*, *41*). Our results indicate that hydrodynamic phonon flow can be observed in a system close to a structural instability (*39*).

## MATERIALS AND METHODS

Single crystals of BP were synthesized under high pressure (*42*). Magnetophonon resonance, which requires reasonable purity, was observed in a sample from the batch used in the present experiments (*7*). A single crystal of P-doped Si was provided by the Institute of Electronic Materials Technology (Warsaw). The electrical resistivity and the thermal conductivity measurements were performed along the *a* and *c* axes of BP. Each sample has a rectangular shape with edges parallel to the three high-symmetry axes: the *a*, *b*, and *c* axes. Lengths of the samples are summarized in table S1. The *a*(*c*)-axis sample has the largest length along the *a*(*c*) axis, respectively. The thickness dependence of the thermal conductivity was investigated by using the same sample. The thickness (number of layers) along the *b* axis was decreased during the course of the investigations by cleaving.

The thermal conductivity was measured by using a home built system. We used a standard one-heater-two-thermometers steady-state method. A very similar setup was used to measure the thermal conductivity of cuprate (*43*) and heavy-fermion samples of a comparable size (*44*, *45*). The measurements were performed in the temperature range between 0.1 and 80 K. The thermometers, the heater, and the cold finger were connected to the sample by gold wires of 25-μm diameter. The gold wires were attached on the BP sample by Dupont 4922N silver paste and were soldered by indium on the P-doped Si sample. The contact resistances were less than 10 ohm at room temperature. The temperature difference generated in the sample by the heater was determined by measuring the local temperature with two thermometers (Cernox resistors in the ^{4}He cryostat and RuO_{2} resistors in the dilution refrigerator). The heat loss along manganin wires connected to the two thermometers and the heater is many orders of magnitude lower than the thermal current along the sample including the thinner ones (fig. S1). The heat loss by radiation is negligible in the temperature range of our study (*T* < 80 K), since it follows a *T*^{4} behavior. The heat loss by residual gas is also negligible because the measurements were carried out in an evacuated chamber with a vacuum level better than 10^{−4} Pa. Moreover, the external surface of the chamber is directly in the helium bath at 4.2 K for the measurement above 2 K, so that residual gas is condensed on the cold wall. The experimental uncertainty in the thermal conductivity arising from the measurements of the thermometer resistances is less than 1% in the whole temperature range. The main source of uncertainty results from uncertainty in the measured thickness of the samples, which is about 5% in the thinnest sample. The increase of thermal conductivity with the sample thickness, however, dominates over the experimental uncertainty.

The same contacts and wires were used for the electrical resistivity measurements by a four-point technique. A DC electric current was applied along the sample using the manganin wire attached to the heater. A Keithley 2182A nanovoltmeter was used for the measurement of electrical voltage. The input impedance of the nanovoltmeter is larger than 10 gigohms, which is well above the resistance of the sample even at the lowest temperatures (*R* ~ 5 megohms).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/6/eaat3374/DC1

Supplementary Text

table S1. Size of samples.

table S2. Phonon velocity along the high-symmetry axes.

fig. S1. Thermal resistance of BP samples and manganin.

fig. S2. Reproducibility of thermal conductivity data.

fig. S3. Specific heat.

fig. S4. Phonon dispersions at low energies.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We acknowledge R. Nomura for fruitful discussions.

**Funding:**This work was supported by the Japan Society for the Promotion of Science Grant-in-Aids KAKENHI 16K05435, 15H05884, and 17H02920; Fonds ESPCI; and the European Research Council grant ERC-319286 QMAC. The calculations were performed at the Swiss National Supercomputing Center, project s575.

**Author contributions:**Y.M. and K.B. conceived the project and planned the research. Y.M. performed transport measurements with the help of K.I. and analyzed data. K.A., A.M., and M.T. carried out the specific heat measurement. Y.A. synthesized samples. A.S. performed theoretical calculations. Y.M., A.S., and K.B. wrote the manuscript. All authors discussed the results and reviewed the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).