Anisotropic in-plane thermal conductivity of black phosphorus nanoribbons at temperatures higher than 100 K

1Thermal transport in nanoscale materials attracts increasing research attention because of both intriguing phonon physics at the nanoscale as well as growing importance of heat management in nanoscale devices. Especially, layered nanomaterials provide a new platform for fundamental study, such as two-dimensional electronic effects at the mono- or few-layer limit, as well as device applications derived from their high flexibility and van der Waals nature of interlayer coupling. It is well known that thermal conductivity in these layered nanomaterials is much lower along the cross-plane direction than the in-plane direction. It is also theoretically predicted that additional anisotropy may exist in their thermal conductivity along different lattice directions in the basal plane. Recent theoretical studies on few-layer graphene nanoribbons expect anisotropic phonon transport along the zigzag (ZZ) and armchair (AC) lattice directions of its honeycomb structure. This anisotropy is not expected in the bulk, because it arises mainly from different strengths of boundary scattering at the nanoribbon edges with different chiralities, that is, specular scattering at the ZZ edge while angle-dependent scattering at the AC edge, an effect that is reduced with increasing nanoribbon width as the material approaches the bulk. However, experimental demonstration of such in-plane anisotropy of thermal conductivity in thin layered materials is lacking, due mostly to technical challenges in sample preparation and measurements.

2On the materials side, a new member, black phosphorus (BP), recently joins the family of layered nanomaterials as a promising candidate for electronic, optical and optoelectronic applications. BP has been shown to have some remarkable properties compared with other two-dimensional materials, such as high hole mobility (∼1,000 cm2 V−1 s−1 in field-effect transistors) and tunable direct bandgap from 0.3 eV (bulk) to >1.4 eV (monolayer). In addition, BP is considered to be a potentially good thermoelectric material. It has been theoretically predicted that BP has opposite anisotropy in thermal and electrical conductivities: electrical conductivity is higher along the AC direction, while thermal conductivity is higher along the ZZ direction. The anisotropic in-plane thermal transport in BP, in stark contrast to few-layer graphene, is believed to be an intrinsic property, that is, caused by anisotropic phonon dispersion and phonon–phonon scattering rate along the ZZ and the puckered AC directions (Fig. 1a). Therefore, the anisotropy exists not only in the monolayer limit or ribbons but also in multilayer or bulk BP. However, experimental demonstration of this intrinsically anisotropic thermal transport is lacking, owing mostly to the highly preferential growth of single-crystal BP only along the ZZ direction. There exists only one report on the directional thermal conductivity of BP, which was measured at room temperature by Raman thermography.

3In this work, we directly measured the in-plane thermal conductivity of single-crystal BP nanoribbons along the ZZ and AC lattice directions. The measurements were carried out in the condition of steady-state longitudinal heat flow, using suspended-pad micro-devices (Fig. 1b), over a wide temperature range from 30 to 350 K. Our results reveal a high anisotropy in thermal conductivity up to a factor of two at temperatures greater than ∼ 100 K. The high anisotropy is attributed mainly to the anisotropic phonon dispersion, and partially to the phonon–phonon scattering. A size effect in the thermal conductivity was also observed from ∼50- to ∼ 300-nm-thick BP nanoribbons in which thinner nanoribbons show lower thermal conductivity. These discoveries not only shed light on phonon physics in this interesting material but also provide important design guidelines in its device applications.

4BP bulk crystals were synthesized following a well-developed process, as described in Supplementary Note 1. The synthesized bulk BP (Supplementary Fig. 1a, inset) is comprised of bundles of long crystallites. Powder X-ray diffraction taken from exfoliated BP flakes on a slide glass confirms orthorhombic (Cmca) crystal structure (Supplementary Fig. 1b,c) without any other phases. Energy dispersive X-ray spectroscopic elemental analysis (Supplementary Fig. 1d) shows only phosphorous element from the crystals. The X-ray diffraction shows exclusively the planes oriented to the [010] direction, proving the layered nature of the crystals. Exfoliated BP flakes were also analysed using a high-resolution transmission electron microscopy with the [010] zone axis (out-of-plane direction), as shown in Fig. 1c. Each flake has identical lattice structure with identical selected area electron diffraction pattern showing the ZZ and AC planes (Fig. 1d), confirming that each flake is a single crystal. The crystallographic orientations agree well with Raman analysis, which utilizes intensity ratio between the Raman peaks of Ag1 and Ag2 phonon modes (Fig. 1e). It is well known that the intensity ratio of the Ag2 mode to the Ag1 mode is the maximum when the laser is polarized in parallel to the AC direction, because the Ag1 and Ag2 modes correspond to atomic vibration along the [010] (out of plane) and [001] (AC, in-plane) directions, respectively.

5The thermal conductivity of BP was measured by a steady-state longitudinal heat flow method, using suspended-pad micro-devices as described in Supplementary Note 2 with Supplementary Figs 2–3, as well as in previous reports. Each pad has two Pt electrodes for four-probe electrical conduction measurements, and a Pt micro-heater/thermometer to control/sense the temperature for thermal conduction measurements. To make BP nanoribbons with proper geometries to bridge the two suspended pads, a top–down micro-fabrication process was used as shown in Fig. 2a–e. First, crystal orientation of BP flakes was identified using Raman analysis; they were then micro-patterned into nanoribbons with desired widths and lengths using electron beam lithography (EBL) and followed by reactive ion etching (see the Method section). Afterwards, a second EBL was applied to define four openings on each BP nanoribbon, followed by in situ Ar+ milling and Ti/Au deposition with electron beam evaporation. The Ar+ milling is essential to remove the surface oxidized layer, which is well known to quickly form when BP is exposed to ambient air. Very mild milling to remove <1 nm of the BP surface was sufficient to obtain a good thermal/electrical contact between the contact metal and the BP nanoribbons. Such contact-ready BP nanoribbons were manually, individually picked up using a sharp probe tip with the aid of a micro-manipulator, and dry transferred onto a suspended-pad micro-device, with careful alignment of the four Ti/Au metal contacts to the four electrodes on the pads. The aligned nanoribbon was subsequently bonded to the electrodes with Pt deposition using a focused ion beam (FIB). Microscopic images for a target BP flake before and after the micro-patterning, and nanoribbons before and after the metal (Ti/Au and Pt) deposition are shown in Supplementary Fig. 4. We note that direct FIB deposition of Pt on the BP, without the pre-formed Ti/Au metal contacts, failed to yield good electrical contacts. To minimize possible damage by the FIB, the Pt deposition was conducted only for connecting the Ti/Au metal contacts to the device electrodes, rather than covering and bonding the entire nanoribbon onto the electrodes. Between each step of the fabrication processes, the BP samples were stored in dry N2 atmosphere to maximally prevent surface oxidation. Accumulated time period of exposure to ambient air during the entire fabrication process was estimated to be <15 min, which warrants minimum oxidation of the BP surface, as evidenced by the nano-Auger electron spectroscopic surface analysis (Supplementary Fig. 5). All the thermal and electrical measurements were performed in vacuum (<10−6 torr). Figure 2f shows the measured thermal resistance as a function of the length of the BP nanoribbons. The linear fitting and extrapolation yield negligible contact thermal resistance (y intercept) between the nanoribbon and the suspended pads for both the ZZ and AC oriented ribbons. Moreover, the Ti/Au/Pt layers make good Ohmic electrical contacts with the nanoribbons, as evidenced by the linear I–V curves (Fig. 2g).

6We first note that the measured total thermal conductivity (κ) in this study is dominated by the lattice (phonon) contribution. Electronic contribution to the thermal conductivity (κe) is estimated to be negligible (<0.1% of total) in the measured temperature range. This is calculated from the electrical conductivity (σ) measured with four-probe I–V and converted using the Wiedemann–Franz law, κe=L0σT, where L0 is the Sommerfeld value of the Lorenz number (Supplementary Fig. 6). We also note that our device geometry and configuration guarantee that the electrical and thermal currents flow along the same crystal direction in the material along the nanoribbon length, which is important to apply the Wiedemann–Franz law.

7The measured κ along the ZZ and AC crystal orientations is compared with each other in Fig. 3a using ZZ and AC nanoribbons with the same thickness and similar widths (<10% difference). It is evident that the ZZ nanoribbon has a higher κ than the AC nanoribbon, by as much as ∼7 Wm−1 K−1 at temperatures above ∼100 K, while both have similar κ at lower temperatures between 30 and 100 K. The anisotropy in κ, that is, the ratio of κ along the ZZ and AC directions, increases with temperature, reaching up to ∼2 at ∼300 K. The increasing anisotropy ratio is attributed to increased contribution from phonon–phonon (Umklapp) scattering, which becomes the dominating scattering mechanism at high temperatures.

8The temperature dependency is also observed from other BP nanoribbons that have different dimensions, as shown in Fig. 3b,c. To investigate the anisotropic effect of phonon scattering on the thermal transport in the high-temperature regime, we evaluate scattering times using the relaxation time approximation model with phonon dispersions obtained using the force constants calculated from density function perturbation theory, as detailed in the Methods part and Supplementary Note 3. The branch-dependent group velocities are extracted from the phonon dispersions (Supplementary Fig. 7). The sound velocities are summarized in Supplementary Table 1.

9In general, the lattice κ along the temperature gradient direction can be expressed as

10where cq is the specific heat contributed by the mode with the wavevector q, is phonon group velocity, is the unit vector in the direction of temperature gradient, τ is the phonon relaxation time and the summation is over the first Brillouin zone of all q modes.

11Total effects of the phonon–phonon scattering, impurity scattering and boundary scattering are commonly summed up by the Matthiessen's rule,. For the phonon–phonon relaxation time (τu), the common expression

12is used, where B1 and B2 are fitting parameters and ω is phonon frequency. For the impurity scattering, the Klemens' expression is used, where Ai is a parameter that can also be evaluated from the density and species of impurities (generally point defects including native defects such as vacancies and interstitials). The boundary scattering is determined by a characteristic size and surface roughness of the specimen.

13The unknown scattering parameters are usually obtained from fitting κ of the bulk, but the bulk κ of BP specifically along the ZZ or AC direction has not been reported. Therefore, we first obtain these parameters by fitting κ of the nanoribbon with the largest thicknesses in our measurements. This corresponds to 310- and 270-nm thick for the ZZ and AC, respectively. The fitting is validated by the fact that in fitting to the measured κ of thinner nanoribbons, good agreements can be obtained by fixing the Ai, B1 and B2 at the same values and adjusting only the boundary size for the thinner ribbons, as shown in Fig. 3b,c. We obtain the following best-fit parameters: B1,ZZ=5.1 × 10−19 sK−1, B2,ZZ=326 K and Ai,ZZ=1.0 × 10−43 s3 for ZZ nanoribbons, and B1,AC=7.8 × 10−19 sK−1, B2,AC=390 K and Ai,AC=1.0 × 10−43 s3 for AC nanoribbons. Among these, the anharmonicity parameter B1, which is the important parameter at high temperatures, is discussed in the next section. Discussion of the other parameters is presented in Supplementary Note 3.

14The minimum flake thickness suitable for our fabrication process is ∼50 nm due to challenges in dry transferring the ultrathin ribbons. BP nanoribbons thinner than 50 nm were not able to be picked up, due to limitation in the sharpness of probe tip. Figure 4 plots room temperature κ of the BP nanoribbons as a function of ribbon thickness. κ shows a clear thickness dependence along both ZZ and AC directions, decreasing from ∼27 to ∼12 Wm−1 K−1 in ZZ, and from ∼15 to ∼5 Wm−1 K−1 in the AC direction, as the thickness varying from ∼300 to ∼50 nm. The anisotropy of κ is about 2 over this thickness range. The thickness dependence implies that in this thickness range, surface or boundary scattering is effective even at such a high temperature. The κ values measured in this study may differ from predicted κ of few-layer BP, due to different heat transport factors, such as evolution of flexural (ZA, out of plane) phonon mode, or more specular boundary scattering with smooth surface in theoretical simulations of few-layer BP. However, the κ anisotropy, that is, larger κ along the ZZ direction than the AC direction, is in agreement with predictions of monolayer (or few layer) BP. For an example, a recent theoretical work predicts κ ∼110 and ∼36 Wm−1 K−1 at room temperature along the ZZ and the AC direction in monolayer BP.

15Compared with the κ (∼12 Wm−1 K−1 at 293 K) of bulk BP reported by Slack in 1965 (ref. ), our thickest samples have slightly higher values. This is understandable given the fact that their bulk BP was polycrystalline, which can significantly reduce κ by grain boundary scattering, cross-layer scattering and averaging of the anisotropy. The crystallite size of the bulk in that work, obtained from fitting to the temperature dependent κ below 10 K, is only 10% of the visual crystallite size of the crystal. Therefore, it is reasonable that our thick crystal BP nanoribbons have a higher κ. We note that the absolute value of κ measured in our study also can be different from a BP single crystal due to possible damage caused during the device fabrication process.

16At high temperatures, most acoustic phonon states are excited, which leads to the saturation of specific heat at the Dulong–Petit limit. If the phonon–phonon scattering rate is the same (scattering rate difference will be discussed later),. Thus, κ in ZZ direction is higher than AC direction (Fig. 3a) due to the higher group velocities (Supplementary Table 1). If we compare the group velocities of TA1 modes (=2,126 m s−1 and =1,567 m s−1, which are close to the results from ultrasonic experiments in literature), we obtain a κ anisotropic ratio of 1.85, which is close to our experimental result of ∼2 at high temperatures.

17When the temperature drops below 100 K, κ becomes close to each other and has a crossover at ∼40 K. The main physical picture for this crossover is freezing out of phonons when temperature drops. At the low-temperature limit, the excited phonon wavevectors (q) are small, where the Debye approximation is valid such that. Comparing the same phonon branch in the ZZ direction to that in the AC direction, the group velocities are all higher in the ZZ direction. Thus, when counting the phonon states in reciprocal q-space, the higher group velocity would lead to lower density of states in ZZ direction. Thus, κ in the ZZ direction could possibly be lower than κ in the AC direction. Thus, in the intermediate temperature, such as between 30 and 75 K, there is a transition regime that the thermal conductivities are close to each other along these two directions.

18To investigate which phonons contribute the most to the heat conduction, we evaluated the contributions of two TA branches (out-of-plane TA1, and in-plane TA2) and LA branch in the 170-nm-thick nanoribbons. We find that the TA1 mode contributes the most to the heat conduction among the three acoustic modes. It contributes 53% (48%) of the total κ of 19.0 Wm−1 K−1 (13.5 Wm−1 K−1) in the ZZ (AC) nanoribbon (Supplementary Table 2). The TA2 mode and LA mode contribute 24% (19%) and 23% (33%), respectively, to the total κ in the ZZ (AC) direction for these nanoribbons. In addition, κ of TA1 and TA2 branches in the ZZ direction are ∼1.6 and ∼1.8 times of those in the AC direction, respectively, which contributes significant to the anisotropy.

19The phonon–phonon scattering, which is the dominant scattering mechanism at high temperatures, are compared by calculating the τu in ZZ and AC directions using the fit parameters. From equation (2), τu, ZZ is ∼1.25 times of τu, AC at room temperature for the identical phonon frequency. Therefore, the different phonon–phonon scattering also contributes in part to the κ anisotropy. However, compared with the difference caused by the group velocities (that is,, Supplementary Table 1), the difference in τu is small. Therefore, we attribute the high anisotropy in κ mostly to the anisotropic phonon dispersion, and partially to the anisotropic phonon–phonon scattering. In addition, it would also be interesting to know the anharmonic difference, such as the Grüneisen number (γ). In general, the anharmonic parameter B1 is given by (ref. ), where M is the average mass of atoms in the crystal, and θ is the Debye temperature, which is proportional to group velocity. As only γ and are direction-sensitive parameters, the ratio of B1,ZZ and B1,AC is, where and are the group velocities in ZZ and AC directions. Using the extracted group velocities of TA1 modes from the calculated dispersion (the TA1 branch that contributes the most to total κ), we obtain. To the best of our knowledge, this is the first experimental evidence of orientation-dependent anharmonicity in BP.

20Due to the relevance of BP for thermoelectrics, we have also measured electrical conductivity (σ) and Seebeck coefficient (S) along the AC and the ZZ directions, as presented in Supplementary Fig. 8. The room temperature Seebeck coefficients in the AC and ZZ directions are ∼320 and ∼270 μV K−1 (p-type conductivity), and close to that of the bulk (∼330 μV K−1). The power factor (=S2σ) of the AC nanoribbon is three times larger than that of the ZZ nanoribbon. ZT (=S2σT/κ) of the AC nanoribbon (t=170 nm) and the ZZ nanoribbon (t=200 nm) are ∼0.0036 and ∼0.0006, respectively, at room temperature. This low ZT is obviously due to the relatively high κ (single-crystalline samples) and low σ (undoped samples).

21In summary, we discovered a strong in-plane anisotropy in thermal conductivity, up to a factor of two, along the ZZ and AC directions of single-crystal black phosphorous. The BP was prepared in the nanoribbon geometry via a top–down micro-fabrication and the measurements were performed with steady-state longitudinal heat flow using suspended-pad micro-devices. At 300 K, the thermal conductivity decreases as the nanoribbon thickness is reduced from ∼300 to ∼50 nm, but the anisotropy ratio stays around two within the thickness range. Detailed analysis from the relaxation time approximation model shows that the anisotropy in the thermal conductivity originates mostly from phonon dispersion, and partially from phonon–phonon scattering rate, that are both orientation dependent. Our results reveal the intrinsic, orientation-dependent thermal conductivity of BP, which would be important for developing high-performance electronic, optoelectronic and thermoelectric devices using BP, as well as for understanding fundamental physical properties of layered materials approaching the few- or monolayer limit.

22BP bulk crystals were synthesized from red phosphorus powders with SnI4 (American Elements, electronic grade 99.995%) and Sn ingot (Sigma Aldrich) promoters. Please see Supplementary Note 1 for more details. BP flakes were mechanically exfoliated from the synthesized BP bulk crystal, and transferred onto SiO2 substrates, using a polydimethylsiloxane stamp (thickness: 1 cm). The morphology and size, crystal structure, composition, phase and crystal orientation of BP flakes were characterized by optical microscopy, SEM, X-ray diffraction, TEM, Energy dispersive X-ray spectroscopic, Auger electron spectroscopic and Raman spectroscopy (with 488 nm laser). For the TEM analysis, selected flakes were transferred from the SiO2 substrate to TEM grid using a sharp probe tip by the aid of a micromanipulator.

23The fabrication of the BP nanoribbons with Ti/Au metal contacts are well described in the main text (Fig. 2a–e). We use the term nanoribbon based on the geometry of the BP samples. For each EBL, PMMA (C4 950) was spin-coated (4,000 r.p.m., 1 min) onto the BP flakes (on SiO2 substrate), and baked at 120 °C for 5 min. After the first EBL, dry etching to remove the exposed BP stripes was conducted via reactive ion etching (in mixed gas, 90% SF6 and 10% O2). The Ti/Au metal was deposited via in situ Ar+ milling and Ti (10 nm)/Au (70 nm) deposition via electron beam evaporation. Lift-off process was conducted with a gentle shaking in acetone, to avoid fracture of nanoribbons. The final nanoribbons were rinsed thoroughly using isopropyl alcohol.

24Suspended-pad micro-devices were used to simultaneously measure the thermal and electrical transport properties of the nanoribbons. Pt lines were patterned on SiNx based, suspended pads and flexural arms. Each pad has two arms to measure the electrical conductance of the nanoribbon bridging the two pads, and four arms to measure the electrical conductance of pre-patterned Pt micro-heater/thermometers on the pads for thermal conductivity measurements. A prepared nanoribbon was transferred to bridge the suspended pads using a sharp probe tip. Pt/C composite was deposited using FIB to bond the Ti/Au coat of the nanoribbon onto the Pt electrode lines for securing the thermal and electrical contacts. To further secure the contacts, the devices were annealed at ∼370 K for 1 h in vacuum chamber before measurements. To minimize the surface oxidation, all the pristine and fabricated BP samples were stored in a desiccator in dry nitrogen stream with over pressure.

25Thermal (K) and electrical (G) conductance of the nanoribbon were measured simultaneously inside a vacuum chamber (<10−6 torr). Global temperature (Tg) was controlled by an external electrical heater and a cryogenic cooler which were connected to the sample holder. G was measured using the four-probe method, and K was determined following a previous reported method. Briefly, K was obtained using the relation,, where Q is the Joule heating power in the micro-heater, and ΔTh and ΔTc are the temperature change of the hot and cool pad, respectively. To generate Q, a d.c. current was applied to one of the two Pt heater/thermometers. ΔT of each pad was evaluated using the resistance change (ΔR) of each Pt heater/thermometer on it. The length, width and thickness of the nanoribbon in the study were measured using SEM or AFM, to ultimately determine the thermal conductivity. More details can be found in the Supplementary Note 2. We note that the measured voltage (Keithley 2420A) has <1% of error, hence the mean error for each data point in I–V curve is negligible compared with the change in voltage, because current was varied in the I–V measurements. Seebeck coefficient (−ΔV/ΔT) of the nanoribbon was measured by reading the voltage change in the nanoribbon with varying the temperature difference of the two suspended pads.

26To obtain the phonon dispersion of black phosphorous, we first calculate its atomic structure and the Hessian matrix using the Vienna ab initio simulation package. The projector augmented-wave potentials are used for the ion-electron interactions, the optB88-vdW exchange-correlation potential is adopted to include long-range dispersion forces, and the energy cutoff is set to 350 eV. In the density functional theory calculation of the bulk structure, a 10 × 14 × 4 k-mesh is used, and the coordinates are fully relaxed until the forces at each atom are <0.001 eV Å−1. The calculated lattice constants of the conventional cell are a=4.348, b=3.329 and c=10.495 Å, respectively, consistent with experimental values. In the following density functional perturbation theory calculation, the Hessian matrix is calculated with a 3 × 3 × 3 supercell of the primitive cell and a 5 × 4 × 3 k-mesh. On the basis of this Hessian matrix, we extract the force constants and calculate the phonon dispersion using Phonopy.

27How to cite this article: Lee, S. et al. Anisotropic in-plane thermal conductivity of black phosphorus nanoribbons at temperatures higher than 100K. Nat. Commun. 6:8573 doi: 10.1038/ncomms9573 (2015).