Giant suppression of phononic heat transport in a quantum magnet BiCu2PO6

Thermal transport of quantum magnets has elucidated the nature of low energy elementary excitations and complex interplay between those excited states via strong scattering of thermal carriers. BiCu2PO6 is a unique frustrated spin-ladder compound exhibiting highly anisotropic spin excitations that contain both itinerant and localized dispersion characters along the b- and a-axes respectively. Here, we investigate thermal conductivity κ of BiCu2PO6 under high magnetic fields (H) of up to 30 tesla. A dip-feature in κ, located at ~15 K at zero-H along all crystallographic directions, moves gradually toward lower temperature (T) with increasing H, thus resulting in giant suppression by a factor of ~30 near the critical magnetic field of Hc ≅ 23.5 tesla. The giant H- and T-dependent suppression of κ can be explained by the combined result of resonant scattering of phononic heat carriers with magnetic energy levels and increased phonon scattering due to enhanced spin fluctuation at Hc, unequivocally revealing the existence of strong spin-phonon coupling. Moreover, we find an experimental indication that the remaining magnetic heat transport along the b-axis becomes almost gapless at the magnetic quantum critical point realized at Hc.

Thermal transport of quasi low-dimensional (D) quantum spin systems has been widely investigated to uncover both the intriguing quantum mechanical nature of magnetism and the low dimensionality of spin networks. For instance, anisotropic thermal conductivity (κ) of quantum magnets along the chain, ladder leg or in-plane directions clearly reflects the kinetics of magnetic quasiparticles that is highly direction dependent 1,2 . Similar to other thermodynamic properties, κ also reflects low-energy quasiparticle excitations of the quantum magnets. The metallic power-law behavior observed in the temperature (T)-dependent thermal conductivity, κ(T), proximity to zero temperature indicates gapless spinon excitation in 1D spin chains 3 or spinon Fermi surfaces in 2D quantum spin liquids 4 .
The magnetic field (H)-dependent behavior of κ, κ(H), further provides fruitful information on the evolution of spin excitations and their coupling with the lattice degree of freedom. H dependence of κ arises for various reasons, such as extra heat conduction due to magnetic quasiparticles or additional scattering of phononic heat carriers via strong spin-phonon coupling [5][6][7][8][9][10][11] . It is thus expected that the H dependence will be most pronounced near the magnetic quantum critical point (QCP) where the spin gap is closed to excite magnetic quasiparticles easily or spin fluctuation becomes maximized. In this respect, the investigations of magneto-thermal conductivity near the critical magnetic field (H c ) for a quantum phase transition can be a promising way to understand the spin-phonon interaction as well as the magnetism itself. BiCu 2 PO 6 is a new kind of spin-ladder compound that forms a unique spin ground state of a spin dimer crystal or valence bond solid. Its crystal structure inherently contains geometric spin frustration within the ladder (Fig. 1a); it is composed of distorted S = 1/2 (Cu 2+ ions) spin chains with antiferromagnetic (AF) interactions J leg (≈ 138 K) and J NNN (≈ 73 K) along the b axis, thereby having frustrated spin interactions. These chains are coupled to one another along the c axis through another AF interaction J rung (≈ 58 K) to form a two-leg spin ladder in the bc plane [12][13][14] . Various thermodynamic, light scattering and nuclear magnetic resonance measurements have indicated anisotropic spin-gap with the excitation energy (Δ /k B ) ranging from 32 to 45 K, without any long-range spin order down to 0.1 K [12][13][14][15][16][17][18][19][20] . The application of magnetic field then leads to reduction in the spin gap (Fig. 1b) and induces cascades of magnetic phase transitions above critical magnetic fields μ 0 H c = 20-24 T, the values of which depend on the direction of the applied field 21,22 .
Inelastic neutron scattering studies on BiCu 2 PO 6 have also uncovered intriguing spin excitations distinct from other conventional spin ladders. The spin excitation branches,  k ( ) i (i = 0, + 1, − 1), become split even at H = 0 because of giant anisotropic spin interactions originating from strong spin-orbit coupling, which can be described by antisymmetric Dzyaloshinskii-Moriya (D-M) and symmetric anisotropic spin exchange terms. Thus, the ground state can no longer be described as singlet only and the excited states exhibit the characteristics of mixed singlet and triplet states. It has thus resulted that spin excitations become highly dispersive particularly along the b* axis and almost dispersionless along the a* axis 19 . Naturally, the minima of the three excitation branches at H = 0 becomes different, giving rise to the three excitation gaps, Δ 0 /k B ≈ 19 K, Δ +1 /k B ≈ 36 K and Δ −1 /k B ≈ 45 K satisfying Δ 0 < Δ +1 < Δ −1 (Fig. 1b) 20 . In the magnetic Raman spectra, in particular, the polarization dependency in the singlet bound modes has been observed, which constitutes a strong evidence of two-triplon bound states with the binding energy of approximately 33 K 18 . Furthermore, the observation of softening in the singlet modes suggested the melting of a spin dimer crystal (or valence bond solid) through the singlet fluctuations 18 . It has been also suggested that the triplons can even become fractionalized and reconstructed into a soliton lattice at higher magnetic fields above H c 17 .
All those unusual spin excitation spectra and postulated high field phases suggest that study of thermal transport could be useful to understand the evolution of spin and phonon excitations of BiCu 2 PO 6 in the high field region where conventional neutron scattering cannot be performed. In particular, it would be interesting to study how the geometrical frustration and strong D-M interaction can manifest themselves in the magnetic and phononic heat transport across the magnetic QCP. Indeed, a recent thermal transport study of BiCu 2 PO 6 in a low-H region has found a pronounced suppression of κ at ~15 K, resulting in a double-peak feature in κ(T) along all crystallographic directions 23 . This behavior is highly unusual, as the double-peak feature has been found only along the leg direction in most other spin ladder compounds because of the separate contributions of phononic and magnetic heat carriers 1,2 . To shed light on the origin of this anomalous behavior, the application of high H across H c should be also useful because it can provide information on the H-dependent variation of the double-peak feature and the intriguing nature of the strong spin-phonon coupling.
Here, we report an investigation of thermal conductivity in BiCu 2 PO 6 single crystals under high magnetic fields up to 30 T. We found that the dip located between the double peaks in κ(T) moves systematically toward lower temperatures as H is increased to H c , which result in the strongest field-induced suppression of κ known among the magnetic materials (by a factor of ~30). Such large T-and H-dependent variations in κ could be explained in terms of the resonant scattering of the phonons induced by magnetic transitions involving two-triplon bound states 8,18 . By developing a theoretical model for phononic κ having resonant scatterings with the three spin excitation branches, we could extract information on the evolution of those magnetic energy levels (a) The schematic crystal structure of the Cu 2+ S = 1/2 frustrated two-leg ladder in BiCu 2 PO 6 . (b) Expected magnetic energy levels of BiCu 2 PO 6 at H = 0 and H ~ H c . Minima of the excitation branches satisfy Δ 0 ≈ 19 K < Δ +1 ≈ 36 K < Δ −1 ≈ 45 K at H = 0 20 . (c) Temperature-dependent thermal conductivity of BiCu 2 PO 6 single crystals at 0 T and along each direction (κ a , κ b and κ c ). The solid black line represents the calculated phonon component along the a axis (κ a,ph ), and dashed black line represents κ a,ph without resonant scattering processes. The inset shows the measured data divided by the boundary scattering limit of the phonon thermal conductivity (κ/κ bd ).
across H c . Furthermore, κ b (T) measured along the ladder leg direction revealed remnant magnetic contribution that scaled almost linearly in temperature at H c ., supporting that the magnetic heat transport nearly becomes gapless at the magnetic QCP.

Results
Zero-field thermal conductivity of BiCu 2 PO 6 . Figure 1c shows the κ i (T) (i = a, b, c) of the BiCu 2 PO 6 single crystals at μ 0 H = 0 T. Note that above 10 K, κ b becomes clearly higher than κ a and κ c , thereby implying that the heat conduction along the leg (b) direction is higher than that along the interlayer (a) or rung (c) direction. The anisotropy observed in Fig. 1c does not appear to originate from the anisotropic phonon thermal conductivity (κ ph ). To estimate the intrinsic anisotropy in κ ph , we rescaled κ in the low-temperature regime with respect to κ bd , i.e., κ in the 'boundary scattering limit' . From the kinetic theory, the boundary scattering limit of phononic thermal conductivity can be expressed as κ bd = c ph v ph l bd /3, where c ph is a phononic specific heat of a crystal approximately given as βT 3 , v ph is an averaged sound velocity determined from the Debye temperature θ D = v(ħ/k B )(6π 2 n/V) 1/3 = 330 K 24 (n/V: number density of atoms), and l bd is the boundary scattering limit of the phonon mean free path determined from the sample geometry as π = l S 2 / bd (S: cross section of the sample) 25 . The calculated κ/κ bd values under the assumption of isotropic phonon dispersion thus reflect the intrinsic anisotropy in the phonon transport, independent of the crystal geometry. The κ/κ bd curves in the inset of Fig. 1c exhibit anisotropy of less than 20%, which is not sufficient to explain the anisotropy in κ above 10 K, which shows a difference of more than 200%. Different strengths of umklapp scattering can also be ruled out as an explanation because it would lead to different peak positions in κ ph (T), which is inconsistent with the results shown in Fig. 1c 26 . We therefore argue that the κ ph along the a axis (κ a,ph ) should be similar to the κ ph along the b axis (κ b,ph ), i.e., κ a,ph ≈ κ b,ph in BiCu 2 PO 6 , regardless of its crystallographic structure. Thus, the extra heat transport along the leg (b) direction could be magnetic in its origin, which indeed reflects the nature of anisotropic magnetic excitation of BiCu 2 PO 6 .
The most conspicuous feature in the shapes of κ i (T), as shown in Fig. 1c, is the presence of a dip feature sandwiched by the double peak structure, which is observed for every crystallographic direction. Such a double-peak structure has been commonly observed in quasi low-dimensional magnets (spin chains or spin ladders), particularly in only one direction where the magnetic thermal transport reaches a maximum 1,2 . It is therefore anomalous to observe the double peaks along every direction in one of spin ladder compounds, BiCu 2 PO 6 . Moreover, none of the peaks in κ(T) appears to be either a solely phononic or a solely magnetic peak; because the magnetic excitation along the reciprocal a axis is known to be non-dispersive 19 and the velocity of magnetic quasiparticles ( ∇ k ( ) k  ) along the a axis should be much smaller than that along any other direction. Therefore, it is unlikely that magnetic thermal conductivity κ mag can produce a separate peak along the a axis. In addition, attempts to estimate κ ph in every direction from the single peak only below 7 K (or above 20 K) resulted in extremely small κ ph at high (low) temperatures. Therefore, it is necessary to consider another mechanism for explaining the double peak feature in κ ph along every direction.

Giant suppression of thermal conductivity under magnetic field.
To better understand the origin of the double-peak feature, we performed κ(H) measurements along the a and b axes for fields of up to μ 0 H = 30 T (Fig. 2). As H is approached to H c , the dip features in both directions systematically shift toward lower temperatures, resulting in very small κ at low temperatures. Above H c where the spin-ordered phases appear, κ at low temperatures around 1.4 K is again enhanced. Figure 3a shows the magneto-thermal conductivity (κ(H)/κ(0)) at T = 1.4, 3 and 5 K, as summarized from the results in Fig. 2. At 1.4 K, it is clearly observed that κ(H)/κ(0) sharply decreases near H c 22 . However, as temperature is increased, the decreasing region becomes widened and gradually shifts toward higher H. It is noteworthy that κ is suppressed by a factor of ~30 at T = 5 K and μ 0 H = 30 T, which is, to our knowledge, the largest magneto-thermal conductivity observed in magnetic materials. For example, SrCu 2 (BO 3 ) 2 , Ba 3 Mn 2 O 8 , and hexagonal HoMnO 3 showed similar field-induced κ suppression by a factor of ~4 at μ 0 H = 17 T 8 , ~10 at 14 T 6 , and ~9 at ~3 T 9 , respectively. Therefore, the observation of such a large κ suppression under magnetic fields is a unique case. In comparison, such large suppression of κ has been also observed in the hexagonal HoMnO 3 as a function of temperature near the Ho 3+ spin ordering as compared with YMnO 3 . This has been attributed to the strong spin-phonon coupling and magnetic fluctuations arising from the spin frustration 27 .

Resonant scattering of phonon by magnetic excitations in BiCu
To understand the anomalous κ(T) behavior and strong H dependence of BiCu 2 PO 6 in a quantitative level, we consider here the model of the resonant scattering of phonons by magnetic excitations, which has been formerly applied in SrCu 2 (BO 3 ) 2 to explain a similar dip feature 8 . SrCu 2 (BO 3 ) 2 represents an archetypal, exactly-solvable spin-frustrated networks, in which its ground state can be described as a bound two-dimer system having two S = 1/2 spins in each dimer 28 . Similarly, the frustrated spin interactions along the leg direction in BiCu 2 PO 6 may provide a similar binding effect between adjacent dimers formed by the rung coupling as evidenced by a recent Raman scattering study 18 .
In such a bound two-dimer system, a phonon-assisted magnetic transition between the excited states with the same S z value, i.e., one of the thermally excited S = 1 one-triplet, one-singlet states and S = 1 bound triplet states, having energy difference of Δ s can be induced, accompanied by absorbing or emitting a phonon of energy ћω ≈ Δ s due to the spin-phonon coupling. The entire two-step process, magnetic excitation and de-excitation, that restores the initial spin state can then be regarded as ω-dependent elastic scattering of the phonons, and the phonon relaxation rate τ −1 will be directly proportional to the transition rate of the magnetic excitation. In the isolated two-dimer system, there is no correlation between the momentum directions of the phonon at initial and final states. Although the original resonant scattering model 8 assumed a flat dispersion of the magnetic excitation spectra, it is likely applicable to any spin-gapped system with dispersive excitations because most excitations  appear within a narrow k region near the band minima. Furthermore, the evolution of the spin excitation gaps under H is likely to modify the transition rate, resulting in the strength change of the resonant phonon scattering (Fig. 1b).
Thermal transport by phonons or magnetic quasiparticles can be generally described by where c k , v k , l k and n k represent the heat capacity, velocity, mean free path and statistical occupation of quasiparticles of wavevector k. Applying the Debye approximation to Eq. (1) leads to the Callaway model for κ ph , the temperature dependence of which is predominantly determined by the phonon relaxation rate τ −1 τ =∑ − ( ) n n 1 26,29 . To calculate the relaxation rate by the resonant scattering (τ − res 1 ), we considered three nondegenerate excitation branches (S z = 0, + 1, − 1), which are split even at zero field by the anisotropic spin interactions, and are subject to change under H by the Zeeman effect 20 . The phonon scattering rate due to such excited levels can be described as follows (Methods): Here, C i is t he overa l l fac tor rel ate d to t he st rengt h of t he spin-phonon coupling , is the statistical occupation of thermally excited states with , Δ i is the H-dependent energy gap between the ground state ( 0, 0 ) and an excited state ( i 1, ), and Δ s,i is the H-independent resonance energy determined by the energy difference between the magnetic excited states of the two-dimer system. In SrCu 2 (BO 3 ) 2 where degenerate triplet states are realized, where B is the binding energy of the two triplons. In addition, are the probabilities of the excitation and de-excitation processes, respectively.
The overall factors C i are here assumed as magnetic field-independent parameters as long as H does not exceed H c 8 . Although there is a possibility that strength of the elastic spin-phonon coupling decreases as the field-induced spin polarization appears even below H c due to the D-M interactions and paramagnetic contributions 14 , its magnetic moment is sufficiently small (less than 0.005 μ B /Cu 2+ below H c 21 ) and thus the effect could be considered as higher-order terms.
Quantitative analyses of thermal transport with magnetic energy levels. The double-peak feature observed in the κ data (Fig. 1c) can be explained in terms of the resonant scattering model for κ ph described in Eqs (1) and (2), with the numerical parameters listed in Table 1. Because the dispersion of the magnetic excitation is nearly flat along the a axis 19 , we apply the resonant scattering model to κ a . At μ 0 H = 0 T, the dip at ~15 K can be quantitatively explained as a consequence of the resonant scattering in κ ph ; the calculated κ (solid line in Fig. 1c) shows good agreement with the experimental κ a in the entire temperature range. In contrast, the predicted curve for κ ph without considering the resonant scattering (dashed line in Fig. 1c) exhibits a typical single-peak feature Table 1. Numerical parameters used to calculate the phonon thermal conductivity at 0 T. θ D = 330 K was taken from the literature 24 , and the other parameters were adjusted to reproduce the experimental data. The resonance energies for each excitation channel, Δ s,i (i = 0, + 1, − 1), were assumed as an identical value.
expected in nonmagnetic, insulating materials 26 . This observation supports the hypothesis that the resonant scattering in κ ph is the origin of the double peaks and the dip observed in Fig. 1c.
One of intriguing findings obtained from the best fit is that the overall constant related to the strength of the resonant scattering for S z = 0 branch (C 0 ) is at least one order of magnitude smaller than the constant for S z = 1 (C 1 ) or S z = − 1 (C −1 ) branches (Table 1). This indicates that the resonant scattering from the transition between S = 0 one-triplet, one-singlet and S = 0 bound triplet could be negligible as compared with other magnetic excitations, which may constitute a clue for developing further microscopic theories on the resonant scattering.
To understand the κ data under magnetic fields (Fig. 2), we apply the same resonant scattering model, assuming that the energy gaps are dependent on H (see Fig. 1b and Table 2). In a recent inelastic neutron scattering study 20 , each magnetic level exhibited the Zeeman splitting up to μ 0 H = 11 T such that the excitation gaps could be roughly described as follows:  Fig. 2a) are in good agreement with the experimental data, particularly for temperatures above 5 K and fields of up to 20 T. This finding clearly indicates that the double-peak feature in κ a (T) and the suppression of the peak at low temperatures can be mostly explained by κ ph based on the resonant scattering model.
The filled squares in Fig. 3b summarize the H dependence of the three energy gaps, Δ i (H), for fields of up to 20 T, as obtained from the best fits to the results presented in Fig. 2a and Table 2. For example, the energy gaps of the three excitation branches at μ 0 H = 0 T were found to be Δ 0 /k B ≈ 25 K, Δ +1 /k B ≈ 42 K and Δ −1 /k B ≈ 57 K. All the H-dependent energy gaps Δ i (H) obtained by fitting κ a (T) (filled squares) are qualitatively consistent with the recent neutron scattering results (open symbols), although they are larger by ~25% than those obtained from the minimum energies of the magnetic excitation spectra 20 . Because magnetic excitations involving phonon scattering are not confined only to the minima of these dispersive branches, it seems reasonable that the effective gap values estimated from the κ ph fits are higher.
Based on the three excitation gaps obtained from our fit (see Table 2), the averaged gap value becomes Δ~41 K. As the resonance energy (Δ s ) obtained from our fit was approximately 32 K, the triplon binding energy B is then roughly estimated to be B = Δ − Δ s ≈ 9 K. Recent Raman scattering result 18 pointed that the binding energy of the bound singlet is 2B ≈ 33 K and thus B ≈ 16 K, which is rather close to the B from our fits although it is a rough estimation based on the averaged spin gap value. Such agreement supports the presence of the triplon bound state in the resonant phonon scattering processes, and justifies the application of the above model into the BiCu 2 PO 6 system.
Along the b axis, however, the experimental κ b data cannot be explained solely by the resonant scattering model in κ ph . Figure 1c shows that the experimental κ b data are overall larger than κ a,ph (≈ κ b,ph ). To better fit the experimental κ b data, an additional contribution from κ mag was indeed essential. This seems consistent with the fact that the inelastic neutron scattering studies found most dispersive spin excitation branches along the reciprocal b* axis 19 . To calculate κ mag explicitly, we have extended Eq. (1) to derive new formula for κ mag when the excited triplet states are nondegenerate (Methods). The result is Here, β = 1/k B T, and l mag,i is the magnetic mean free path of the i-th branch. The solid lines in Fig. 2b show the calculated κ b ( = κ b,ph + κ b,mag ) results for fields of up to μ 0 H = 20 T, in which κ b,ph is assumed to be the same as κ a,ph . As a result of including the additional κ mag , we could achieve good agreement between the experimental κ b data and the calculated ones, particularly at temperatures above 5 K and fields of up to μ 0 H = 20 T. The excitation gap values obtained from the fit are also plotted in Fig. 3b as filled circles, which overall show good agreement with the excitation gaps obtained from the fit along the a axis. However, when 20 T < H < H c , the theoretically predicted values of both κ a (= κ a,ph ) and κ b (= κ b,ph + κ b,mag ) clearly overestimated the experimental data in a wide temperature window below 10 K. For example, the experimental data at μ 0 H = 22 T are far more suppressed below ~10 K than the theoretically predicted κ a and κ b , even after applying a minimum Δ +1 (H) ≈ 0, which yields the most significant suppression of κ ph in the resonant scattering model (dashed lines in Fig. 2). This observation directly points out that the additional scattering of heat carriers, beyond the resonant scattering model adopted here, should exist particularly near H c . Two possible   origins could be postulated. First, when the spin energy gap becomes sufficiently small, the resonant scattering model, which assumes flat dispersion, may not fully account for the effect of dispersion in the magnetic energy levels. However, it is unlikely that the dispersion effect will bring up enhanced resonant scattering at H ≈ H c . More importantly, beyond the resonant scattering, the additional nonresonant scattering of phononic heat carriers due to the spin fluctuation seems to play a role because the spin fluctuation most likely increases at the magnetic quantum critical point.

Field dependence of magnetic thermal conductivity near the critical magnetic field. The very
large suppression of κ ph at H ≤ H c in both directions suggests the possibility that κ b,mag can be extracted approximately at low temperatures and at H ≈ H c where κ ph becomes sufficiently small due to the enhanced phonon scattering 30 . Assuming that the relation κ a,ph ≈ κ b,ph is valid at T > 7 K or H ≈ H c , we subtracted the κ a from the κ b data to plot (κ b − κ a )/T, as shown in Fig. 4. The (κ b − κ a )/T at T > 7 K, μ 0 H = 0 T shows a clear peak at T~13 K, which is a typical feature of κ mag in the gapped spin ladder systems 1,30 . The κ mag peak is almost insensitive to the magnetic field, suggesting that Δ 0 is indeed magnetic field-independent and magnetic mean free path is not strongly subject to the resonant phonon scattering. Interestingly, the most constant behavior is observed at μ 0 H c = 23.5 T in the low-temperature limit down to 0.8 K, whereas the neighboring data at 23 T or 24 T show clear decreases. This observation implies that the magnetic thermal conductivity becomes nearly linear in the low-temperature limit at the critical field μ 0 H c = 23.5 T.
However, it was difficult to analyze the experimental data at μ 0 H ≤ 22 T and ~1 K < T < 7 K in a more quantitative level because anomalous peak-or dip-like structures were superimposed on top of the exponential tail of κ mag (e.g. gray solid line in Fig. 4). Although the anisotropy ratio between κ a,ph and κ b,ph (= κ a,ph /κ b,ph ) is expected to be less than 20%, the anisotropy of κ ph that is presumably related to the sample geometry and/or directional dependence of the spin-phonon coupling could be still non-negligible to produce such additional peak-or dip-like structures in the (κ b − κ a )/T data.
With the caveat of having such additional structures, to understand whether κ mag can still explain the low temperature experimental data of the (κ b − κ a )/T near H ≈ H c , we calculated κ mag,i /T (i = 0, + 1, − 1) individually with fitted parameters listed in Tables 1 and 3 and plotted κ mag /T ≡ ∑ i = 0,±1 κ mag,i /T, κ mag,0 /T and κ mag,+1 /T in Fig. 4. Here, the energy gaps were assumed to be Δ 0 /k B ≈ 25 K, Δ +1 /k B ≈ 0 K and Δ −1 /k B ≈ 99 K. The calculated κ mag /T results well reproduced the overall shape of (κ b − κ a )/T at μ 0 H = 0 T and 23.5 T, including the magnetic field-independent broad peak at 11 K~13 K. [The hump at ~3 K may be attributable to the anisotropic phonon component that does not exactly cancel out in the (κ b − κ a )/T plot.] Moreover, the results for κ mag,+1 /T with Δ +1 /k B ≈ 0 K (at H = H c ) can explain the constant behavior observed in the low-temperature limit, whereas the κ mag,−1 /T data contribute only very little. This finding indicates that a gapped magnetic excitation to the S z = 0 branch and a nearly-gapless magnetic excitation to S z = + 1 coexist at H = H c . When a finite energy gap of Δ +1 /k B = 1 K or 2 K is used, the calculated κ mag /T results decrease at low temperatures, showing a clear deviation from the constant κ mag /T behavior. Therefore, the energy gap Δ +1 /k B at μ 0 H c = 23.5 T, if any, could be less than 1 K to reproduce constant κ mag /T behavior for temperatures down to 0.8 K. Moreover, the energy gap Δ +1 /k B at 23 T is expected to be approximately 2 K. Based on the experimental data in this work, we estimate that the excitation gaps are close to Δ 0 /k B ≈ 25 K, Δ +1 /k B ≈ 0 K and Δ −1 /k B ≈ 99 K at H = H c . Triangles in Fig. 3b and Table 2 summarize those estimated gap values around H c . The results in Fig. 3b exhibiting the continuous evolution of the three excited levels corroborates that a magnetic quantum phase transition indeed occurs at H = H c due to the crossing of the energy levels between the S z = + 1 excited state and the magnetic ground state. Almost negligible energy gap in the S z = + 1 branch around H c expected from the T dependence of the magnetic thermal conductivity is rather peculiar, because the spin gap in the spin-dimer systems with strong D-M interaction is not expected to be closed completely even at H c due to the mixing of singlet and triplet states. It may indicate that opening of the spin gap by the anisotropic D-M interaction is comparable or less than the energy scale of ~0.8 K 17 , or the D-M interaction is strong enough so that it cannot be treated as a higher-order effect and different physics might be need to understand our observation. Investigations of the critical behavior of BiCu 2 PO 6 at lower temperatures and related theoretical works will be necessary for understanding its magnetism around the critical field.

Discussion
One fundamentally interesting question is why the phonon transport in BiCu 2 PO 6 is subject to strong resonant and nonresonant scatterings with the magnetic energy levels and then result in giant magneto-thermal conductivity. A common feature in SrCu 2 (BO 3 ) 2 and BiCu 2 PO 6 where the resonant scattering process becomes conspicuous, is the existence of spin frustration and anisotropic spin interaction described by the Dzyaloshinskii-Moriya interaction 20,31 . Hence, the spin-orbit coupling represented by the Dzyaloshinskii-Moriya interaction in the frustrated spin network might be one of necessary conditions for having the dominant resonant phonon scattering. Moreover, the conventional phonon scattering due to spin fluctuation effect, beyond the level described by the resonant scattering model, seems to be also significant particularly near H c and at low temperatures in BiCu 2 PO 6 . It gives rise to the additional suppression of κ by a factor of ~3 while the resonant scattering leads to the κ suppression typically by a factor of ~10 according to our fit results. Similar kinds of large κ suppression have been observed in multiferroics with strong magnetoelectric (or spin-phonon) coupling and attributed to the enhanced spin-phonon scattering near the magnetic phase transitions 9-11 . Related to this, it is noteworthy that the spin dimer (valence bond solid) state in BiCu 2 PO 6 lies in the distorted, non-coplanar network of CuO 4 squares that share their edges and corners along the rung and leg directions of the ladders, respectively. Hence, both spin couplings along the rung and leg directions should be greatly affected by the Cu-O-Cu superexchange interaction, which is inherently sensitive to the local Cu-O-Cu boding angle. This structural characteristic seems to provide one natural route to favor strong spin-phonon coupling in BiCu 2 PO 6 . It is further noteworthy that the field-induced electric polarization and magnetoelectric effect have been recently found in the magnetically-ordered, high-field phases of BiCu 2 PO 6 , clearly supporting the existence of strong spin-phonon coupling (B.-G. Jeon et al., unpublished data). Therefore, the presence of both resonant and nonresonant phonon-scattering processes via the strong spin-orbit and spin-phonon coupling in frustrated spin networks seems to be the salient feature for causing the giant suppression of κ near the magnetic quantum phase transition in BiCu 2 PO 6 .
In summary, we have observed the unprecedentedly large suppression of the phononic thermal conductivity in the vicinity of magnetic quantum critical point in BiCu 2 PO 6 single crystals. In contrast to other quasi-low-dimensional magnets, in which both phononic and magnetic heat conduction exist almost independently, both resonant scattering of phononic heat carriers with the magnetic excitation levels and the additional phonon scattering due to the magnetic fluctuation were found to be dominant sources of the κ suppression. We claim the importance of spin-orbit and spin-phonon coupling in the frustrated spin network to boost up such resonant and nonresonant scatterings of phononic heat carriers with the spin degree of freedom. Finally, our work also points to an interesting possibility of using the quantum magnet as a field-tunable thermal insulator under high magnetic field environment while superconducting materials have often been used as a thermal insulator at low field regions below their upper critical fields.

Methods
Sample preparation. High-quality single crystalline BiCu 2 PO 6 samples were grown using the floating zone method. After the crystallographic axes had been aligned using a Laue camera, the single crystals were carefully cut into three small pieces with typical dimensions of 1.2 × 0.2 × 0.2 mm 3 to have the longest directions parallel to L sd (10 −7 m) 1.0 ± 0.1 Table 3. Parameters for calculating κ b,mag − T at 0 T. We assumed that ⁎ T 1 = 70 K (marked as †) based on the energy scales of ⁎ T 0 and − ⁎ T 1 , since κ is nearly insensitive with respect to the ⁎ T 1 and thus it was difficult to determine its value from the fitting.
Scientific RepoRts | 6:36970 | DOI: 10.1038/srep36970 a, b, and c axes. The sample size was kept nearly constant to ensure that the boundary scattering effect was similar for all samples. We note that the chemical formula for this compound should be BiCu 2 (PO 4 )O 2 , in which (PO 4 ) 3− acts as a polyatomic ion group, while the notation in this manuscript was used to follow the former literature.
Thermal conductivity measurements in high magnetic fields. A one-heater, two-thermometer configuration was used to measure κ in a Physical Property Measurement System (PPMS TM ) for fields of up to 14 T. The calibration curves of the Cernox TM thermometers at each magnetic field were measured in prior to the κ(T) measurement. A similar configuration with custom-made, vacuum-sealed 3 He probe was inserted in a resistive magnet for fields of up to 35 T at the National High Magnetic Field Laboratory (NHMFL) in Tallahassee, USA. The magnetic field dependence of the Cernox TM thermometers used in here has been calibrated at fixed temperatures down to 0.6 K and up to 50 T, by using a 65 T short-pulse magnet at NHMFL in Los Alamos, USA. The resistance change of the resistance thermometers due to the temperature gradient inside the sample was measured with the lock-in amplifiers (Stanford Research SR830), while applying the heat with the dc current source (Keithley 6221) and measuring the voltage difference of the heater with the nanovoltmeter (Keithley 2182). A high vacuum of well below 10 −5 mbar was maintained throughout the measurements. For these measurements of κ under high magnetic fields, H was applied along the a axis.
Kinetic theory of phonon thermal conductivity with resonant scattering. The phonon thermal conductivity of an insulating system can be explained with the equation derived from the transport theory 26 with the phonon scattering rate of Calculation based on the concept of the resonant scattering in the isolated S = 1/2 dimer system described above shows that the resonant scattering rate is 8 . Δ s is the resonance energy. When the energy gap between the ground state and three excited states are split, the resonant scattering process for each excited configuration should be considered separately. This can be achieved phenomenologically by decomposing the Eq. (6) into three different terms. Although the Shastry-Sutherland lattice is considered while establishing details of this model 8 , any kind of the spin gap system can be applicable in principle since the model does not require direct information about geometry of the Shastry-Sutherland system. Theory of magnetic thermal conductivity. The calculation for magnetic thermal conductivity also starts from the kinetic theory of the thermal conductivity given in Eq. (1). Now we need to consider the one-dimensional dispersion of the magnetic excitations in BiCu 2 PO 6 . As we mentioned above, the energy gap between ground state ( 0, 0 ) and excited states ( i 1, ) are split due to the anisotropic interactions 20 . Thus, one should rewrite the Eq. (1) for three different branches (ε k,0 , ε k,+1 and ε k,−1 ). Now the formula for the magnetic thermal conductivity can be calculated as follows. For simplicity, one can assume that dispersion of the energy bands ε k,i are same except constant shift (i.e. ε k,1 = ε k,0 + (Δ 1 − Δ 0 ) and ε k,−1 = ε k,0 + (Δ −1 − Δ 0 )). Then n k,i become an explicit function of ε k,i . We also assume one-dimensional heat conduction along the N ladders in a unit area. By considering v k = ħ −1 dε/dk, the summation over k can be replaced to the integration with respect to ε. The integral range for each excitation branch is from Δ i (band minima) to ∆ i max (band maxima). Further approximation such as a constant mean free path l k = l mag can be possible since the quasiparticle excitation for magnetic heat conduction dominates just close to the band minima at k = k 0 1 . By considering all of above, the formula for magnetic thermal conductivity will be