# Lone pair rotation and bond heterogeneity leading to ultralow thermal conductivity in aikinite

Virginia Carnevali, Shriparna Mukherjee, David Voneshen, Krishnendu Maji, Emmanuel Guilmeau, Anthony Powell, Paz Vaqueiro, Marco Fornari

► To cite this version:

Virginia Carnevali, Shriparna Mukherjee, David Voneshen, Krishnendu Maji, Emmanuel Guilmeau, et al.. Lone pair rotation and bond heterogeneity leading to ultralow thermal conductivity in aikinite. Journal of the American Chemical Society, 2023, 145 (16), pp.9313-9325. 10.1021/jacs.3c02536. hal-04237441

HAL Id: hal-04237441

https://hal.science/hal-04237441v1

Submitted on 11 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Lone pair rotation and bond heterogeneity leading to ultralow thermal
conductivity in aikinite

Virginia Carnevali,¹†§ Shriparna Mukherjee,²§ David J. Voneshen,³,⁴ Krishnendu Maji,⁵
Emmanuel Guilmeau,⁵ Anthony V. Powell,² Paz Vaqueiro,²* Marco Fornari¹*

¹ Department of Physics and Science of Advanced Materials Program, Central Michigan
University, Mt. Pleasant, MI 48859, USA

²Department of Chemistry, University of Reading, Whiteknights, Reading, RG6 6DX, United Kingdom.

³ISIS Pulsed Neutron and Muon Source, Rutherford Appleton Laboratory, Chilton, Didcot,
Oxon OX11 0QX, United Kingdom

⁴Department of Physics, Royal Holloway University of London, Egham, TW20 0EX, United Kingdom

⁵CRISMAT, CNRS, Normandie Univ, ENSICAEN, UNICAEN, 14000 Caen, France

†Present address: Laboratory of Computational Chemistry and Biochemistry, École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland

§ These authors contributed equally to this work.

* Corresponding authors: E-mail addresses: [p.vaqueiro@reading.ac.uk](mailto:p.vaqueiro@reading.ac.uk) (Paz Vaqueiro), [marco.fornari@cmich.edu](mailto:marco.fornari@cmich.edu) (Marco Fornari)
**Abstract**

Understanding the relationship between crystal structure, chemical bonding and lattice
dynamics is crucial for the design of materials with low thermal conductivities, which are
essential in fields as diverse as thermoelectrics, thermal barrier coatings, and optoelectronics.
The bismuthinite-aikinite series, Cu<sub>1-x</sub>[xPb<sub>1-x</sub>Bi<sub>1+x</sub>S<sub>3</sub> (0 &le; x &le; 1, where x represents a
vacancy), has recently emerged as a family of n-type semiconductors with exceptionally low
lattice thermal conductivities. We present a detailed investigation of the structure, electronic
properties and the vibrational spectrum of aikinite, CuPbBiS<sub>3</sub> (x = 0), in order to elucidate
the origin of its ultralow thermal conductivity (0.48 W m<sup>-1</sup> K<sup>-1</sup> at 573 K), which is close to the
calculated minimum for amorphous and disordered materials, despite its polycrystalline
nature. Inelastic neutron scattering data reveal an anharmonic optical phonon mode at ca. 30
cm<sup>-1</sup>, attributed mainly to the motion of Pb<sup>2+</sup> cations. Analysis of neutron diffraction data,
together with *ab initio* molecular dynamics simulations, shows that the Pb<sup>2+</sup> lone pairs are
rotating and that, with increasing temperature, Cu<sup>+</sup> and Pb<sup>2+</sup> cations, which are separated at
distances of ca. 3.3 Å, exhibit significantly larger displacements from their equilibrium
positions than Bi<sup>3+</sup> cations. In addition to bond heterogeneity, a temperature-dependent
interaction between Cu<sup>+</sup> and the rotating Pb<sup>2+</sup> lone pair is a key contributor to the scattering
effects that lower the thermal conductivity in aikinite. This work demonstrates that coupling
of rotating lone pairs and vibrational motion is an effective mechanism to achieve ultralow
thermal conductivity in crystalline materials.
INTRODUCTION

Thermal transport is of paramount importance for a broad range of technological applications. For instance, thermally insulating materials can be used to prevent heat from damaging critical components,¹ and the performance, safety and lifetime of batteries in electric vehicles can be improved with optimised thermal management.² Finding materials that combine the desired thermal conductivity with other properties, required for specific technological applications, can be exceptionally challenging. This is the case for thermoelectric devices, which enable conversion of a temperature difference into electrical power. For thermoelectric applications, materials with low thermal conductivity are required in order to limit parasitic heat transfer. However, these materials also need to be excellent electrical conductors; this is a conflicting requirement because heat is transported by electrons as well as by phonons. A variety of extrinsic and intrinsic strategies have been proposed in order to minimize phonon transport (which determines the lattice thermal conductivity, κ<sub>L</sub>) in thermoelectric materials.³ While extrinsic effects such as multiscale hierarchical structuring,⁴ grain-boundary engineering⁵ and nanoprecipitates⁶ are effective at lowering κ<sub>L</sub>, they can also adversely influence other properties, such as mechanical and thermal stability, as well as the charge-carrier mobility. Recent research efforts have focused on intrinsic mechanisms that reduce κ<sub>L</sub> by tuning the structure and bonding of materials. These mechanisms include complex crystal structures with a large number of atoms per unit cell,⁷ order-disorder phenomena,⁸,⁹ liquid-like ionic mobility (phonon-liquid-electron-crystal),¹⁰ rattling,¹¹ resonant bonding,¹² anharmonicity induced by lone pairs,¹³ and bonding heterogeneity.¹⁴,¹⁵

The efficiency of thermoelectric energy recovery is related to the figure of merit, $ZT = \sigma S^2 T / (\kappa_e + \kappa_L)$ (where $\sigma$, $S$, $T$, $\kappa_e$, and $\kappa_L$ are the electrical conductivity, the Seebeck coefficient, the operating temperature, and the electronic and lattice components of the thermal conductivity, respectively) of the thermoelectric materials found in the device.¹⁶ Among the chalcogenides, sulfide-based minerals are attractive as potential thermoelectric materials due to the large terrestrial abundance and availability of sulfur when compared to selenium and tellurium.¹⁷ A feature, common to many of the best performing thermoelectric sulfides, is their low thermal conductivity, the origin of which is not fully understood. Promising p-type sulfide minerals, with a thermoelectric figure of merit, $ZT$, approaching unity at moderate temperatures and low thermal conductivity, include tetrahedrites, $Cu_{12+x}Sb_4S_{13}$¹⁸,¹⁹ and $Cu_{12-x}M_xSb_4S_{13}$ ($M = Zn, Ni$),²⁰,²¹,²² colusites, $Cu_{26}T_2M_6S_{32}$ ($T = Cr$,
Mo, W and M = Ge, Sn)²³,²⁴,²⁵,²⁶ and bornite, Cu₅FeS₄.²⁷ By contrast, progress on the corresponding *n*-type sulfides has been limited,²⁸ with several bismuth-containing sulfides being amongst the best *n*-type candidates for thermoelectric applications at moderate temperatures. For instance, bismuthinite, Bi₂S₃, reaches *ZT* ≈ 0.6 at 773 K when doped with chlorine,²⁹ while CdPb₂Bi₄S₉, which is a member of the pavonite-homologous series, M*n*+1Bi₂Q*n*+5 (*n* = 4), exhibits a figure of merit, *ZT* = 0.53 at 775 K.³⁰

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Comparison of the structures of (a) Bi₂S₃ and (b) CuPbBiS₃ (space group *Pnma*); the black rectangle highlights the Bi₄S₆ ribbons in Bi₂S₃ and the Bi₂Pb₂S₆ ribbons in aikinite.</figcaption></figure>

Very recently, it has been reported that the quaternary sulfide CuPbBi₅S₉ exhibits ultralow
thermal conductivity, κ ≈ 0.5 W m⁻¹ K⁻¹, and can reach ZT = 0.43 at 700 K upon doping.³¹,³²
This material is a member (x = 2/3) of the bismuthinite-aikinite series, Cu₁₋ₓPb₁₋ₓBi₁₊ₓS₃ (0
≤ x ≤ 1); all of which exhibit closely-related crystal structures (Figure 1).³³,³⁴ The structure of
bismuthinite, Bi₂S₃ (x = 1), which is highly anisotropic, contains Bi₄S₆ ribbons arranged in a
herringbone pattern. In aikinite (x = 0), half of the Bi³⁺ cations are replaced with Pb²⁺, with
Cu⁺ cations filling tetrahedral holes between the Bi₂Pb₂S₆ ribbons. In the bismuthinite-
aikinite series, the aikinite end member (x = 0), as well as kruptaite, CuPbBi₃S₆ (x = 0.5),
adopt unit cells based on the prototype illustrated in Figure 1(b). In natural specimens,
minerals with intermediate compositions exhibit superstructures based on the ordered
intergrowth of blocks of bismuthinite, aikinite, and kruptaite.³⁵ However, cation ordering is
extremely slow, and synthetic materials with intermediate compositions adopt a copper-
deficient aikinite structure, which contains disordered vacancies at the tetrahedral copper
site.³⁶
Despite the large number of minerals in the aikinite-bismuthinite series, little is known about the electrical and thermal transport properties of these materials, with exception of those of Bi₂S₃ and CuPbBi₅S₉,²⁹,³¹ which are *n*-type semiconductors. Here, we present a detailed study of the structure and transport properties of aikinite, CuPbBiS₃, from an experimental and theoretical point of view. We demonstrate that CuPbBiS₃ is a crystalline *p*-type semiconductor with ultralow lattice thermal conductivity, close to the minimum lattice thermal conductivity calculated on the basis of Cahill's model³⁷ for amorphous and disordered materials. Although neutron diffraction data indicates that in the crystal structure of aikinite the Cu⁺, Pb²⁺ and Bi³⁺ cations are fully ordered, our sound velocity measurements reveal that, at room temperature, the phonon mean-free path is only *ca.* 5 Å, which is approximately twice the interatomic spacing. The temperature dependence of the vibrational spectrum has been investigated by *ab initio* molecular dynamics (AIMD), while inelastic neutron scattering (INS) data has been exploited to estimate the lifetime of the low-frequency Pb²⁺ mode, which is only *ca.* 0.4 ps. Our analysis of the calculated and experimental vibrational density of states (vDOS), which are in excellent agreement, provides clear evidence for the presence of Cu⁺ rattling-like modes, together with anharmonic low-energy modes arising from weakly bonded Pb²⁺ cations. *Ab initio* molecular dynamics simulations reveal that the intrinsic mechanism responsible for the ultralow thermal conductivity in aikinite entails the cooperative interaction between the rotating lone pair on the Pb²⁺ cations and the Cu⁺ cations.

## EXPERIMENTAL AND COMPUTATIONAL SECTION

**Material synthesis:** CuPbBiS₃ was prepared by mechanical alloying of stoichiometric amounts of Cu (Sigma Aldrich, powder 425 µm, 99.5%), Pb (Goodfellow, rods 3.2 mm dia., 99.95%), Bi (Alfa Aesar, needles, 99.99%), and S (Sigma Aldrich, flakes, 99.99%). The reagents were handled under an Ar atmosphere, inside a glovebox. Prior to loading the elements into a 45 ml ball-milling jar made of stainless steel, the Pb rods were cut into small pieces and Bi needles were ground into a powder using a pestle and mortar. Eighteen stainless-steel balls, each with a diameter of 10 mm and a weight of 4 g, were added to the ball-milling jar. A powder to ball weight ratio of 1:12 was used. Milling was carried out using a Fritsch Pulverisette 6 Planetary Ball Mill at 500 rpm, for 60 hours and stopped for 10 minutes every 10 minutes of milling. Following milling, the resulting powder was sealed into an evacuated (< 10⁻⁴ mbar) fused-silica ampoule. The sealed ampoule was heated to 573 K (at a rate of 1 K min⁻¹), held for 48 hours at this temperature and subsequently cooled to room
temperature (at a rate of 1 K min⁻¹). The annealed powder was hand ground in air and consolidated into a densified pellet of ca. 13 mm diameter, by hot pressing under N₂ at 473 K under 80 MPa for 1 hour. The pressure was released, and the hot press was then cooled down for 1 hour to room temperature. The density of the hot-pressed pellet was determined by the Archimedes' method, using an AE Adam PW 184 balance. The pellet has a density more than 96% that of the crystallographic density of aikinite.

**X-ray and neutron diffraction data collection and analysis:** Powder X-ray diffraction data were collected on a Bruker D8 Advance Powder X-ray diffractometer equipped with a LynxEye detector and operating with monochromatic Cu Kα1 (λ=1.54046 Å) radiation. Data collections of 8 hours over the range 10 ≤ 2θ/° ≤ 120 were used. Lattice parameters were determined by the Rietveld method, carried out using GSAS.³⁸ A shifted Chebyschev polynomial with ten coefficients was used to model the background and a pseudo-Voigt function was used to model the peak shape. High-resolution neutron powder diffraction data were collected on the time-of-flight POWGEN beamline³⁹ at the Spallation Neutron Source (Oak Ridge National Laboratory, US). Data were collected at room temperature for 3 h using the center wavelength setting of 1.5 Å, with a *d*-spacing coverage of 0.5 ≤ *d*/Å ≤ 11.8. The powder sample was loaded into a vanadium can with 8 mm inner diameter and sealed with a copper gasket and aluminum lid. Rietveld refinements using neutron data were carried out using GSAS-II.⁴⁰ A logarithmic interpolation function with twenty terms was used to model the background contribution. Lattice parameters, atomic coordinates, profile parameters and the phase fraction for both aikinite and PbS were refined. A single isotropic atomic displacement parameter (Uiso) was used for all the sulfur atoms in the aikinite phase and refined along with those of copper, lead and bismuth atoms. Different structural models for the aikinite phase were explored as detailed in the results section.

**Optical measurements:** A diffuse reflectance measurement was carried out over the wavelength range 200 ≤ λ/nm ≤ 2500, with a step size of 1 nm, and using an Agilent Cary 7000 spectrophotometer equipped with a diffuse reflectance accessory. The reflectance data were transformed to the corresponding absorption spectra using the Kubelka-Munk function,

$F(R_∞) = \frac{K}{S} = \frac{(1-R_∞)^2}{2R_∞}$ where R_∞ is the reflectance of an infinitely thick specimen, K and S are the absorption and scattering coefficients respectively. The band gap (E_g) was estimated, using the Tauc method,⁴¹ from the linear fit of $(F(R_∞)hv)^{1/n} = A(hv - E_g)$ vs. hv, where A is a
proportionality constant, *h* is Planck’s constant, *v* is the photon frequency, and *n* = 1/2 or 2 for direct and indirect transitions respectively.

**Thermal and electrical property measurements:** Differential scanning calorimetry (DSC) data over the temperature range $300 \le T/K \le 575$ with a heating rate of 10 K min$^{-1}$ were collected under a flowing $N_2$ atmosphere using a TA-Q2000 DSC instrument. Thermogravimetric analysis (TGA) was carried out over the temperature range $300 \le T/K \le 970$ using a TA-TGA Q50 under a $N_2$ atmosphere; a heating rate of 10 K min$^{-1}$ was used. The electrical conductivity and Seebeck coefficient were measured simultaneously using a Linseis LSR 3 instrument, using a 4-probe configuration, under a helium atmosphere. Measurements were carried out on pellets with a diameter of 12.7 mm and thickness of ~1.5-2 mm. The electrical conductivity and Seebeck coefficient were measured over the temperature range $423 \le T/K \le 573$. A current of 20 mA was used for the conductivity measurements and a maximum gradient of 50 K was maintained between the upper and lower electrodes for the measurement of the Seebeck coefficient. The instrument was calibrated using a constantan reference. Thermal diffusivity (*D*) measurements were made on graphite-coated circular pellets with a diameter of 12.7 mm and thickness of ~1.5-2 mm, using a Netzsch LFA 447 NanoFlash system, over the temperature range $273 \le T/K \le 573$. Data were analyzed using Cowan’s model with a pulse correction applied. The thermal conductivity (κ) was then calculated from the relation, $κ= DdC_p$ where *d* is the density of the material and $C_p$ is the specific heat capacity. The Dulong-Petit limit for $C_p$, which for $CuPbBiS_3$ is 0.259 J g$^{-1}$ K$^{-1}$, has been used. The uncertainties in the values of the electrical resistivity, Seebeck coefficient and thermal conductivity are 5%, 5% and 10% respectively. The calculations of the minimum thermal conductivity are described in the SI. The Lorenz number L was determined using the relation $L= 1.5+\exp(-|S|/116)$,⁴² where *S* is the temperature dependent Seebeck coefficient. Using the Wiedemann-Franz relation, the electronic part of thermal conductivity $κ_e$ was calculated. Hall effect measurements, to determine the charge carrier concentration, were carried out at room temperature, using a Physical Properties Measurement System (PPMS, Quantum Design) under applied magnetic fields of up to 9 T.

**Sound velocity measurements:** The transverse ($v_T$) and longitudinal velocities ($v_L$) of sound were measured on a disc-shaped sample of ca. 12.5 mm diameter and ca. 2 mm thickness using an Olympus ultrasonic flaw detector (model 38DL plus) with a transducer frequency of 5 MHz.
**Inelastic neutron scattering data:** INS data were collected using the LET spectrometer (ISIS Neutron and Muon Source, UK).<sup>43</sup> The powder sample was loaded into an annular aluminum can. Data were collected at four temperatures, 10, 100, 200 and 300 K. The LET choppers were set up to use incident energies of 24.93, 9.01, 4.60 and 2.79 meV. Identical measurements were carried out for the empty aluminum can and the instrumental background was subtracted using the Mantid package.<sup>44,45</sup> The INS data were integrated from 3 to 5 Å<sup>-1</sup> in Q-space. The neutron-weighted phonon density of states was normalized using a custom Python script. The phonon energies were analyzed using the DAVE package.<sup>46</sup> Peaks were modelled using individual Gaussian functions.

**Electronic bands and transport:** Ab initio electronic calculations were performed using Quantum ESPRESSO<sup>47</sup> as incorporated in the high-throughput infrastructure AFLOWπ.<sup>48</sup> We used ultrasoft PBE pseudopotentials, well converged basis sets corresponding to an energy cutoff of 60 Ry for the wave functions and 600 Ry for the charge density, and the ACBN0 approach<sup>49</sup> to self-consistently determine the values for the Hubbard corrections for each atomic species of the material (U(Pb) = 0.003 eV, U(Cu) = 3.403 eV, U(Bi) = 0.014 eV, and U(S) = 1.589 eV). Hubbard corrections were applied to the 3d orbitals of Cu, 6p orbitals of Pb and Bi, and 3p orbitals of S. Spin-orbit coupling was included in the calculation. To integrate over the Brillouin zone, we used a 4 × 8 × 4 Monkhorst-Pack k-point mesh.<sup>50</sup> The optimized theoretical lattice parameters used for the ab initio simulations are *a* = 11.943 Å, *b* = 4.058 Å and *c* = 11.321 Å. The effective masses have been computed with the method developed by Supka *et al.*<sup>51</sup>

**Born-Oppenheimer ab initio molecular dynamics:** AIMD simulations with the mixed Gaussian and plane wave (GPW) method as implemented in the CP2K package<sup>52</sup> were performed to compute the vibrational and structural properties of aikinite as function of temperature. Valence electrons were expanded as a double-ζ Gaussian basis set with polarization functions (DZVP).<sup>53</sup> The energy cutoff for the electron density expansion in the GPW method was 400 Ry. The temperature was controlled by the velocity rescaling thermostat of Bussi *et al.*<sup>54</sup> with a time constant of 1.0 fs. The system was first equilibrated to 300 K for 10 ps in the isothermal-isobaric ensemble (NPT) with P = 1 atm. The system was then equilibrated in the microcanonical (NVE) ensemble and statistics were gathered for the last 10 ps of the production run. Maximally localized Wannier functions (MLWF)<sup>55</sup> and their centers (MLWFC) were obtained using CP2K, minimizing the MLWF spreads as
explained by Beghold et al.⁵⁶ The Pb²⁺ lone pair dynamics has been characterized through the rotational time correlation function (TCF)⁵⁷

$$C_{rot}(t) = \langle P_2(\mathbf{dip}(t) \cdot \mathbf{dip}(0)) \rangle$$

where **dip**(t) is the Pb²⁺ MLWFC dipole moment at time t and P₂(x) is the second order Legendre polynomial.

**Vibrational properties:** Phonon dispersions and the vDOS for aikinite were reported by Maji et al.²⁶ and serve as a starting point for the computation of the lattice thermal conductivity in the quasi-harmonic approximations.⁵⁸,⁵⁹ The transverse and longitudinal sound velocities have been derived in two different ways: from the phonon dispersion and from the computed elastic constants (see SI). All AIMD simulations were performed at the Γ-point in a 2×5×2 supercell. The lattice parameters of the triclinic simulation box were allowed to relax for 10 ps in the NPT ensemble with P = 1 atm. The system was then equilibrated over 5 ps in the isothermal-isochoric ensemble (NVT) using the supercell volume obtained by the NPT runs; 20 ps of simulation were used to estimate temperature dependent effects at 100 K, 200 K, and 300 K. The vDOS, D(ω), was calculated as a Fourier transform of the velocity autocorrelation function (VACF) as:

$$D(\omega) = \frac{1}{3Nk_B T} \int_0^\infty \frac{\langle \mathbf{v}(0) \cdot \mathbf{v}(t) \rangle}{\langle \mathbf{v}(0) \cdot \mathbf{v}(0) \rangle} e^{i\omega t} dt$$

where ⟨v(0) · v(t)⟩ is the VACF computed over the productions run trajectory; ω is the frequency, N is the number of atoms, k_B is the Boltzmann constant, T is the absolute temperature. We found good agreement between the Quantum ESPRESSO and the CP2K results when comparisons were possible. From the calculated vDOS, we have also computed the Helmholtz free energy F, the internal energy E, the entropy S and the specific heat C_v at zero pressure (SI, Figure S1).

## RESULTS AND DISCUSSION

**Crystal structure and bonding:** Powder X-ray diffraction data collected at room temperature can be indexed on an orthorhombic unit cell (space group *Pnma*); Rietveld refinement yielded lattice parameters of a = 11.6136(2) Å, b = 4.0433(1) Å and c = 11.3675(2) Å, which are in good agreement with those previously reported for aikinite (CuPbBiS₃).⁶⁰ This refinement also indicates that the sample is essentially a single phase (ca.
99 wt.%), with only trace amounts of PbS (ca. 1 wt.%) present. Changes to the synthetic procedures, including sealed-tube synthesis instead of ball milling, different hot-pressing temperatures, and replacement of elemental Pb with PbS, were attempted to remove the traces of PbS. These attempts did not increase the weight percentage of aikinite above 99%.

According to the X-ray structural model for aikinite,⁶¹ in which differences in bond lengths were used to allocate crystallographic sites to bismuth and lead, the Pb²⁺ and Bi³⁺ cations are fully ordered into two distinct crystallographic sites, M1 (0.332, ¼, 0.448) and M2 (0.0185, ¼, 0.681), respectively. As it is impossible to distinguish between isoelectronic Pb²⁺ and Bi³⁺ by X-ray diffraction methods, we collected room-temperature powder neutron diffraction data, which provides contrast between this pair of elements (b<sub>Pb</sub> = 9.4; b<sub>Bi</sub> = 8.5 fm). Rietveld refinements using neutron data were carried out considering three alternative scenarios for the cation distribution: (1) all Pb²⁺ cations located on M1 and all Bi³⁺ cations on M2, (2) all Pb²⁺ cations located on M2 and all Bi³⁺ cations on M1, and (3) Pb²⁺ and Bi³⁺ cations disordered between the M1 and M2 sites.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Rietveld refinement using neutron diffraction data for CuPbBiS<sub>3</sub>; (a) data over the range 0.9 ≤ *d* / Å ≤ 3.5 and (b) 3 ≤ *d* / Å ≤ 8.5. The refinement corresponds to structural model (1) in which all Pb²⁺ cations are located on the M1 site and all Bi³⁺ cations on the M2 site.</figcaption></figure>

<table><caption>**Table 1.** Refined parameters for CuPbBiS<sub>3</sub> (space group *Pnma*) obtained from a Rietveld refinement using neutron diffraction data collected at room temperature.</caption><thead><tr><th>Atom</th><th>Wyckoff site</th><th>x</th><th>y</th><th>z</th><th>SOF</th><th>U<sub>iso</sub> (Å<sup>2</sup>)</th></tr></thead><tbody><tr><td>Pb (M1)</td><td>4c</td><td>0.33288(7)</td><td>0.25</td><td>0.48930(7)</td><td>1</td><td>2.29(2)</td></tr><tr><td>Bi (M2)</td><td>4c</td><td>0.01707(6)</td><td>0.25</td><td>0.68216(7)</td><td>1</td><td>1.43(2)</td></tr><tr><td>Cu</td><td>4c</td><td>0.23464(8)</td><td>0.25</td><td>0.20908(8)</td><td>1</td><td>2.10(2)</td></tr><tr><td>S1</td><td>4c</td><td>0.0470(2)</td><td>0.25</td><td>0.1391(2)</td><td>1</td><td>1.28(3)</td></tr></tbody></table>
<table><tr><td>S2</td><td>4c</td><td>0.3790(2)</td><td>0.25</td><td>0.0558(2)</td><td>1</td><td>1.28(3)</td></tr><tr><td>S3</td><td>4c</td><td>0.2143(2)</td><td>0.25</td><td>0.8013(2)</td><td>1</td><td>1.28(3)</td></tr></table>

R<sub>wp</sub>= 2.87%, GOF= 2.39

Lattice parameters: a = 11.61722(7) Å, b = 4.044072(23) Å, c = 11.37165(7) Å

Weight fraction: CuPbBiS<sub>3</sub>= 98.99(5)%; PbS= 1.01(5)%

<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Distorted octahedral coordination of Bi<sup>3+</sup> (purple sphere) with sulfur (yellow sphere). (b) One-dimensional [BiS<sub>4</sub>]<sup>-</sup> chains directed along the *b*-axis. (c) Bi-S radial pair distribution function (RDF) computed from the NVT *ab initio* molecular dynamics trajectory. (d) Distorted tetrahedral coordination of Cu<sup>+</sup> (blue sphere) with sulfur. (e) Chains of corner-sharing [CuS<sub>4</sub>]<sup>7-</sup> tetrahedra along the *b*-axis. (f) Cu-S RDF computed from the NVT *ab initio* molecular dynamics trajectory. (g) Capped octahedral coordination of Pb<sup>2+</sup> (grey sphere) with sulfur. (h) Ribbons, of stoichiometry [PbS<sub>3</sub>]<sup>-</sup>, directed along the *b*-axis. (i) Pb-S RDF computed from the NVT *ab initio* molecular dynamics trajectory. Legend for the RDF plots: blue, orange and green lines correspond to simulations at 100, 200 and 300 K respectively; the red shaded area highlights the experimental distances.</figcaption></figure>

Clear mismatches in intensities are observed when the Pb<sup>2+</sup> cations are located on the M2 site and the Bi<sup>3+</sup> cations are located at the M1 site (SI, Figure S2 (a) and (b)), leading to a significantly higher R<sub>wp</sub> factor. Structural disorder associated with the M1 and M2 sites
(model (3)) was introduced by refining the site occupancy factors (SOFs), with the constraint
that each site remained fully occupied and that the overall stoichiometry was maintained.
This, however, did not improve R<sub>wp</sub> (SI, Figure S2 (c) and (d)). The lowest value of R<sub>wp</sub> is
found when Pb<sup>2+</sup> cations and Bi<sup>3+</sup> cations are fully ordered on M1 and M2 sites respectively
(Figure 2). Refined parameters are presented in Table 1, while selected bond lengths and the
corresponding bond valence sums are presented in SI, Tables S1-S4. The latter are consistent
with the formal oxidation states of Pb<sup>2+</sup>, Bi<sup>3+</sup> and Cu<sup>+</sup>.

Neutron diffraction confirms that the structure of aikinite (Figure 1) contains three
crystallographically distinct cation sites, which are occupied, in a fully ordered fashion, by
Bi<sup>3+</sup>, Pb<sup>2+</sup> and Cu<sup>+</sup> cations. The Bi<sup>3+</sup> cation adopts a highly distorted octahedral coordination
(Figure 3(a)), and forms one-dimensional [BiS<sub>4</sub>]<sup>-</sup> chains of edge-sharing of [BiS<sub>6</sub>]<sup>3-</sup>
octahedra, oriented parallel to the *b*-axis (Figure 3(b)). The highly distorted octahedral
coordination of the [BiS<sub>6</sub>]<sup>3-</sup> octahedra, with three shorter and three longer Bi-S distances,
arises from the displacement of the central Bi<sup>3+</sup> cation towards one of the octahedral faces.
The marked polyhedral distortion is reflected in large values of the bond-angle variance<sup>62</sup> as
well as non-zero values of bond-length distortion<sup>62</sup> (SI, Table S5). In the structurally-related
Bi<sub>2</sub>S<sub>3</sub>, the displacement of the cation from the ideal center of the octahedron has been
attributed to the effect of the 6s<sup>2</sup> lone pair.<sup>63</sup> The short Bi-S distances, which vary between
2.662(2) and 2.759(2) Å, are comparable to the sum of the covalent radii<sup>64</sup> for Bi and S (ca.
2.53 Å). By contrast, the three longer Bi-S distances, which range between 2.962(2) and
3.145(2) Å, are larger than the sum of ionic radii for Bi<sup>3+</sup> and S<sup>2-</sup>, which is ca. 2.87 Å (r(Bi<sup>3+</sup>)
= 1.03 Å and r(S<sup>2-</sup>) = 1.84 Å).<sup>65</sup> This indicates that the bonding environment around the Bi<sup>3+</sup>
cation is heterogeneous, comprising both weaker and stronger bonds (see bond valence sums,
SI Table S4). The AIMD simulations are in excellent agreement with the structural analysis:
the Bi-S radial pair distribution function (RDF) obtained from NVT AIMD shows a first peak
centered on the shorter bond lengths and a shoulder related to the longer Bi-S distances
(Figure 3(c)); the asymmetry of this peak is consistent with a stereochemically active lone
pair, which is evident in the electron localization function (ELF) (Figure 4). The magnitude
of the Bader charges (*vide infra*) also indicates substantial deviations from a purely ionic
picture of the bonding.

The Cu⁺ cations adopt a distorted tetrahedral coordination (Figure 3(d)), with corner-sharing
[CuS₄]⁷⁻ tetrahedra forming chains parallel to the *b*-axis (Figure 3(e)). The Cu-S distances
vary between 2.320(2) and 2.419(2) Å, which is consistent with the computed Cu-S RDF (Figure 3(f)). The significant deviation of the Bader charges (vide infra) for Cu⁺ and S²⁻ from the formal oxidation states suggests a large degree of covalency for the Cu-S bonds. Each Cu⁺ cation is surrounded by three Pb²⁺ cations at a distance of ca. 3.3 Å, which is smaller than the sum of their van der Waals radii⁶⁶ (ca. 3.8 Å) and suggests a possible interaction between the 6s² lone pairs of Pb²⁺ and the Cu⁺ cations. Although longer than the Pb-Cu distances, there are also Cu-Bi distances of ca. 3.5 and 3.7 Å, which are below the sum of the van der Waals radii⁶⁶ for Bi and Cu (ca. 3.9 Å). These cation-cation distances are also evident in the Cu-Bi and Cu-Pb RDFs (SI, Figure S3). The Pb²⁺ cation is coordinated to seven S²⁻ anions, forming a capped octahedron (Figure 3(g)). Each [PbS₇]¹²⁻ capped octahedron shares faces with four other capped octahedra, to form ribbons, with stoichiometry [PbS₃]⁻, oriented along the *b*-axis (Figure 3(h)). The Pb-S distances, which range from 2.885(2) to 3.274(2) Å, are comparable to the sum of the ionic radii, which is ca. 3.07 Å (r (Pb²⁺) = 1.23 Å),⁶⁵ and in agreement with the AIMD RDF (Figure 3(i)). In addition, each Pb²⁺ cation has three neighboring Cu⁺ cations at approximately 3.3 Å (Movie S1).

<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Contour of the electron localization function (ELF) on the [011], [001] (top) and [-211] (bottom) planes. Values for the ELF range between 0 and 1: ELF = 0.5 (blue contours) indicates free electron behavior and ELF = 1.0 (red contour) indicates perfect localization. Values smaller than 0.5 are less significant and usually point to small local electron density. (b) *Ab initio* molecular dynamics (AIMD) snapshot highlighting the Pb²⁺</figcaption></figure>
(black) – Cu⁺ (blue) coordination environment. The white sphere indicates the center of the maximally localized Wannier functions (MLWFC) associated with the Pb²⁺ lone pair. The dipole moment points from the Pb²⁺ position to the MLWFC. The θ, φ angles are used to indicate the dipole orientation; φ is in the plane in which the Pb²⁺ and one Cu⁺ lie, θ is in the perpendicular plane. (c) Rotational time correlation function (TCF) for the MLWFC dipole moment of the Pb²⁺ lone pair.

Analysis of the electron localization function (ELF) (Figure 4(a)), which can be used to determine if the bonding interactions involve shared electrons (e.g. covalent bonding) or unshared electrons (e.g. ionic bonding),⁶⁷ reveals the presence of directional covalent bonds between Bi³⁺ and S²⁻, as well as that of a lone pair on the Bi³⁺ cations. By contrast, the distribution around the Pb²⁺ cations in the ELF contours is reasonably isotropic and spherical. This is consistent with weak electrostatic interactions between the Pb²⁺ cations and neighboring atoms, and suggests that Pb²⁺ could be considered to be in a quasi-liquid state. The isotropic ELF can arise from dynamic behavior of the Pb²⁺ lone pair, which would give rise to rotational motion, as recently found in halide perovskites.⁶⁸

The presence of different types of bonding (bond heterogeneity) has been identified as a characteristic feature that can lead to increased phonon scattering, resulting in low thermal conductivity.⁶⁹⁻⁷¹ Moreover, in agreement with a previous single-crystal X-ray diffraction study,⁶⁰ the atomic displacement parameters (Uᵢₐ) for the Cu⁺ and Pb²⁺ cations are larger than those of Bi³⁺ and S²⁻ (Table 1), and the RDFs involving Cu⁺ and Pb²⁺ (Figure 3(f) and (i) and Figure S3), exhibit marked peak broadening with increasing temperature. The effect of temperature has been investigated by examining the displacements of the cations with respect to their equilibrium positions (SI, Figure S4). This shows a larger response of the Cu⁺ and Pb²⁺ cations to increases in temperature. A large atomic displacement parameter has been related to weak interatomic bonding and ‘rattling’-like vibrations,¹⁹,⁷² or to an underlying distortion at the local scale, arising from uncorrelated lone pair stereochemical activity.⁷³ For aikinite, analysis of the AIMD trajectories performed with Wannier functions provides clear evidence of a dynamical effect arising from rotation of the Pb²⁺ lone pair (Figure 4 and SI, Figure S5 plus Movie S2). This analysis suggests that the large atomic displacement parameters for Cu⁺ and Pb²⁺ arise from the incoherent rotation of the Pb²⁺ lone pair, which is accompanied by cooperative displacements of the Pb²⁺ cation towards one of the three neighboring Cu⁺ cations, due to attractive electrostatic interaction between the lone pair and the Cu⁺ cation. From the statistical analysis of the AIMD trajectories, we determined the maxima of the angular amplitude for the Pb²⁺ lone pair MLWFC dipole moment; the values
are $\phi_{max} = 360^{\circ}$ and $\theta_{max} = 134^{\circ}$. These angles are consistent with rotation of the lone pair between the three Cu$^{+}$ cations.

**Electronic structure and transport properties:** Aikinite is a semiconductor with a theoretical indirect band gap of 0.8 eV (Figure 5), which is in close agreement with the indirect band gap estimated from optical measurements, *ca.* 0.91 eV (SI, Figure S6). The computed partial density of states (Figure 5, right panel) shows that the major contribution at the top of the valence manifold arises, as expected, from the Cu$^{+}$ and S$^{2-}$ ions, which form corner-sharing tetrahedral chains. The top of the valence band has a multivalley character with large effective masses (Table S6) at the $\Gamma$, $X$, and $Z$ high-symmetry points of the Brillouin zone (Figure 5, left panel). The effective masses, computed over the entire band structure (Table S6), are indicative of a high degree of anisotropy, which is not immediately evident from Figure 5, which shows the band structure along the high symmetry path of the Brillouin zone of aikinite. The dispersive valence band indicates significant deviation from ionic character, as demonstrated by the Bader charges, which are +1.3, +1.0, +0.5, and -1.0 for Bi, Pb, Cu, and S respectively at T = 0 K (interestingly the values at cT = 300 K are +1.5/+2.0, +1.0, +1.0, and -1.0/-1.5). The features of the band structure suggest a large Seebeck coefficient for *p*-type transport, arising from the large anisotropic effective masses at $Z$, and high electrical resistivity.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Electronic band structure (left panel) and atom-projected electron density of states (right panel) of aikinite. Spin-orbit coupling has been included in the calculation. Hubbard U corrections are included. The top of the valence band is set to 0 eV.</figcaption></figure>
The conclusions from the band structure calculations are consistent with the experimentally determined Seebeck coefficient (S) (Figure 6(a)), which shows large, positive, and monotonically decreasing values between 423 and 573 K, indicating that aikinite is a non-degenerate p-type semiconductor. Measurements also indicate that CuPbBi<sub>5</sub>S<sub>9</sub> is highly resistive (1.27 × 10<sup>3</sup> Ω cm at 423 K, Figure 6(a)), with a very low intrinsic charge carrier concentration, ca. 10<sup>14</sup> cm<sup>-3</sup>, comparable to those determined for other undoped members of the aikinite-bismuthinite series (4.5 × 10<sup>12</sup> cm<sup>-3</sup> for CuPbBi<sub>5</sub>S<sub>9</sub>, 3.7 × 10<sup>16</sup> for Bi<sub>2</sub>S<sub>3</sub>).<sup>31,29</sup> Unlike the intermediate composition CuPbBi<sub>5</sub>S<sub>9</sub> (x = 2/3) in the bismuthinite-aikinite series Cu<sub>1-x</sub>☐<sub>x</sub>Pb<sub>1-x</sub>Bi<sub>1+x</sub>S<sub>3</sub>,<sup>31</sup> which exhibits n-type conductivity probably because of sulfur vacancies, aikinite CuPbBi<sub>5</sub>S<sub>9</sub> (x = 0), as prepared following the synthesis route described herein, retains p-type conductivity even after repeated heated-cooling cycles (SI, Figure S7). TGA measurements (SI, Figure S8) confirm the stability of this material, which is thermally stable up to 800 K under an inert atmosphere, whilst DSC data are consistent with the absence of phase transitions (SI, Figure S9). While in aikinite there is a continuous one-dimensional network of corner-sharing [CuS<sub>4</sub>]<sup>7-</sup> tetrahedra which facilitates p-type electrical conduction, other members of the bismuthinite-aikinite series, Cu<sub>1-x</sub>☐<sub>x</sub>Pb<sub>1-x</sub>Bi<sub>1+x</sub>S<sub>3</sub>, with x > 0 contain vacant sites in the one-dimensional [Cu<sub>1-x</sub>☐<sub>x</sub>S<sub>3</sub>]<sup>(5+x)-</sup> chains, instead of a continuous network. Therefore, in copper-deficient materials such as CuPbBi<sub>5</sub>S<sub>9</sub>, n-type electrical conduction involves the Bi-S network, as previously discussed by Maji *et al*.<sup>31</sup>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) Temperature dependence of the electrical resistivity (ρ) and the Seebeck coefficient (S) for aikinite, consistent with p-type semiconducting behaviour. (b) Temperature dependence of the total thermal conductivity (κ) of aikinite. The inset shows a linear fit (black line) of the lattice thermal conductivity (κ<sub>L</sub>) plotted as a function of inverse temperature.</figcaption></figure>
A study of the formation of defects in Bi₂S₃ (x = 1) has shown that for the S rich limit, Bi vacancies (□''') and S anti-site defects (S''''') are dominant, while in the Bi rich limit, the dominant defects are S vacancies (□•) and Bi anti-site defects (Bi•••••).<sup>74</sup> The situation for CuPbBiS₃ is more complex. Acceptor defects that can occur in aikinite include copper vacancies (□'Cu) usually observed in chalcogenides, lead vacancies (□''Pb), antisite defects of Cu at Pb site (Cu'Pb) or Pb on Bi site (Pb'Bi). In addition, bismuth vacancies (□'''), owing to the volatility of bismuth, could also form, although under the synthesis conditions (low temperature annealing and hot-pressing) used here this is considered unlikely. A detailed study of the energetics of defect formation is ongoing.

**Lattice thermal conductivity, elastic and vibrational properties:** The total thermal conductivity of CuPbBiS₃ (Figure 6(b)) is extremely low (~0.5 W m⁻¹ K⁻¹), and is dominated by the lattice component, κL (SI, Figure S10(a)), which constitutes approximately 99% of the total thermal conductivity. κL follows a T⁻¹ law only at temperatures above 450 K, indicating a significant contribution from Umklapp type scattering above this temperature (Figure 6(b), inset). Despite the differences in cation ordering between CuPbBi₅S₉, in which there is disorder between Cu⁺ and vacancies on the copper site and between Pb²⁺ and Bi³⁺ cations on the M1 site,³¹ and that of CuPbBiS₃, in which Cu⁺, Pb²⁺ and Bi³⁺ cations are fully ordered, the thermal conductivities of CuPbBiS₃ and CuPbBi₅S₉ are very similar. By contrast, Bi₂S₃ (x = 1) exhibits significantly larger thermal conductivities; for highly oriented ingots, the thermal conductivity is ~ 1.3 W m⁻¹ K⁻¹ along the *b*-axis and ~ 0.9 W m⁻¹ K⁻¹ along *a* and *c*,²⁹ while polycrystalline Bi₂S₃ exhibits a total thermal conductivity of ~ 0.87 W m⁻¹ K⁻¹ at room temperature.⁷⁴ Therefore, the reduction in thermal conductivity that occurs in Cu₁₋ₓBiₓPb₁₋ₓBi₁+xS₃ for x < 1 suggests that Cu⁺ and Pb²⁺ cations play a key role in the heat transport. To investigate the origin of the ultralow thermal conductivity in aikinite, sound velocity measurements were performed, from which elastic properties were derived (Table S7). The sound velocities obtained computationally are in good agreement with the experimental values (Table 2). Both the transverse (1560 m s⁻¹) and the longitudinal sound velocities (2771 m s⁻¹), which can be related to the group velocities of the heat-carrying acoustic phonons, are low and comparable to those for CuPbBi₅S₉ (x = 2/3).³¹,³² Since it has been shown⁷⁵ that, above the Debye temperature, κL is directly proportional to the cube of the average sound velocity (κL = A vₐ³/T, where A is a proportionality constant), a low sound velocity is expected to result in low thermal conductivity. For aikinite, the minimum value of κL, according to the Cahill-Watson-Pohl (CWP) model,³⁷ which describes the limit for
amorphous and strongly disordered materials, is ca. 0.41 W m⁻¹ K⁻¹. This is similar to the experimental value of 0.48 at 573 K determined here. The minimum $κ_L$ can also be calculated by considering a diffusive mechanism in which non-propagating (i.e. not phononic) atomic vibrations, known as diffusons, carry heat by diffusion.⁷⁶ This leads to an estimate of the diffuson thermal conductivity, $κ_{diff}$ ~ 0.26 W m⁻¹ K⁻¹. This represents the limit for entirely diffusive mediated transport, and is significantly lower than the values found for aikinite, providing a strong indicator that phonons contribute to heat transport in aikinite. Taking into account that $κ_L = \frac{1}{3}C_p v_a l$ (where $C_p$ is the heat capacity per unit volume and $l$ is the phonon mean free path), the estimated phonon mean free path is ca. 5 Å, which is approximately twice the interatomic spacing in aikinite, and comparable to the b lattice parameter.

The Grüneisen parameter derived from the sound velocity measurements is large, γ~1.59, indicating a high degree of anharmonicity, and comparable to values found for other thermoelectric materials with low thermal conductivities (SI, Table S8). The Grüneisen parameter obtained here is similar to that reported for other members of the Cu₁⁻ₓ□ₓPb₁⁻ₓBi₁+xS₃ series (SI, Table S9). Anharmonicity enhances phonon-phonon scattering processes, hence lowering the lattice thermal conductivity. The Young's modulus (E), extracted from the sound velocities (Table 2), is rather low (SI, Table S8), and similar to values reported for other members of the Cu₁⁻ₓ□ₓPb₁⁻ₓBi₁+xS₃ series (SI, Table S9). As the Young's modulus is related to the stiffness of the atomic bonds, this is indicative of weaker interatomic bonding. This supports the conclusions drawn from the structural analysis.

<table><caption>**Table 2.** Experimentally and computationally determined sound velocities and elastic properties, derived from the sound velocities, for aikinite.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Sound velocity (m s<sup>-1</sup>)</th><th colspan="3">Derived parameters</th></tr><tr><th>Transverse<br>ν<sub>T</sub></th><th>Longitudinal<br>ν<sub>L</sub></th><th>Average<br>sound<br>velocity<br>ν<sub>a</sub></th><th>Poisson's<br>ratio</th><th>Young's<br>modulus<br>(GPa)</th><th>Grüneisen<br>parameter</th><th>Debye<br>temperature<br>θ<sub>D</sub><br>(K)</th></tr></thead><tbody><tr><td>Exp</td><td>1560</td><td>2771</td><td>1736</td><td>0.27</td><td>42.7</td><td>1.59</td><td>183</td></tr><tr><td>Th (phonon)<sup>a</sup></td><td>1599<br>1865</td><td>3123</td><td>1929</td><td>0.28</td><td></td><td>1.64</td><td></td></tr><tr><td>Th (elastic)<sup>b</sup></td><td>1733</td><td>3199</td><td>1933</td><td>0.29</td><td></td><td>1.73</td><td></td></tr></tbody></table>

<footnote>a Determined using the phonon dispersion; b determined using the elastic constants.</footnote>
**Figure 7** presents the atom-resolved phonon density of states (vDOS) computed at 100, 200 and 300 K by AIMD, compared with the vDOS at 0 K³¹ calculated using Quantum ESPRESSO. These results should also be compared with the dispersion curves presented in Ref. 31 which show optical modes with very low frequencies (~20-50 cm⁻¹). It has been shown that a low cutoff frequency of acoustic phonons, which can be ascribed to weak bonding and a correspondingly low sound velocity, is a good indicator for low thermal conductivity.⁷⁷ Moreover, the presence of low frequency optical modes close to the acoustic mode frequencies, as is the case here, can lead to phonon scattering and affect thermal transport.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Total (red line) and atom-resolved vDOS of aikinite at different temperatures. Blue, black, green and purple lines correspond to Cu⁺, Pb²⁺, S²⁻ and Bi³⁺ ions respectively. The vDOS at 0 K was computed by DFT, and those between 100 and 300 K by AIMD. (b) Neutron-weighted total calculated vDOS convolved with the LET instrumental resolution function, as a function of temperature. Arrows indicate shifts in peak positions.</figcaption></figure>
Low frequency optical modes in the range 20 - 50 cm⁻¹ and Einstein-like modes centered around 100 cm⁻¹ are evident in the vDOS of aikinite (Figure 7). The atom-resolved vDOS shows that the main contribution at low frequencies is from Pb²⁺, with a smaller contribution from Bi³⁺. It is interesting to note that the difference between the Bi³⁺ and Pb²⁺ modes at low frequencies cannot be accounted for in terms of the variation in atomic mass (208.9 and 207.2 amu respectively), but is consistent with the weak bonding of Pb²⁺ highlighted by the structural analysis. In addition, analysis of the vDOS computed by AIMD as a function of temperature (Figure 7 (a)) indicates that, with increasing temperature, the Cu⁺ vibrational modes spread and shift to lower frequencies, overlapping more with phonons with a large Pb²⁺ component, while the Bi³⁺ contribution remains largely unchanged (Figure 7 (a)). The shift, with increasing temperature, of the Cu⁺ vibrational modes to lower frequencies, and of the peak at ca. 31 cm⁻¹ towards higher frequencies is also evident in the neutron-weighted total calculated vDOS convolved with the LET instrumental resolution function (Figure 7(b)).

In order to confirm, experimentally, the presence of the low energy phonon modes of Pb²⁺ and the Einstein-like modes arising from the Cu⁺ vibrations described above, temperature-dependent INS data were collected. The calculated neutron-weighted vDOS (convolved with a Gaussian function approximating the instrumental resolution) and the experimental vDOS are in very good agreement (Figure 8(a) and SI, Figure S11). The slight shift in frequencies between the experimental and calculated vDOS is due to the PBE functional used in the calculations, which is known to underestimate bond strengths, and the slight differences in peak intensities are related to approximations made when convolving the INS instrumental resolution with the calculated vDOS. Comparison with the partial calculated vDOS (Figure 7(a)), allows us to assign the first peak in the experimental vDOS, which at 10 K is centered at 31.7 cm⁻¹, to primarily Pb²⁺-based vibrations. Of the other five peaks present in the INS data, most have contributions from multiple atoms, and only that at 110.7 cm⁻¹ can be assigned to the Einstein-like modes of Cu⁺. Examination of the experimental vDOS collected as a function of temperature (Figure 8(b)), shows that peaks shift significantly with temperature (Table S10). In particular, with increasing temperature (Figure 8(c)), the Cu⁺ mode softens, in agreement with the findings of the computed vDOS. By contrast, the phonon peak corresponding to the Pb²⁺-based mode shifts to higher frequency (Figure 8(c)). Such ‘hardening’ of the Pb²⁺ mode with increasing temperature is indicative of anharmonic behavior and is therefore likely be a key contributor to the ultralow thermal conductivity of
aikinite. To evaluate the role of Pb²⁺ in lowering the lattice thermal conductivity, we used the quasi-harmonic approximation to estimate the degree of anharmonicity through the mode-resolved Grüneisen parameter, γᵢₐ, and its atom projection as a function of the frequency (Figure 9). This analysis provides clear evidence for significant anharmonicity. While in aikinite anharmonicity could be associated with the presence of the 6s² lone pair of electrons found in Pb²⁺ and Bi³⁺ cations, examination of Figure 9 reveals that at low frequencies, over the range 20 - 50 cm⁻¹, the magnitude of mode-resolved Grüneisen parameter is significantly higher for Pb²⁺ than for the other atoms. Finally, the rotational TCF (Figure 4(c)) associated with the Pb²⁺ lone pair shows a rotational time scale of ~ 80 fs; this corresponds to a frequency of ~ 66 cm⁻¹, comparable with the low-frequency mode of Cu⁺. This observation is consistent with the cooperative interaction between the rotation of the Pb²⁺ lone pair and the Cu⁺ cations discussed earlier.

Other members of the bismuthinite-aikinite series, Cu₁₋ₓ☐ₓPb₁₋ₓBi₁₊ₓS₃, are isostructural with aikinite (x = 0), and when x < 1, also contain Cu⁺ and Pb²⁺ cations at relatively short distances, which could facilitate coupling of rotating lone pairs and the vibrational motion of the Cu⁺ cations. Moreover, based on the limited data available in the literature, members of this series with x < 1 exhibit similarly low values of thermal conductivity to that of aikinite. This ultralow thermal conductivity has previously been ascribed to the complex crystal structure with a large unit cell and heavy atoms,³² or to disorder.⁷⁸ However, given that aikinite has a ultralow thermal conductivity even with complete ordering of the Cu⁺, Pb²⁺ and Bi³⁺ cations, we can effectively discount disorder as the origin of the ultralow thermal conductivity in other members of the Cu₁₋ₓ☐ₓPb₁₋ₓBi₁₊ₓS₃ series with x < 1. Furthermore, the larger thermal conductivity of Bi₂S₃ (x = 1) when compared to other members of the series indicates that the Cu⁺ and Pb²⁺ cations play a key role in the heat transport. Together these observations suggest that the mechanism of lone pair rotation we have unveiled in aikinite, which we have demonstrated is a key contributor to the scattering effects that lower the thermal conductivity, is at play in other members of the bismuthinite-aikinite series, Cu₁₋ₓ☐ₓPb₁₋ₓBi₁₊ₓS₃.
<figure><img src="image_10.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c)</figcaption></figure>

<caption>(d)</caption>
<figure><img src="image_12.png" /><figcaption>**Figure 8:** (a) The calculated neutron-weighted total vDOS (gⁿ(E)) (convolved with the instrumental resolution) at 0 K (blue line) and experimental neutron-weighted vDOS (black) at 10 K. Data have been normalized. (b) Neutron scattering intensities as a function of temperature. Gaussian fits for the six peaks are shown. The Gaussian corresponding to the Pb²⁺-based vibration is shown as a black line and that for the Cu⁺ rattling vibration is shown as a blue line. Experimental data have been offset along the y axis with respect to the Gaussian fits for clarity. Red dotted line shows the fitted background. Arrows highlight the change in peak position for the Pb²⁺-based and Cu⁺ rattling vibrations with temperature. (c) Percentage change in the energy of the phonon mode between 10 K and a temperature T. The grey line corresponds to the Pb²⁺ mode (ω = 32 cm⁻¹) and the blue line corresponds to the Cu⁺ mode (ω = 110 cm⁻¹). (d) Temperature dependence of the phonon lifetime of the Pb²⁺ mode.</figcaption></figure>

It is extremely challenging to extract lifetimes from INS data collected on a powder, owing to the averaging over Brillouin zones that occurs in such an experiment.⁶⁹ This is especially true when, with increasing temperature, modes are shifting very significantly in frequencies and in different directions, with corresponding changes to the underlying dispersion. As a result, peaks might sharpen or broaden independently of broadening from phonon-phonon scattering. In the INS data for aikinite, the full width at half maxima (FWHM) of the majority of the peaks decreases with increasing temperature (SI Table S11), indicating that the modes sharpen at higher temperatures. Such a sharpening could be related to a flattening of the dispersion and thereby a reduction in the optical mode group velocities. However, for Pb²⁺, the low-energy mode broadens with increasing temperature (Table 3 and SI Table S11). While it is not possible to disentangle the likely changes to the dispersion from phonon scattering, if we assume that the broadening is all from phonon-phonon scattering, we can extract a rather low lifetime of 0.41(3) ps at 300 K (Table 3, Figure 8(d)). This compares favourably with the values of ~3 ps, ~0.66 ps and ~2ps for PbTe,⁷⁹ TlInTe₂⁶⁹ or Na₀.₈CoO₂,⁸⁰ respectively. While the lifetime determined here should be considered a lower bound, since it is likely that at least some of the broadening is due to changes in the dispersion, this lifetime
would correspond to a phonon mean-free path ($l = v_a \tau$) of ca. 7 Å, which is in reasonable agreement with the estimated mean-free path based on the lattice thermal conductivity.

<figure><img src="image_13.png" /></figure>

**Figure 9.** Contribution of each atomic species to the total Grüneisen parameter as a function of the mode frequency. Red, green, blue and yellow points correspond to Bi$^{3+}$, Cu$^{+}$, Pb$^{2+}$ and S$^{2-}$ contributions, respectively.

<table><caption>**Table 3.** Phonon lifetime ($\tau$) for the Pb$^{2+}$ mode, calculated using the relation $\tau = \frac{1}{\pi(\Delta_{FWHM})}$ where $\Delta_{FWHM} = \sqrt{(FWHM_T)^2 - (FWHM_{10K})^2}$ $^{69,81}$</caption><thead><tr><th>T (K)</th><th>Peak centre (cm<sup>-1</sup>)</th><th>FWHM (cm<sup>-1</sup>)</th><th>Δ<sub>FWHM</sub> (cm<sup>-1</sup>)</th><th>τ (ps)</th></tr></thead><tbody><tr><td>10</td><td>31.7(4)</td><td>11.6(5)</td><td></td><td></td></tr><tr><td>100</td><td>36.5(4)</td><td>24.7(8)</td><td>21.7(4)</td><td>0.48(7)</td></tr><tr><td>200</td><td>36.9(3)</td><td>26.4(5)</td><td>23.7(3)</td><td>0.44(5)</td></tr><tr><td>300</td><td>37.3(2)</td><td>28.3(4)</td><td>25.8(3)</td><td>0.41(3)</td></tr></tbody></table>

## CONCLUSIONS

In aikinite, the combination of bond heterogeneity and the presence of heavy metal cations with stereochemically-active lone pairs, leads to an exceptionally low thermal conductivity.
Moreover, the Cu⁺ and Pb²⁺ cations have large atomic displacement parameters and contribute to the same low-frequency vibrational manifold. With increasing temperature, the Cu⁺ contributions in the low-frequency region increase; this is consistent with a thermally activated interaction between Cu⁺ and Pb²⁺ cations. Using *ab initio* molecular dynamics and Wannier functions analysis, we have characterized the weak and isotropic bonding of the Pb²⁺ cations that facilitates the rotation of the 6s² lone pair. In turn, such rotations influence the Cu⁺ dynamics and reduce the thermal conductivity. The interaction mechanism can be rationalized by considering the electrostatic attraction between the thermally activated lone pairs on the Pb²⁺ cations, which are rotating, and the Cu⁺ cations. Given that synthetic samples in the bismuthinite-aikinite series, Cu₁₋ₓ☐ₓPb₁₋ₓBi₁₊ₓS₃, are isostructural with aikinite (x = 0), it is highly likely that the ultralow thermal conductivity observed for this family of materials is also a consequence of the cooperative interaction between the rotating lone pair on the Pb²⁺ cations and the Cu⁺ cations.

Lone pair rotation in solids is an emergent phenomenon, recently reported in halide perovskites, although it has been suggested that such rotational motion may also occur in other solids containing lone-pair electrons.⁵⁷ The work presented here provides clear evidence for the dynamical behavior of the Pb²⁺ lone pair and its contribution to lowering thermal conductivity, in a family of thermoelectric sulfides, Cu₁₋ₓ☐ₓPb₁₋ₓBi₁₊ₓS₃. We have demonstrated, for the first time, that the coupling of rotating lone pairs and vibrational motion is an effective mechanism to achieve ultralow thermal conductivity in crystalline materials. We suggest that this mechanism may also occur in other families of materials, including halide perovskites such as CsSnBr₃₋ₓIₓ, which also exhibit ultralow thermal conductivities.⁸² Therefore, the results presented here offer new insights for the search of materials with ultralow thermal conductivity.

## ACKNOWLEDGMENTS

The authors acknowledge The Leverhulme Trust for Research Project Grant RPG-2019-288, which funded this work. We would also like to acknowledge the STFC ISIS Neutron and Muon Source for provision of beam time on the LET beamline under proposal RB2210076. This research used resources at the Spallation Neutron Source (experiment IPTS-29028.1), a DOE Office of Science User Facility operated by the Oak Ridge National Laboratory. We also thank the Chemical Analysis Facility (CAF) at the University of Reading, Reading, UK, for access to powder X-ray diffraction and thermal analysis instrumentation. Calculations
were performed at the High-Performance Computing Center at Michigan State University.
We thank Rajendra Paudel for useful discussions.

## SUPPORTING INFORMATION

The Supporting Information is available free of charge at http://******

Details of the methodology, Helmholtz free energy, internal energy, entropy and specific heat, Rietveld refinements using different models, bond distances, bond valence sums and distortion parameters, RDF plots, ion displacement plots, diffuse reflectance, TGA, DSC, heat capacity, electrical and thermal transport data, calculated effective masses, computed, vibrational DOS determined using inelastic neutron scattering data and derived parameters, computed elastic constants, AIMD movie showing Cu-Pb distances, AIMD movie showing Pb²⁺ lone pair rotation, Pb²⁺ lone pair Wannier function isosurface.

## REFERENCES

(1) Padture, N. P.; Gell, M.; Jordan, E. H. Thermal Barrier Coatings for Gas-Turbine Engine Applications. Science 2002, 296 (5566), 280-284. https://doi.org/10.1126/science.1068609.

(2) Wang, Q.; Jiang, B.; Li, B.; Yan, Y. A Critical Review of Thermal Management Models and Solutions of Lithium-Ion Batteries for the Development of Pure Electric Vehicles. Renew. Sustain. Energy Rev. 2016, 64, 106-128. https://doi.org/10.1016/j.rser.2016.05.033.

(3) Ghosh, T.; Dutta, M.; Sarkar, D.; Biswas, K. Insights into Low Thermal Conductivity in Inorganic Materials for Thermoelectrics. J. Am. Chem. Soc. 2022, 144 (23), 10099-10118. https://doi.org/10.1021/jacs.2c02017.

(4) Biswas, K.; He, J.; Blum, I. D.; Wu, C.-I.; Hogan, T. P.; Seidman, D. N.; Dravid, V. P.; Kanatzidis, M. G. High-Performance Bulk Thermoelectrics with All-Scale Hierarchical Architectures. Nature 2012, 489 (7416), 414-418. https://doi.org/10.1038/nature11439.

(5) Kim, S. I.; Lee, K. H.; Mun, H. A.; Kim, H. S.; Hwang, S. W.; Roh, J. W.; Yang, D. J.; Shin, W. H.; Li, X. S.; Lee, Y. H.; Snyder, G. J.; Kim, S. W. Dense Dislocation Arrays Embedded in Grain Boundaries for High-Performance Bulk Thermoelectrics. Science 2015, 348 (6230), 109-114. https://doi.org/10.1126/science.aaa4166.

(6) Hsu, K. F.; Loo, S.; Guo, F.; Chen, W.; Dyck, J. S.; Uher, C.; Hogan, T.; Polychroniadis, E. K.; Kanatzidis, M. G. Cubic AgPb<sub>m</sub>SbTe<sub>2+m</sub>: Bulk Thermoelectric Materials with High Figure of Merit. Science 2004, 303 (5659), 818-821. https://doi.org/10.1126/science.1092963.

(7) Suekuni, K.; Kim, F. S.; Nishiate, H.; Ohta, M.; Tanaka, H. I.; Takabatake, T. High-Performance Thermoelectric Minerals: Colusites Cu<sub>26</sub>V<sub>2</sub>M<sub>6</sub>S<sub>32</sub> (M = Ge, Sn). Appl. Phys. Lett. 2014, 105 (13), 132107. https://doi.org/10.1063/1.4896998.
(8) Kumar, V. P.; Lemoine, P.; Carnevali, V.; Guélou, G.; Lebedev, O. I.; Raveau, B.; Al Rahal Al Orabi, R.; Fornari, M.; Candolfi, C.; Prestipino, C.; Menut, D.; Malaman, B.; Juraszek, J.; Suekuni, K.; Guilmeau, E. Local-Disorder-Induced Low Thermal Conductivity in Degenerate Semiconductor Cu<sub>22</sub>Sn<sub>10</sub>S<sub>32</sub>. Inorg. Chem. 2021, 60 (21), 16273–16285.
https://doi.org/10.1021/acs.inorgchem.1c02105.

(9) Guélou, G.; Pavan Kumar, V.; Carnevali, V.; Lebedev, O. I.; Raveau, B.; Couder, C.; Prestipino, C.; Lemoine, P.; Malaman, B.; Juraszek, J.; Candolfi, C.; Lenoir, B.; Al Rahal Al Orabi, R.; Fornari, M.; Guilmeau, E. Long-Range Cationic Order Collapse Triggered by S/Cl Mixed-Anion Occupancy Yields Enhanced Thermoelectric Properties in Cu<sub>5</sub>Sn<sub>2</sub>S<sub>7</sub>. Chem. Mater. 2021, 33 (23), 9425–9438. https://doi.org/10.1021/acs.chemmater.1c03434.

(10) Liu, H.; Shi, X.; Xu, F.; Zhang, L.; Zhang, W.; Chen, L.; Li, Q.; Uher, C.; Day, T.; Snyder, G. J. Copper Ion Liquid-like Thermoelectrics. Nat. Mater. 2012, 11 (5), 422–425.
https://doi.org/10.1038/nmat3273.

(11) Lai, W.; Wang, Y.; Morelli, D. T.; Lu, X. From Bonding Asymmetry to Anharmonic Rattling in Cu<sub>12</sub>Sb<sub>4</sub>S<sub>13</sub> Tetrahedrites: When Lone-Pair Electrons Are Not So Lonely. Adv. Funct. Mater. 2015, 25 (24), 3648–3657. https://doi.org/10.1002/adfm.201500766.

(12) Lee, S.; Esfarjani, K.; Luo, T.; Zhou, J.; Tian, Z.; Chen, G. Resonant Bonding Leads to Low Lattice Thermal Conductivity. Nat. Commun. 2014, 5 (1), 3525.
https://doi.org/10.1038/ncomms4525.

(13) Nielsen, M. D.; Ozolins, V.; Heremans, J. P. Lone Pair Electrons Minimize Lattice Thermal Conductivity. Energy Environ. Sci. 2013, 6 (2), 570–578. https://doi.org/10.1039/C2EE23391F.

(14) Sato, N.; Kuroda, N.; Nakamura, S.; Katsura, Y.; Kanazawa, I.; Kimura, K.; Mori, T. Bonding Heterogeneity in Mixed-Anion Compounds Realizes Ultralow Lattice Thermal Conductivity. J. Mater. Chem. A 2021, 9 (39), 22660–22669. https://doi.org/10.1039/D1TA04958E.

(15) Dutta, M.; Pal, K.; V. Waghmare, U.; Biswas, K. Bonding Heterogeneity and Lone Pair Induced Anharmonicity Resulted in Ultralow Thermal Conductivity and Promising Thermoelectric Properties in N-Type AgPbBiSe<sub>3</sub>. Chem. Sci. 2019, 10 (18), 4905–4913.
https://doi.org/10.1039/C9SC00485H.

(16) Vaqueiro, P. Chapter 1:Synthesis and Property Measurements of Thermoelectric Materials. In Inorganic Thermoelectric Materials; 2021; pp 1–52. https://doi.org/10.1039/9781788019590-00001.

(17) Powell, A. V. Recent Developments in Earth-Abundant Copper-Sulfide Thermoelectric Materials. J. Appl. Phys. 2019, 126 (10), 100901. https://doi.org/10.1063/1.5119345.

(18) Vaqueiro, P.; Guélou, G.; Kaltzoglou, A.; Smith, R. I.; Barbier, T.; Guilmeau, E.; Powell, A. V. The Influence of Mobile Copper Ions on the Glass-Like Thermal Conductivity of Copper-Rich Tetrahedrites. Chem. Mater. 2017, 29 (9), 4080–4090.
https://doi.org/10.1021/acs.chemmater.7b00891.

(19) Long, S. O.; Powell, A. V.; Hull, S.; Orlandi, F.; Tang, C. C.; Supka, A. R.; Fornari, M.; Vaqueiro, P. Jahn–Teller Driven Electronic Instability in Thermoelectric Tetrahedrite. Adv. Funct. Mater. 2020, 30 (12), 1909409. https://doi.org/10.1002/adfm.201909409.
(20) Suekuni, K.; Tsuruta, K.; Kunii, M.; Nishiate, H.; Nishibori, E.; Maki, S.; Ohta, M.; Yamamoto, A.; Koyano, M. High-Performance Thermoelectric Mineral Cu<sub>12-x</sub>Ni<sub>x</sub>Sb<sub>4</sub>S<sub>13</sub> Tetrahedrite. J. Appl. Phys. 2013, 113 (4), 043712. https://doi.org/10.1063/1.4789389.

(21) Lu, X.; Morelli, D. T.; Xia, Y.; Zhou, F.; Ozolins, V.; Chi, H.; Zhou, X.; Uher, C. High Performance Thermoelectricity in Earth-Abundant Compounds Based on Natural Mineral Tetrahedrites. Adv. Energy Mater. 2013, 3 (3), 342–348.
https://doi.org/10.1002/aenm.201200650.

(22) Barbier, T.; Lemoine, P.; Gascoin, S.; Lebedev, O. I.; Kaltzoglou, A.; Vaqueiro, P.; Powell, A. V.; Smith, R. I.; Guilmeau, E. Structural Stability of the Synthetic Thermoelectric Ternary and Nickel-Substituted Tetrahedrite Phases. J. Alloys Compd. 2015, 634, 253–262.
https://doi.org/10.1016/j.jallcom.2015.02.045.

(23) Bourgès, C.; Bouyrie, Y.; Supka, A. R.; Al Rahal Al Orabi, R.; Lemoine, P.; Lebedev, O. I.; Ohta, M.; Suekuni, K.; Nassif, V.; Hardy, V.; Daou, R.; Miyazaki, Y.; Fornari, M.; Guilmeau, E. High-Performance Thermoelectric Bulk Colusite by Process Controlled Structural Disordering. J. Am. Chem. Soc. 2018, 140 (6), 2186–2195.
https://doi.org/10.1021/jacs.7b11224.

(24) Pavan Kumar, V.; Guélou, G.; Lemoine, P.; Raveau, B.; Supka, A. R.; Al Rahal Al Orabi, R.; Fornari, M.; Suekuni, K.; Guilmeau, E. Copper- Rich Thermoelectric Sulfides: Size- Mismatch Effect and Chemical Disorder in the [TS<sub>4</sub>]Cu<sub>6</sub> Complexes of Cu<sub>26</sub>T<sub>2</sub>Ge<sub>6</sub>S<sub>32</sub> (T = Cr, Mo, W) Colusites. Angew. Chem. Int. Ed. 2019, 58 (43), 15455–15463.
https://doi.org/10.1002/anie.201908579.

(25) Pavan Kumar, V.; Supka, A. R.; Lemoine, P.; Lebedev, O. I.; Raveau, B.; Suekuni, K.; Nassif, V.; Al Rahal Al Orabi, R.; Fornari, M.; Guilmeau, E. High Power Factors of Thermoelectric Colusites Cu<sub>26</sub>T<sub>2</sub>Ge<sub>6</sub>S<sub>32</sub> (T = Cr, Mo, W): Toward Functionalization of the Conductive “Cu–S” Network. Adv. Energy Mater. 2019, 9 (6), 1803249. https://doi.org/10.1002/aenm.201803249.

(26) Guélou, G.; Lemoine, P.; Raveau, B.; Guilmeau, E. Recent Developments in High-Performance Thermoelectric Sulphides: An Overview of the Promising Synthetic Colusites. J. Mater. Chem. C 2021, 9 (3), 773–795. https://doi.org/10.1039/D0TC05086E.

(27) Long, S. O. J.; Powell, A. V.; Vaqueiro, P.; Hull, S. High Thermoelectric Performance of Bornite through Control of the Cu(II) Content and Vacancy Concentration. Chem. Mater. 2018, 30 (2), 456–464. https://doi.org/10.1021/acs.chemmater.7b04436.

(28) Freer, R.; Ekren, D.; Ghosh, T.; Biswas, K.; Qiu, P.; Wan, S.; Chen, L.; Han, S.; Fu, C.; Zhu, T.; Shawon, A. K. M. A.; Zevalkink, A.; Imasato, K.; Snyder, G. J.; Ozen, M.; Saglik, K.; Aydemir, U.; Cardoso-Gil, R.; Svanidze, E.; Funahashi, R.; Powell, A. V.; Mukherjee, S.; Tippireddy, S.; Vaqueiro, P.; Gascoin, F.; Kyratsi, T.; Sauerschnig, P.; Mori, T. Key Properties of Inorganic Thermoelectric Materials—Tables (Version 1). J. Phys. Energy 2022, 4 (2), 022002.
https://doi.org/10.1088/2515-7655/ac49dc.

(29) Biswas, K.; Zhao, L.-D.; Kanatzidis, M. G. Tellurium-Free Thermoelectric: The Anisotropic n-Type Semiconductor Bi<sub>2</sub>S<sub>3</sub>. Adv. Energy Mater. 2012, 2 (6), 634–638.
https://doi.org/10.1002/aenm.201100775.
(30) Zhao, J.; Hao, S.; Islam, S. M.; Chen, H.; Tan, G.; Ma, S.; Wolverton, C.; Kanatzidis, M. G. Six Quaternary Chalcogenides of the Pavonite Homologous Series with Ultralow Lattice Thermal Conductivity. Chem. Mater. 2019, 31 (9), 3430–3439.
https://doi.org/10.1021/acs.chemmater.9b00585.

(31) Maji, K.; Lemoine, P.; Renaud, A.; Zhang, B.; Zhou, X.; Carnevali, V.; Candolfi, C.; Raveau, B.; Al Rahal Al Orabi, R.; Fornari, M.; Vaqueiro, P.; Pasturel, M.; Prestipino, C.; Guilmeau, E. A Tunable Structural Family with Ultralow Thermal Conductivity: Copper-Deficient Cu₁⁻ₓ□ₓPb₁⁻ₓBi₁⁺ₓS₃. J. Am. Chem. Soc. 2022, 144 (4), 1846–1860.
https://doi.org/10.1021/jacs.1c11998.

(32) Liang, H.; Guo, J.; Zhou, Y.-X.; Wang, Z.-Y.; Feng, J.; Ge, Z.-H. CuPbBi₅S₉ Thermoelectric Material with an Intrinsic Low Thermal Conductivity: Synthesis and Properties. J. Materiomics 2022, 8 (1), 174–183. https://doi.org/10.1016/j.jmat.2021.03.016.

(33) Mumme, W. G.; Welin, E.; Wuensch, B. J. Crystal Chemistry and Proposed Nomenclature for Sulfosalts Intermediate in the System Bismuthinite-Aikinite (Bi₂S₃-CuPbBiS₃). Am. Mineral. 1976, 61 (1–2), 15–20.

(34) Harris, D. C.; Chen, T. T. Crystal Chemistry and Re-Examination of Nomenclature of Sulfosalts in the Aikinite-Bismuthinite Series. Can. Mineral. 1976, 14 (2), 194–205.

(35) Chen, T. T.; Kirchner, E.; Paar, W. Friedrichite, Cu₅Pb₅Bi₇S₁₈, a New Member of the Aikinite-Bismuthinite Series. Can. Mineral. 1978, 16 (2), 127–130.

(36) Pring, A. Annealing of Synthetic Hammarite, Cu₂Pb₂Bi₄S₉, and the Nature of Cation-Ordering Processes in the Bismuthinite-Aikinite Series. Am. Mineral. 1995, 80 (11–12), 1166–1173.
https://doi.org/10.2138/am-1995-11-1207.

(37) Cahill, D. G.; Watson, S. K.; Pohl, R. O. Lower Limit to the Thermal Conductivity of Disordered Crystals. Phys. Rev. B 1992, 46 (10), 6131–6140.
https://doi.org/10.1103/PhysRevB.46.6131.

(38) Toby, B. H. EXPGUI, a Graphical User Interface for GSAS. J. Appl. Crystallogr. 2001, 34 (2), 210–213. https://doi.org/10.1107/S0021889801002242.

(39) Huq, A.; Kirkham, M.; Peterson, P. F.; Hodges, J. P.; Whitfield, P. S.; Page, K.; Húgle, T.; Iverson, E. B.; Parizzi, A.; Rennich, G. POWGEN: Rebuild of a Third-Generation Powder Diffractometer at the Spallation Neutron Source. J. Appl. Crystallogr. 2019, 52 (5), 1189–1201.
https://doi.org/10.1107/S160057671901121X.

(40) Toby, B. H.; Von Dreele, R. B. GSAS-II: The Genesis of a Modern Open-Source All Purpose Crystallography Software Package. J. Appl. Crystallogr. 2013, 46 (2), 544–549.
https://doi.org/10.1107/S0021889813003531.

(41) Tauc, J.; Grigorovici, R.; Vancu, A. Optical Properties and Electronic Structure of Amorphous Germanium. Phys. Status Solidi B 1966, 15 (2), 627–637.
https://doi.org/10.1002/pssb.19660150224.
(42) Kim, H. S.; Gibbs, Z. M.; Tang, Y.; Wang, H.; Snyder, G. J. Characterization of Lorenz Number with Seebeck Coefficient Measurement. APL Mater. 2015, 3 (4), 1–6.
https://doi.org/10.1063/1.4908244.

(43) Vaqueiro, P.; Mukherjee, S.; Powell, A. V.; Voneshen, D. J. Ultralow thermal conductivity in aikinite, CuPbBiS₃. https://doi.org/10.5286/ISIS.E.RB2210076 (accessed 2022-08-22).

(44) Arnold, O.; Bilheux, J. C.; Borreguero, J. M.; Buts, A.; Campbell, S. I.; Chapon, L.; Doucet, M.; Draper, N.; Ferraz Leal, R.; Gigg, M. A.; Lynch, V. E.; Markvardsen, A.; Mikkelson, D. J.; Mikkelson, R. L.; Miller, R.; Palmen, K.; Parker, P.; Passos, G.; Perring, T. G.; Peterson, P. F.; Ren, S.; Reuter, M. A.; Savici, A. T.; Taylor, J. W.; Taylor, R. J.; Tolchenov, R.; Zhou, W.; Zikovsky, J. Mantid—Data Analysis and Visualization Package for Neutron Scattering and µSR Experiments. Nucl. Instrum. Methods Phys. Res. Sect. Accel. Spectrometers Detect. Assoc. Equip. 2014, 764, 156–166. https://doi.org/10.1016/j.nima.2014.07.029.

(45) Akeroyd, F.; Ansell, S.; Antony, S.; Arnold, O.; Bekasovs, A.; Bilheux, J.; Borreguero, J.; Brown, K.; Buts, A.; Campbell, S.; Champion, D.; Chapon, L.; Clarke, M.; Cottrell, S.; Dalgliesh, R.; Dillow, D.; Doucet, M.; Draper, N.; Fowler, R.; Gigg, M. A.; Granroth, G.; Hagen, M.; Heller, W.; Hillier, A.; Howells, S.; Jackson, S.; Kachere, D.; Koennecke, M.; Le Bourlot, C.; Leal, R.; Lynch, V.; Manuel, P.; Markvardsen, A.; McGreevy, R.; Mikkelson, D.; Mikkelson, R.; Miller, R.; Nagella, S.; Nielsen, T.; Palmen, K.; Parker, P. G.; Pascal, M.; Passos, G.; Perring, T.; Peterson, P. F.; Pratt, F.; Proffen, T.; Radaelli, P.; Rainey, J.; Ren, S.; Reuter, M.; Sastry, L.; Savici, A.; Taylor, J.; Taylor, R. J.; Thomas, M.; Tolchenov, R.; Whitley, R.; Whitty, M.; Williams, S.; Zhou, W.; Zikovsky, J. Mantid: Manipulation and Analysis Toolkit for Instrument Data., 2013. https://doi.org/10.5286/SOFTWARE/MANTID.

(46) Azuah, R. T.; Kneller, L. R.; Qiu, Y.; Tregenna-Piggott, P. L. W.; Brown, C. M.; Copley, J. R. D.; Vimeo, R. M. DAVE: A Comprehensive Software Suite for the Reduction, Visualization, and Analysis of Low Energy Neutron Spectroscopic Data. J.Res.Natl.Inst.Stan.Technol 2009, 114 (341).

(47) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M. QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials. J. Phys. Condens. Matter 2009, 21 (39), 395502.
https://doi.org/10.1088/0953-8984/21/39/395502.

(48) Supka, A. R.; Lyons, T. E.; Liyanage, L.; D'Amico, P.; Al Rahal Al Orabi, R.; Mahatara, S.; Gopal, P.; Toher, C.; Ceresoli, D.; Calzolari, A.; Curtarolo, S.; Nardelli, M. B.; Fornari, M. AFLOWπ: A Minimalist Approach to High-Throughput Ab Initio Calculations Including the Generation of Tight-Binding Hamiltonians. Comput. Mater. Sci. 2017, 136, 76–84.
https://doi.org/10.1016/j.commatsci.2017.03.055.

(49) Agapito, L. A.; Curtarolo, S.; Buongiorno Nardelli, M. Reformulation of DFT + U as a Pseudohybrid Hubbard Density Functional for Accelerated Materials Discovery. Phys. Rev. X 2015, 5 (1), 011006. https://doi.org/10.1103/PhysRevX.5.011006.
(50) Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13 (12), 5188-5192. https://doi.org/10.1103/PhysRevB.13.5188.

(51) Supka, A.; Mecholsky, N. A.; Buongiorno Nardelli, M.; Curtarolo, S.; Fornari, M. Two-Layer High-Throughput: Effective Mass Calculations Including Warping. Engineering 2022, 10, 74–80. https://doi.org/10.1016/j.eng.2021.03.031.

(52) Kühne, T. D.; Iannuzzi, M.; Del Ben, M.; Rybkin, V. V.; Seewald, P.; Stein, F.; Laino, T.; Khaliullin, R. Z.; Schütt, O.; Schiffmann, F.; Golze, D.; Wilhelm, J.; Chulkov, S.; Bani-Hashemian, M. H.; Weber, V.; Borštnik, U.; Taillefumier, M.; Jakobovits, A. S.; Lazzaro, A.; Pabst, H.; Müller, T.; Schade, R.; Guidon, M.; Andermatt, S.; Holmberg, N.; Schenter, G. K.; Hehn, A.; Bussy, A.; Belleflamme, F.; Tabacchi, G.; Glöß, A.; Lass, M.; Bethune, I.; Mundy, C. J.; Plessl, C.; Watkins, M.; VandeVondele, J.; Krack, M.; Hutter, J. CP2K: An Electronic Structure and Molecular Dynamics Software Package - Quickstep: Efficient and Accurate Electronic Structure Calculations. J. Chem. Phys. 2020, 152 (19), 194103. https://doi.org/10.1063/5.0007045.

(53) VandeVondele, J.; Hutter, J. Gaussian Basis Sets for Accurate Calculations on Molecular Systems in Gas and Condensed Phases. J. Chem. Phys. 2007, 127 (11), 114105. https://doi.org/10.1063/1.2770708.

(54) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126 (1), 014101. https://doi.org/10.1063/1.2408420.

(55) Marzari, N.; Mostofi, A. A.; Yates, J. R.; Souza, I.; Vanderbilt, D. Maximally Localized Wannier Functions: Theory and Applications. Rev. Mod. Phys. 2012, 84 (4), 1419-1475. https://doi.org/10.1103/RevModPhys.84.1419.

(56) Berghold, G.; Mundy, C. J.; Romero, A. H.; Hutter, J.; Parrinello, M. General and Efficient Algorithms for Obtaining Maximally Localized Wannier Functions. Phys. Rev. B 2000, 61 (15), 10040-10048. https://doi.org/10.1103/PhysRevB.61.10040.

(57) Remsing, R. C.; Klein, M. L. Lone Pair Rotational Dynamics in Solids. Phys. Rev. Lett. 2020, 124 (6), 066001. https://doi.org/10.1103/PhysRevLett.124.066001.

(58) Zhang, Y. First-Principles Debye-Callaway Approach to Lattice Thermal Conductivity. J. Materiomics 2016, 2 (3), 237-247. https://doi.org/10.1016/j.jmat.2016.06.004.

(59) Zhang, Y.; Skoug, E.; Cain, J.; Ozoliņš, V.; Morelli, D.; Wolverton, C. First-Principles Description of Anomalously Low Lattice Thermal Conductivity in Thermoelectric Cu-Sb-Se Ternary Semiconductors. Phys. Rev. B 2012, 85 (5), 054306. https://doi.org/10.1103/PhysRevB.85.054306.

(60) Kohatsu, I.; Wuensch, B. J. The Crystal Structure of Aikinite, PbCuBiS$_3$. Acta Crystallogr. Sect. B 1971, 27 (6), 1245-1252. https://doi.org/10.1107/S0567740871003819.

(61) Ohmasa, M.; Nowacki, W. A redetermination of the crystal structure of aikinite [BiS$_2$/S/CuIVPbVII]. Z. Für Krist. - Cryst. Mater. 1970, 132 (1-6), 71-86. https://doi.org/10.1524/zkri.1970.132.16.71.
(62) Baur, W. H. The Geometry of Polyhedral Distortions. Predictive Relationships for the Phosphate Group. Acta Crystallogr. Sect. B 1974, 30 (5), 1195–1215.
https://doi.org/10.1107/S0567740874004560.

(63) Lundegaard, L. F.; Makovicky, E.; Boffa-Ballaran, T.; Balic-Zunic, T. Crystal Structure and Cation Lone Electron Pair Activity of Bi₂S₃ between 0 and 10 GPa. Phys. Chem. Miner. 2005, 32 (8), 578–584.
https://doi.org/10.1007/s00269-005-0033-2.

(64) Cordero, B.; Gómez, V.; Platero-Prats, A. E.; Revés, M.; Echeverría, J.; Cremades, E.; Barragán, F.; Alvarez, S. Covalent Radii Revisited. Dalton Trans. 2008, No. 21, 2832–2838.
https://doi.org/10.1039/B801115J.

(65) Shannon, R. D. Revised Effective Ionic Radii and Systematic Studies of Interatomic Distances in Halides and Chalcogenides. Acta Crystallogr. Sect. A 1976, 32 (5), 751–767.
https://doi.org/10.1107/S0567739476001551.

(66) Batsanov, S. S. Van Der Waals Radii of Elements. Inorg. Mater. 2001, 37 (9), 871–885.
https://doi.org/10.1023/A:1011625728803.

(67) Koumpouras, K.; Larsson, J. A. Distinguishing between Chemical Bonding and Physical Binding Using Electron Localization Function (ELF). J. Phys. Condens. Matter 2020, 32 (31), 315502.
https://doi.org/10.1088/1361-648X/ab7fd8.

(68) Remsing, R. C.; Klein, M. L. A New Perspective on Lone Pair Dynamics in Halide Perovskites. APL Mater. 2020, 8 (5), 050902.
https://doi.org/10.1063/5.0001908.

(69) Dutta, M.; Samanta, M.; Ghosh, T.; Voneshen, D. J.; Biswas, K. Evidence of Highly Anharmonic Soft Lattice Vibrations in a Zintl Rattler. Angew. Chem. 2021, 133 (8), 4305–4311.
https://doi.org/10.1002/ange.202013923.

(70) Luu, S. D. N.; Supka, A. R.; Nguyen, V. H.; Vo, D.-V. N.; T. Hung, N.; Wojciechowski, K. T.; Fornari, M.; Vaqueiro, P. Origin of Low Thermal Conductivity in In₄Se₃. ACS Appl. Energy Mater. 2020, 3 (12), 12549–12556.
https://doi.org/10.1021/acsaem.0c02489.

(71) Jana, M. K.; Pal, K.; Warankar, A.; Mandal, P.; Waghmare, U. V.; Biswas, K. Intrinsic Rattler-Induced Low Thermal Conductivity in Zintl Type TlInTe₂. J. Am. Chem. Soc. 2017, 139 (12), 4350–4353.
https://doi.org/10.1021/jacs.7b01434.

(72) Bhui, A.; Dutta, M.; Mukherjee, M.; Rana, K. S.; Singh, A. K.; Soni, A.; Biswas, K. Ultralow Thermal Conductivity in Earth-Abundant Cu₁.₆Bi₄.₈S₈: Anharmonic Rattling of Interstitial Cu. Chem. Mater. 2021, 33 (8), 2993–3001.
https://doi.org/10.1021/acs.chemmater.1c00659.

(73) Laurita, G.; Seshadri, R. Chemistry, Structure, and Function of Lone Pairs in Extended Solids. Acc. Chem. Res. 2022, 55 (7), 1004–1014.
https://doi.org/10.1021/acs.accounts.1c00741.

(74) Chmielowski, R.; Péré, D.; Bera, C.; Opahle, I.; Xie, W.; Jacob, S.; Capet, F.; Roussel, P.; Weidenkaff, A.; Madsen, G. K. H.; Dennler, G. Theoretical and Experimental Investigations of the Thermoelectric Properties of Bi₂S₃. J. Appl. Phys. 2015, 117 (12), 125103.
https://doi.org/10.1063/1.4916528.
(75) Hanus, R.; Agne, M. T.; Rettie, A. J. E.; Chen, Z.; Tan, G.; Chung, D. Y.; Kanatzidis, M. G.; Pei, Y.; Voorhees, P. W.; Snyder, G. J. Lattice Softening Significantly Reduces Thermal Conductivity and Leads to High Thermoelectric Efficiency. Adv. Mater. 2019, 31 (21), 1900108. https://doi.org/10.1002/adma.201900108.

(76) Agne, M. T.; Hanus, R.; Snyder, G. J. Minimum Thermal Conductivity in the Context of Diffusion-Mediated Thermal Transport. Energy Environ. Sci. 2018, 11 (3), 609–616. https://doi.org/10.1039/C7EE03256K.

(77) Lin, S.; Li, W.; Li, S.; Zhang, X.; Chen, Z.; Xu, Y.; Chen, Y.; Pei, Y. High Thermoelectric Performance of Ag<sub>0.6</sub>GaSe<sub>0.4</sub> Enabled by Low Cutoff Frequency of Acoustic Phonons. Joule 2017, 1 (4), 816–830. https://doi.org/10.1016/j.joule.2017.09.006.

(78) Balijapelly, S.; Hauble, A.; Sundaramoorthy, S.; Watts, J. L.; Kauzlarich, S. M.; Chernatynskiy, A.; Choudhury, A. Ultralow Lattice Thermal Conductivity in the Aikinite Structure Family, Cu<sub>x</sub>Pb<sub>x</sub>Bi<sub>2-x</sub>S<sub>3</sub>, and Thermoelectric Properties of Cu<sub>0.14</sub>Pb<sub>0.14</sub>Bi<sub>1.86</sub>S<sub>3</sub>. ACS Appl. Energy Mater. 2022, 5 (11), 14222–14230. https://doi.org/10.1021/acsaem.2c02790.

(79) Shiga, T.; Shiomi, J.; Ma, J.; Delaire, O.; Radzynski, T.; Lusakowski, A.; Esfarjani, K.; Chen, G. Microscopic Mechanism of Low Thermal Conductivity in Lead Telluride. Phys. Rev. B 2012, 85 (15), 155203. https://doi.org/10.1103/PhysRevB.85.155203.

(80) Voneshen, D. J.; Refson, K.; Borissenko, E.; Krisch, M.; Bosak, A.; Piovano, A.; Cemal, E.; Enderle, M.; Gutmann, M. J.; Hoesch, M.; Roger, M.; Gannon, L.; Boothroyd, A. T.; Uthayakumar, S.; Porter, D. G.; Goff, J. P. Suppression of Thermal Conductivity by Rattling Modes in Thermoelectric Sodium Cobaltate. Nat. Mater. 2013, 12 (11), 1028–1032. https://doi.org/10.1038/nmat3739.

(81) Nilsson, G.; Nelin, G. Phonon Dispersion Relations in Ge at 80 K. Phys. Rev. B 1971, 3 (2), 364–369. https://doi.org/10.1103/PhysRevB.3.364.

(82) Xie, H.; Hao, S.; Bao, J.; Slade, T. J.; Snyder, G. J.; Wolverton, C.; Kanatzidis, M. G. All-Inorganic Halide Perovskites as Potential Thermoelectric Materials: Dynamic Cation Off-Centering Induces Ultralow Thermal Conductivity. J. Am. Chem. Soc. 2020, 142 (20), 9553–9563. https://doi.org/10.1021/jacs.0c03427.

Table of Contents (TOC)

<figure><img src="image_14.png" /></figure>