Research ArticlePHYSICS

Observation of Poiseuille flow of phonons in black phosphorus

See allHide authors and affiliations

Science Advances  22 Jun 2018:
Vol. 4, no. 6, eaat3374
DOI: 10.1126/sciadv.aat3374


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.


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 (69), few studies have been devoted to its thermal transport (1015). They are restricted to relatively high temperatures and did not explore the low-temperature range, the focus of the present study.

Fig. 1 Crystal structure and resistivity.

(A) Side view of the crystal structure of BP. (B) Top view of single-layer BP. (C) Temperature dependence of the electrical resistivity ρ of BP along the a and c axes in a logarithmic scale. The upper inset shows an activated behavior of ρ for T > 20 K, where solid lines represent an Arrhenius behavior. The lower inset shows the ratio ρac, which quantifies charge transport anisotropy. At low temperatures, VRH governs, as shown in the bottom panels where ρ is plotted against T−1/2 (D), T−1/3 (E), and T−1/4(F).

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 T3. 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.


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(Eg/2kBT) between 40 and 20 K, one finds an activation energy of Eg ~ 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/T0)−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, ρac (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 (1215). 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.

Fig. 2 Thermal conductivity.

(A) Temperature dependence of the thermal conductivity along the two high-symmetry axes for the samples with a comparable thickness (solid circles). The a-axis thermal conductivity measurement is extended down to 0.1 K using the thicker sample (open circles). We compared our data to three sets of previous reports (10, 12, 13). We also show the thermal conductivity of P-doped Si obtained by using the same experimental setup. Top inset illustrates a schematic of the measurement setup for the thermal conductivity. Bottom inset shows κ divided by T3 as a function of temperature. By plotting the thermal conductivity as a function of T3, one can see how the ballistic regime, where κ ∝ T3, evolves to the Poiseuille regime where κ changes faster than T3. Such a behavior can be seen in BP (B) and in the literature data of Bi (20) (C), but it is absent in P-doped Si (D). (E) Schematic representation of Poiseuille flow: the temperature dependence of mean free path of normal scattering lN (dashed-dotted line) and resistive scattering lR (solid line) and the sample dimension d in the three regimes of thermal transport. Poiseuille flow of phonons takes place in the limited temperature range between ballistic and diffusive transport when the inequality of lNdlR is fulfilled. In the center of samples, the phonon mean free path (dashed line) becomes lphd2/lN in these conditions. The three small panels represent three distinct scattering mechanisms. Normal scattering between two phonons does not lead to any loss of momentum and does not degrade the heat current. In contrast, Umklapp scattering and collision with impurities lead to a loss in heat current and amplify thermal resistance. In any finite-size sample, phonons are also scattered by the boundaries.

In the absence of anisotropy, the lattice thermal conductivity can be represented by the following simple expressionκ=13Cvlph(1)where C is the phonon specific heat, 〈v〉 is the mean phonon velocity, and lph 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 T3 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 T3 just below the peak temperature in BP. This can be seen in Fig. 2B. As a consequence, κ divided by T3 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 T3 behavior of κ (saturation of κ/T3) observed in the low-temperature data of the thickest sample (open circles in Fig. 2A). As a consequence, the effective phonon mean free path, lph, 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 lph using 〈va〉 = 0.536 × 104 m/s, 〈vc〉 = 0.354 × 104 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, lph 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 T3 evolution of κ (Fig. 2C) and a peak in the temperature dependence of lph just below the maximum thermal conductivity (Fig. 3E). Similar features have also been found in crystalline 4He (21), 3He (22), and H (23).

Fig. 3 Thickness dependence of thermal conductivity.

Temperature dependence of the thermal conductivity in samples with different thicknesses along the a axis (A) and along the c axis (B). We also show the data from the thin BP flakes (14). Thermal conductivity shows a thickness dependence in the whole temperature range. The extracted phonon mean free path lph in samples with different thicknesses along the a axis (C) and along the c axis (D). The local maximum and minimum are present in all samples. The insets of the (C) and (D) show plots of lph versus T in a logarithmic scale for BP including the thicker sample and for P-doped Si, respectively. (E) The effective phonon mean free path of Bi for various average diameters from (20).

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 lph 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, lph 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 inequalitylNdlR(2)

Here, d is a characteristic sample dimension, and lN and lR 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 lph ~ d2/lN between successive collisions with the boundaries, and in the ideal case, lph can reach a value larger than the characteristic sample dimension d. If lN increases with decreasing temperature as T−5, then since κ ~ Cvd2/lN, one can expect thermal conductivity to follow as T8 (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 T3 thermal evolution of κ and a nonmonotonic lph(T) (2023), 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, lph is accompanied by a mean-free-path minimum, dubbed Knudsen minimum, which occurs when d/lN ~ 1. This is where the diffuse boundary scattering rate is effectively increased because of normal scattering (Fig. 4A).

Fig. 4 Thickness dependence of the Knudsen minimum.

(A) Schematic illustration of the Knudsen minimum adapted from (32). The diffuse boundary scattering rate is effectively increased by normal scattering when the mean free path of normal scattering lN represented by the dashed circle is comparable with the sample dimension d, producing a local minimum in the effective phonon mean free path. With increasing thickness, the minimum shifts to lower temperatures since the fraction of phonons that suffer numerous normal scattering (phonons in the red region) is larger in the thicker sample. The phonon mean free path normalized by the value at the Knudsen minimum in samples with different thicknesses along the a axis (B) and the c axis (C) and Bi adapted from (20) (D).

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 lph is longer. Moreover, as seen in the case of Bi (Fig. 3E) (20), both the local maximum and the local minimum in lph(T) persist (Fig. 3, C and D). In addition, note that despite the anisotropic thermal conductivity, lph attains values comparable to the sample dimensions at low temperatures, irrespective of the orientation of the heat current, meaning that lph 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 lph 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.

Fig. 5 T−1 dependence of thermal conductivity and relative decrease of thermal conductivity with thickness.

(A) κ as a function of T−1 in samples with different thicknesses. Note the gradual decrease in slope with thickness, as shown in the inset. The temperature dependence of the relative decrease in κ, κratio, when the thickness decreases from 20(30) μm to 6(15) μm along the a(c) axis (B). We defined κratio as (κ(t1) − κ(t2))/κ(t2) × t2/(t1t2), where t1,2 denotes different thicknesses. One expects this to become of the order of unity in the ballistic regime and zero in the high temperature when heat transport is purely diffusive. These features are seen both in Bi (20) (when the average diameter decreases from 0.51 to 0.26 cm) and in Ge (34) (when the width was decreased from 4 to 1 mm). In BP samples of a typical size of 6 to 140 μm, κratio persists (while slowly decreasing) in the diffusive regime.

The absence of superlinear size dependence can be traced back to the expression for phonon kinematic viscosity, proposed by Gurzhi (1), ν = vTlN, which implies that the phonon fluid is non-Newtonian. Here, vT is the phonon thermal velocity whose magnitude is comparable to the velocity of the second sound, so that vT ~ 〈v〉/3 = 2000 to 3000 m/s and lN is estimated to be of the order of 10 μm at the Knudsen minimum. This yields ν ~ 0.02 to 0.03 m2/s, several orders of magnitude larger than the kinematic viscosity of water. This comparison is to be contrasted with the case of electrons in PdCoO2 (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 lph 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 lph normalized by its minimum value lphmin, 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.


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 sp3 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.

Fig. 6 Phonon dispersions.

Calculated phonon dispersions of BP (A) and Si (D). (B) and (E) show the blow-up in the region below 100 cm−1. (C) and (F) show the Brillouin zones along with the high-symmetry points. The out-of-plane direction in BP corresponds to the path along Γ-Y. (G) Calculated PDOS of BP and Si. (H) PDOS below 100 cm−1. The acoustic modes of BP are less dispersive than that of Si, leading to higher PDOS in BP in the low-frequency region. At low temperatures, BP should have much more phase space for the momentum-conserving phonon scattering processes among the linearly dispersive acoustic modes that are relevant for Poiseuille heat flow.

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 T3 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).


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 4He cryostat and RuO2 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 T4 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 material for this article is available at

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.

References (4655)

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.


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.
View Abstract

Stay Connected to Science Advances

Navigate This Article