High frequency atomic tunneling yields ultralow and glass-like thermal conductivity in chalcogenide single crystals

Crystalline solids exhibiting glass-like thermal conductivity have attracted substantial attention both for fundamental interest and applications such as thermoelectrics. In most crystals, the competition of phonon scattering by anharmonic interactions and crystalline imperfections leads to a non-monotonic trend of thermal conductivity with temperature. Defect-free crystals that exhibit the glassy trend of low thermal conductivity with a monotonic increase with temperature are desirable because they are intrinsically thermally insulating while retaining useful properties of perfect crystals. However, this behavior is rare, and its microscopic origin remains unclear. Here, we report the observation of ultralow and glass-like thermal conductivity in a hexagonal perovskite chalcogenide single crystal, BaTiS3, despite its highly symmetric and simple primitive cell. Elastic and inelastic scattering measurements reveal the quantum mechanical origin of this unusual trend. A two-level atomic tunneling system exists in a shallow double-well potential of the Ti atom and is of sufficiently high frequency to scatter heat-carrying phonons up to room temperature. While atomic tunneling has been invoked to explain the low-temperature thermal conductivity of solids for decades, our study establishes the presence of sub-THz frequency tunneling systems even in high-quality, electrically insulating single crystals, leading to anomalous transport properties well above cryogenic temperatures.


I. ADDITIONAL SAMPLE CHARACTERIZATION
A. Rotational x-ray diffraction (XRD) map To identify the crystallographic directions of the crystal, we extended the method of outof-plane Bragg scan in a thin film x-ray diffractometer to confirm the 6-fold rotation axis of the hexagonal structure. The rotational XRD map shown in Fig. S1(c) was constructed by a series of consecutive out-of-plane 2θ/ω scans taken as the crystal plate was rotating around the hypothesized c-axis. The illustrative schematic for the x-ray and crystal orientation projected along a-and c-axis are shown in Fig. S1(a) and (b), respectively. It is essentially a map of reciprocal space that is perpendicular to the rotating axis. The reciprocal lattice points are extended to long streaks along ψ due to the relatively poor divergence of the x-ray 3 beam in this direction in the instrument. All the reflections in Fig. S1(c) are (hk0) type, and we can see a repetition of {100} reflections when rotated 60 degrees away, confirming that the rotating axis is indeed the c-axis. We can also see {110} reflections when the crystal is rotated 30 degrees away from the 100 facet. These observations agree well with the crystal symmetry. The absence of other reflections and the narrow dispersion of observed reflections indicate high crystalline quality.
B. Single crystal x-ray diffraction

D. Electrical resistivity measurement
We measured the temperature dependence of electrical resistivity of a BaTiS 3 (BTS) crystal. Four contacts were made by directly painting silver epoxy onto the crystal surface, as shown in the inset picture in Fig. S4(a). Four-probe resistance was obtained by passing current through two end contacts and measuring voltage across the two inside contacts.
Linear I-V curves were measured in the temperature range shown. The sample dimensions were measured in a microscope and electrical resistivity was calculated. The sample clearly shows semiconducting behavior with increased resistivity when cooled to lower temperatures, as shown in Fig. S4(a). The electronic contribution of thermal conductivity was estimated using the Wiedemann-Franz law. The largest electronic contribution to the total thermal conductivity measured by TDTR was less than 5% at all temperatures, as shown in Fig.S4(b). Hence, the electronic contribution to the thermal conductivity is negligible and the thermal transport is dominated by phonons.

E. Rutherford Backscattering Spectrometry (RBS) for compositional analysis
We use RBS for quantitative determination of the composition of our BaTiS 3 samples.
In RBS analysis, samples are bombarded with He 2+ ions at an energy in MeV range and the energy distribution and yield of the backscattered ions at a given angle are detected.
The energies of backscattered ions depend both on the mass of atoms from which they are scattered as well as the the depth in the sample at which the collisions occur. The number of backscattered ions is directly proportional to the concentration of a given element. For BaTiS 3 , we estimate that the uncertainty of RBS measurement can be as good as 2%. In

B. Steady state temperature rise
Steady state temperature rise (∆T ss ) from pulse laser heating is a major concern when measuring bulk materials with low thermal conductivity using TDTR, which can easily reach tens or even hundreds of K depending on the thermal conductivity of substrate, laser beam spot size and absorbed laser power. As the fitting procedure in TDTR is highly sensitive to the temperature-dependent input thermal properties (volumetric heat capacity, thermal conductivity) of Al transducer and substrate, we normally limit ∆T ss to less than 10 K, or 20 K at most for materials with ultralow thermal conductivity such as BaTiS 3 . Temperature increases larger than this value can lead to inaccurate thermal conductivity measurement.
∆T ss can be estimated using: [1] where P is incident laser power, R is the reflectance of sample surface, w 0 is 1/e 2 radius of laser beam and κ is the thermal conductivity of substrate. This equation is an estimation for the case of quasi-1D heat flow where lateral heat spreading in Al transducer is not important.
However, in our case, as BaTiS 3 has ultralow thermal conductivity, lateral heat spreading in Al transducer is significant thus cannot be ignored.
We thus calculated ∆T ss by setting the heating frequency to zero (no periodic component) in the thermal model for multilayer sample [1]. For a typical experiment, ∆T ss is about half of that estimated using equation 1, indicating the effect of lateral heat spreading in Al transducer. The laser power we used was optimized as a compromise between signal to noise ratio and ∆T ss . We used a total power of 4 mW (pump 2 mW and probe 2 mW), which gave us ∆T ss < 20 K while yielding a good signal to noise ratio of 50. The 1/e 2 beam radii for pump and probe were 10 microns.

C. Heat capacity of BaTiS 3
The heat capacity of BaTiS 3 was determined using first-principles calculations as well as measured using PPMS. We show the results from both measured values and calculated values in Fig. S7. The two show a maximum deviation of 5%. In the uncertainty analysis, we put 5% as the uncertainty of heat capacity.

D. Sensitivity and Uncertainty Analysis
To evaluate the uncertainty of TDTR measurements, we first calculate the sensitivity of all parameters that are sensitive to TDTR measurements. The sensitivity parameter S α is defined as: where R is the absolute value of the ratio of in-phase and out-of-phase of the lock-in amplifier in TDTR measurements and α is the parameter in the thermal model, such as thermal conductivity, heat capacity, or thickness of each layer.
After all the sensitivity parameters are calculated, the measurement uncertainty of thermal conductivity can be calculated using the following equation: where S α is the sensitivity calculated using Eq.2 to parameters used in TDTR modeling, ∆α α is the uncertainty for α, S κ is the sensitivity to thermal conductivity of BaTiS 3 , δφ is the phase uncertainty for lock-in detection,  and 5% for the uncertainty for laser beam spot size. From Equation 3, TDTR measurement uncertainty for the thermal conductivity of BaTiS 3 is about 6% at temperatures from 70 K to 400 K.

E. Transducer thickness
A thin Al transducer is preferred for the current parameters for TDTR measurements, as suggested by the sensitivity plot shown in Fig. S9. In our experiments, we use a thin Al transducer of thickness ∼ 40 nm. The thickness was determined by both picosecond acoustics and AFM.
In picosecond acoustics, the thickness of a film can be determined by measuring the roundtrip time of a thermoelastically generated stress pulse. An echo after the pulse was generated is observed after round trips through the Al film, and thus the Al transducer thickness can be calculated from time interval between echoes t echo and the sound velocity of the film. From Fig. S10(a), we determined that thickness of Al film is 38.6 nm. We also employed AFM to verify the thickness of Al film. As shown in Fig. S10(b), Al thickness was measured to be 38.9 nm, which agrees the result from picosecond acoustics.

III. NEUTRON SCATTERING
The fit to the pair distribution function (PDF) calculated using PDFgui [2] follows a procedure very similar to a refinement of a diffraction pattern and involves starting with a model structure (already known from previous x-ray studies), instrument parameters, and then refining adjustable parameters including the lattice parameter and ADPs to obtain a best fit. The step by step procedure is described in the PDFgui user guide found here https://www.diffpy.org/doc/pdfgui/pdfgui.pdf. The fits were performed for the patterns over ranges of both 10Å and 30Å to ensure consistency, and larger range is used for the final results. All results were consistent. Using the results of these fits the PDF can also be broken down into pair distribution functions for each atomic pair, as shown in Fig. S11.
The partials allow the determination of which atomic pairs contribute to individual features in the measured PDF. However, with increasing pair distances the peaks exhibit more and more overlap making identification by direct inspection more challenging. Also, because Ti has a negative scattering strength it produces negative peaks in the PDF, which when positioned adjacent to positive peaks can produce the appearance of larger positive peaks than if no negative peaks were present. This situation occurs with the Ba-Ba, Ba-S peaks at r = 4.84Å. The apparent height of this peak is partially due to the adjacent negative Ti-Ba and Ti-S peaks. Assuming a single site for the Ti atom the lowest R w value (weighted profile residual) obtained was 0.13. Splitting the Ti atomic position into two sites along the c-axis lowered R w further to 0.126. This amounts to changing the football shaped ellipsoids for Ti in Fig. 2 into more of an hourglass shape. This helps explain why the ADPs from single site fits appear to increase on cooling.  Figure S12 shows the fits to all of the inelastic x-ray scattering data used to determine the phonon dispersion curves of BaTiS 3 . Figure S13 shows elastic x-ray scattering in several regions of the Brillouin zone, demonstrating that diffuse elastic scattering that would arise from static disorder is not observed.