# Ultralow lattice thermal transport and considerable wave-like phonon
tunneling in chalcogenide perovskite BaZrS3

Yu Wu<EMAIL_ADDRESS>Yangtze Delta Region Institute (Huzhou), University of
Electronic Science and Technology of China, Huzhou 313001, China  Ying Chen
School of Information Science and Technology, Fudan University, Shanghai
200433, China  Qiaoqiao Li  School of Physics, University of Electronic
Science and Technology of China, Chengdu 610054, China  Kui Xue  School of
Physics, University of Electronic Science and Technology of China, Chengdu
610054, China  Hezhu Shao  College of Electrical and Electronic Engineering,
Wenzhou University, Wenzhou, 325035, China  Hao Zhang<EMAIL_ADDRESS>School of Information Science and Technology, Fudan University, Shanghai
200433, China  Yiwu Research Institute of Fudan University, Chengbei Road,
Yiwu City, Zhejiang 322000, China  Liujiang Zhou<EMAIL_ADDRESS>Yangtze
Delta Region Institute (Huzhou), University of Electronic Science and
Technology of China, Huzhou 313001, China  School of Physics, University of
Electronic Science and Technology of China, Chengdu 610054, China

###### Abstract

Chalcogenide perovskites provide a promising avenue for non-toxic, stable
thermoelectric materials. Here, thermal transport and thermoelectric
properties of BaZrS3 as a typical orthorhombic perovskite are investigated. An
extremely low lattice thermal conductivity $\kappa_{L}$ of 1.84 W/mK at 300 K
is revealed for BaZrS3, due to the softening effect of Ba atoms on the lattice
and the strong anharmonicity caused by the twisted structure. We demonstrate
that coherence contributions to $\kappa_{L}$, arising from wave-like phonon
tunneling, leading to a 18 % thermal transport contribution at 300 K. The
increasing temperature softens the phonons, thus reducing the group velocity
of materials and increasing the scattering phase space. However, it
simultaneously reduces the anharmonicity, which is dominant in BaZrS3 and
ultimately improves the particle-like thermal transport. Further, by replacing
S atom with Se and Ti-alloying strategy, $ZT$ value of BaZrS3 is significantly
increased from 0.58 to 0.91 at 500 K, making it an important candidate for
thermoelectric applications.

## Introduction

The search for materials with low thermal conductivity is of great
significance in the field of thermoelectric conversion. Perovskite materials
typically exhibit low lattice thermal conductivity due to weak bonding and
strong anharmonicity1, 2. Among them, organic$-$inorganic halide perovskites
show great promise, but suffer from poor thermal and moisture stability3,
along with the toxicity of Pb4. The rise of halide perovskites encouraged
researchers to explore oxide and chalcogenide perovskites in the electronic
applications. Oxygen atoms have high electronegativity, leading to the large
electronegativity difference along the B$-$X bond in ABX3, which may enhance
the ionization and electron localization. Hence, oxide perovskites usually
have a wide band gap and poor electrical transport performance which is
unfavorable to the thermoelectric properties of the material3.

The lower electronegativity difference of the B$-$X bond in chalcogenide
perovskites makes them potential semiconductors for solar cell and other
electronic applications. BaZrS3 is a representative chalcogenide perovskite.
It is a semiconductor with direct band gap of $\sim 1.7-1.85\;\rm{eV}$ and
shows strong absorption in the visible spectrum5, 3. The band gap of BaZrS3
can be easily tuned under pressure or by atomic substitution6, 7. BaZrS3 is
also shown to be stable against oxidation and moisture8. Different from the
ideal cubic perovskite, BaZrS3 is isostructural with GdFeO3 and adopts the
twisted perovskite structure with orthorhombic phase9. Studies have shown that
the reduction of symmetry can significantly reduce the phonon degeneracy and
enhance the avoid-crossing behavior of phonon branches, which result in an
increase in the phonon$-$phonon scattering phase space to limit lattice
thermal transport10, 11, 12. Using T=0 K force constants, the lattice thermal
conductivity of BaZrS3 at 300 K is calculated to be as low as 1.16 W/mK base
on the iterative solution of Boltzmann Transport Equation (BTE)13. However, in
this scheme, the temperature effect only affects the Bose-Einstein
distribution of phonons, and its effects on phonon dispersion, phonon-phonon
interaction and electron-phonon interaction are completely ignored, which may
lead to misestimation of thermal conductivity14. In addition, the traditional
particle-like treatment of phonons may not be sufficient for describing
thermal transport properties in the complex BaZrS3 structure. Sun et al. found
large wave-like phonon tunneling effect through elastic and inelastic
scattering measurements in hexagonal chalcogenide perovskite BaTiS315. We
expect that this phenomenon will also occur in orthorhombic BaZrS3 and
contribute to the total lattice thermal conductivity. Without considering the
heat loss from lattice thermal transport, Osei-Agyemang et al. obtained an
upper limit of 1.0 for the $ZT$ value of BaZrS3 at 300 K16, indicating its
potential application in the thermoelectric field. The accurate lattice
thermal conductivity considering temperature and wave-like phonon tunneling
effects is necessary for the further reliable estimation of $ZT$ value.

Here, we unveil the lattice thermal transport and thermoelectric properties of
BaZrS3 by means of first principles calculations. The conventional particle-
like behavior and wave-like behavior of phonons are all taken into account in
lattice thermal transport calculations. We show that the temperature effect
can reduce the phonon group velocity, increase the phase space of phonon
scattering and reduce the anharmonicity of BaZrS3. Their influence on the
particle-like thermal conductivity $\kappa_{p}$ is competitive with each
other. Because the anharmonicity is more sensitive to temperature, the
conventional fixed force constants method underpredicts the $\kappa_{p}$. We
find significant contributions from wave-like transport channels in BaZrS3 due
to the closely spaced phonon branches. Further, we apply the lattice thermal
conductivity including wave-like thermal transport to thermoelectric
evaluation for the first time. By atom substitution method, we modulate the
thermoelectric performance of BaZrS3 and find that it is significantly
improved after replacing S atom with Se or Ti-alloying, indicating their
potential application in the thermoelectric field.

## Results and Discussion

The configurations of BaZrS3 are shown in Fig 1. BaZrS3 has a twisted
perovskite structure with a space group of $Pnma$. The framework of BaZrS3 is
composed of corner-shared ZrS6 octahedra with Zr atom at the center and S
atoms at the corners. The Ba atoms are located in the hollow within the
perovskite framework. The optimized lattice parameters are listed in Table 1
which agree well with previous results16, 9. The twisted structures are
usually accompanied by small group velocities, large phase space and lattice
anharmonicity10, which will together limit the thermal transport of materials.
From different directions, the ZrS6 octahedra in BaZrS3 have different
distortion forms. Along the $b$ axis, the adjacent octahedral layers have in-
phase rotation and are completely overlapping, as shown in Fig 1(a). The
distortion along [101] direction is out-of-phase with the adjacent octahedral
layers rotating in opposite directions, as shown in Fig 1(c). For in-phase
distortion, the twisted angle is the half of the difference between
$180^{\circ}$ and the dihedral angle of faces labeled in Fig 1(b). For out-of-
phase distortion, the twisted angle is the half of two Zr-Zr-S angles labeled
in Fig 1(d). The extracted in-phase and out-of-phase twisted angles are
$9.07^{\circ}$ and $6.98^{\circ}$, respectively, close to the results from the
experiment of $9.00^{\circ}$ and $7.03^{\circ}$9. There are two kinds of S
atoms in each unit cell of BaZrS3, denoted as S1 and S2, respectively. There
are two different Zr$-$S1 bonds, i.e. Zr$-$S1 and Zr$-$S${}_{1}^{*}$ with the
lengths of $2.541\;\rm{\AA}$ and $2.530\;\rm{\AA}$. The bond length between Zr
and S2 is relatively short, with a value of $2.524\;\rm{\AA}$. The bonding
between Ba and S plays an important role in the lattice distortion. In BaZrS3,
Ba atoms distort eight S atoms from the neighboring eight adjacent ZrS6
octahedra, with Ba-S bond lengths ranging from $3.160$ to $4.260\;\rm{\AA}$.
The complicated twisted angles and bond lengths indicate a poor phonon
transport properties of BaZrS3. Moreover, the second-order interaction-force
constants (2nd-order IFCs) of atomic pairs are extracted to reflect the
bonding strength as shown in Fig S1. The greater the IFC of a atomic pair, the
stronger the bonding strength. The Zr-S2 bonding is the strongest among these
atomic pairs, with the 2nd-order IFCs of $2.24\;\rm{eV\AA^{-2}}$. The Ba-S
bonding is much weaker than Zr-S2 bonding but comparable to Zr-S1 and
Zr-S${}_{1}^{*}$. The difference of 2nd-order IFCs between Zr-S and Ba-S pairs
is smaller than that between Sn-S and Ba-S pairs in BaSnS217. This is because
compared with Sn atoms, the electronegativity of Zr atoms is closer to that of
Ba atoms. Hence, the Ba atoms in BaZrS3 will have a more significant effect on
the distortion of the skeleton than in BaSnS2. The bonding between Ba atoms
and the Zr-S framework would significantly soften the lattice and limit the
phonon transport.

Figure 1: (a, b) Crystal structure of BaZrS3 viewed along $b$ axis with adjacent octahedral layers rotating in-phase. (c, d) Crystal structure of BaZrS3 viewed along [101] direction with adjacent octahedral layers rotating out-of-phase. Table 1: The calculated lattice parameters (Å), bong length (Å) and band gap (eV) of BaZrS3, compared with other calculated and experimental results. | $a$ | $b$ | $c$ | $\rm{Zr-S_{1}}$ | $\rm{Zr-S_{1}^{*}}$ | $\rm{Zr-S_{2}}$ |  $E_{g}$ (eV)
---|---|---|---|---|---|---|---
$this\,work$ | 7.084 | 9.944 | 6.960 | 2.541 | 2.530 | 2.524 | 1.79
$Cal.$ | 7.06616 | 9.92316 | 6.97216 | / | / | / | 1.7916
$Exp.$ | 7.0569 | 9.9629 | 6.9969 | 2.5399 | 2.5329 | 2.5289 | 1.70-1.857, 8, 5

The phonon dispersion and atom-resolved density of states for BaZrS3 at 100 K
and 500 K are presented in Fig 2. In BaZrS3, most phonon branches are much
flatter than those in ZrS3 and the maximum frequency of acoustic phonon modes
for BaZrS3 is 1.5 THz at 100 K much smaller than that of 3.9 THz in ZrS318. It
indicates that the phonon modes in BaZrS3 have lower group velocity due to the
softening effect of Ba atoms on the lattice. Moreover, the optical phonon
branches in the range of 2-4 THz show larger slopes than acoustic phonon
branches, which indicates higher phonon group velocities. The optical phonon
branches are close to each other. On one hand, it will increase the
phonon$-$phonon scattering phase space by increasing the scattering channels.
On the other hand, it gives optical phonons the opportunity to couple with
each other and benefits to the wave-like tunneling transport. From the atom-
resolved density of states, it can be seen that the acoustic and low frequency
optical branches are mainly contributed by Ba atoms, while the high frequency
optical branches are mainly contributed by Zr and S atoms. With the increase
of temperature, the phonon frequencies decrease overall and the density of
states tend to be compressed.

Figure 2:  (a) Calculated temperature-dependent phonon dispersions of BaZrS3
at 100 K and 500 K. (b) Atom-resolved phonon density of states of BaZrS3 at
100 K and 500 K.

In the framework of the Wigner formalism19, the lattice thermal conductivity
is generalized to an expression including both the populations’ and
coherences’ contributions20:

$\kappa_{L}^{\alpha\beta}=\kappa_{\mathrm{p}}^{\alpha\beta}+\kappa_{\mathrm{c}}^{\alpha\beta}$
(1)

where $\kappa_{\mathrm{p}}^{\alpha\beta}$ is related to the particle-like
propagation of phonon wavepackets discussed by Peierls’s semiclassical
picture, written as

$\kappa_{\mathrm{p}}^{\alpha\beta}=\frac{1}{VN}\sum_{\mathrm{q}\mathrm{s}}C_{\mathrm{q}}^{s}v_{\mathrm{q},\alpha}^{s}v_{\mathrm{q},\beta}^{s}\tau_{\mathrm{q}}^{s}$
(2)

$\kappa_{\mathrm{c}}^{\alpha\beta}$ is associated to the wave-like tunnelling
and loss of coherence between different phonon branches $s$ and $s^{\prime}$,
written as

$\displaystyle\kappa_{\mathrm{c}}^{\alpha\beta}=$
$\displaystyle\frac{\hbar^{2}}{k_{B}T^{2}VN}\sum_{\mathrm{q}}\sum_{s\neq
s^{\prime}}\frac{\omega_{\mathrm{q}}^{s}+\omega_{\mathrm{q}}^{s^{\prime}}}{2}v_{\mathrm{q},\alpha}^{s,s^{\prime}}v_{\mathrm{q},\beta}^{s,s^{\prime}}$
(3)
$\displaystyle\times\frac{\omega_{\mathrm{q}}^{s}n_{\mathrm{q}}^{s}\left(n_{\mathrm{q}}^{s}+1\right)+\omega_{\mathrm{q}}^{s^{\prime}}n_{\mathrm{q}}^{s^{\prime}}\left(n_{\mathrm{q}}^{s^{\prime}}+1\right)}{4\left(\omega_{\mathrm{q}}^{s}-\omega_{\mathrm{q}}^{s^{\prime}}\right)^{2}+\left(\Gamma_{\mathrm{q}}^{s}+\Gamma_{\mathrm{q}}^{s^{\prime}}\right)^{2}}\times\left(\Gamma_{\mathrm{q}}^{s}+\Gamma_{\mathrm{q}}^{s^{\prime}}\right)$

In these expressions, $\alpha$ and $\beta$ are the Cartesian indices, $V$ is
the the volume of the unit cell, $N$ is the number of sampled phonons in the
first Brillouin zone. $C_{\mathrm{q}}^{s}$, $v_{\mathrm{q}}^{s}$,
$\tau_{\mathrm{q}}^{s}$, $\omega_{\mathrm{q}}^{s}$,
$v_{\mathrm{q}}^{s,s^{\prime}}$, $\Gamma_{\mathrm{q}}^{s}$ are the heat
capacity, the group velocity, the lifetime, the phonon frequency, the
interband velocity matrix element and the scattering rate of a phonon mode
with wavevector $\mathrm{q}$. $n_{\mathrm{q}}^{s}$ is the equilibrium
Bose–Einstein distribution.

Fig 3(a) shows the anisotropic $\kappa_{\mathrm{p}}$ of BaZrS3 using the fixed
force fonstants (FFC) at 100 K and temperature dependent force constants
(TDFC) considering three-phonon scattering. The calculated
$\kappa_{\mathrm{p}}$ is 1.25 W/mK, 1.28 W/mK and 1.18 W/mK along the $a-$,
$b-$ and $c-$axis at 300 K using FFC method, which are close to the previous
results of Eric $et\,al$.13. The decrease of thermal conductivity with
temperature becomes slow when considering TDFC, making $\kappa_{\mathrm{p}}$
increase to 1.47 W/mK, 1.64 W/mK and 1.41 W/mK in the corresponding
directions. This phenomenon has also been reported in UO2 21, but the opposite
trend was found in LaWN3 recently2, which we will discuss later. As shown in
Fig 3(b), $\kappa_{\mathrm{c}}$ increases with the increasing temperature. The
$\kappa_{\mathrm{c}}$ is largest along the $b$ direction, with 0.38 W/mK at
300 K using TDFC. It can be seen in Fig 2 that the phonon dispersion is denser
and flatter along $\Gamma$-Z direction (related to $b$ direction), making it
easy for two phonons coupling with each other. Fig 3(c) shows the average
thermal conductivity by calculating the arithmetic average along three axes.
Since $\kappa_{\mathrm{p}}$ and $\kappa_{\mathrm{c}}$ show an opposite trend
with temperature, the contribution of $\kappa_{\mathrm{c}}$ to $\kappa_{L}$
gradually increases with the increase of temperature. At 300 K, $\kappa_{L}$
is 1.84 W/mK and $\kappa_{\mathrm{c}}$ reaches 18% of the $\kappa_{L}$.
Moreover, $\kappa_{\mathrm{c}}$ makes the temperature dependence of thermal
conductivity weaker, easing the decay from $T^{-0.82}$ to $T^{-0.59}$. Fig
3(c) shows the thermal conductivity of BaZrS3 compared with some other
inorganic perovskites. Our results show that the low thermal conductivity of
BaZrS3 has a competitive advantage, indicating its potential application in
the field of thermoelectric. The lattice thermal conductivity can be reduced
by reducing the symmetry of the material or by introducing heavy atoms22. On
the basis of BaZrS3 structure, by atomic substitution method, three materials
of BaZrSe3, BaZr0.75Ti0.25S3 and BaZr0.5Ti0.5S3 are obtained and their lattice
thermal conductivity are calculated. Compared with BaZrS3, they all have lower
lattice thermal conductivity. For example, by replacing the S atom with the
heavier Se, the lattice thermal conductivity of the obtained BaZrSe3 can
decrease to 1.48 W/mK.

Figure 3:  Calculated temperature-dependent (a) $\kappa_{\mathrm{p}}$ and (b)
$\kappa_{\mathrm{c}}$ using FFC and TDFC considering three-phonon interaction.
The result from reference is also shown with green line13. (c) The
temperature-dependent average $\kappa_{\mathrm{p}}$, $\kappa_{\mathrm{c}}$ and
$\kappa_{L}$. (d) $\kappa_{L}$ from other literature of some other inorganic
perovskites at 300 K.23, 2, 24, 15 $\kappa_{L}$ of LaWN3, BaZrS3, BaZrSe3,
BaZr0.75Ti0.25S3 and BaZr0.5Ti0.5S3 contains $\kappa_{\mathrm{c}}$
contribution, which is represented by the part above the solid red line. The
prefixes $c-$, $o-$, $t-$ and $h-$ represent cubic, orthorhombic, trigonal and
hexagonal, respectively. The resolved $\kappa_{\mathrm{c}}$ with various pairs
of phonon frequencies ($\omega^{s}$ and $\omega^{s^{\prime}}$) at 500 K using
(e) FFC and (f) TDFC methods.

In Fig 3(e, f), we depicted the $\kappa_{\mathrm{c}}$ terms resolved by
various pairs of phonon frequencies at 500 K. The phonon pairs around the
diagonal ($\omega^{s}=\omega^{s^{\prime}}$) have a larger contribution to
$\kappa_{\mathrm{c}}$. Compared with FFC method, the $\kappa_{\mathrm{c}}$ is
more concentrated when considering TDFC due to the softening of phonon
dispersion with the increase of temperature.

To investigate the contribution of phonons with different frequencies to
$\kappa_{\mathrm{p}}$, Fig 4(a) shows the frequency dependent cumulative
$\kappa_{\mathrm{p}}$ and the differential $\kappa_{\mathrm{p}}$ at 300 K. The
most significant increase in thermal conductivity is observed in two regions
with the frequency ranges of 0-2 THz and 2-4 THz. The former mainly contains
acoustic phonons. The latter related to low-frequency optical phonons which
contributes 42% to the total $\kappa_{\mathrm{p}}$ along the $a$-axis. As seen
in Fig 4(b), the contribution of optical phonons to the average
$\kappa_{\mathrm{p}}$ maintains around 75% in the temperature ranging from 100
K to 500 K, close to 68% in BaSnS2 and 71% in SnS17. Atomic vibrations of
$\Gamma$-point phonon modes corresponding to V1, V2 and V3 in Fig 4(a) are
shown in Fig 4(c). For the V1 mode, where Ba vibrations dominate in the $a-c$
plane, the neighboring Ba atoms along $b$ and [101] direction both vibrate in
opposite direction. The V2 mode is also mainly contributed by the vibration of
Ba atoms, but vibrating along $b$ direction. The neighboring Ba atoms along
$b$ direction vibrate in opposite direction but in the same direction along
[101] direction. The V3 mode is mainly contributed by S2 vibrations. The S2
atoms in one octahedra vibrate in the same direction but vibrate in a vertical
direction in the adjacent octahedra. The vibrations along different directions
lead to significant deformation for the lattice, which will further bring huge
anharmonicity.

Figure 4:  (a) Calculated cumulative $\kappa_{\mathrm{p}}$ and differential
$\kappa_{\mathrm{p}}$ as a function of phonon frequencies. (b) Temperature
dependence of $\kappa_{\mathrm{p}}$ contributed from three acoustic phonon
branches (TA1, TA2 and LA) and all the optical branches (OP). (c) Atomic
vibrations of $\Gamma$-point phonon modes corresponding to the peaks with the
most $\kappa_{\mathrm{p}}$ contribution.

To further understand the underlying mechanisms of the low lattice thermal
conductivities and the dominant contribution by optical phonons for BaZrS3, we
also calculate the corresponding group velocity $V_{g}$, phonon lifetime
$\tau$, and Grüneisen parameter $\gamma$. As shown in Fig 5(a), most of the
phonons have group velocity less than 2 km/s due to the flat band properties
seen in Fig 2(a). Moreover, the group velocity of optical phonons is generally
higher than that of acoustic phonons, especially for optical phonons in the
2-4 THz frequency range. These phonons with large group velocity will mainly
contribute to the thermal conductivity of the material. The average group
velocity can be defined as

$\bar{V_{g}}=\sqrt{\frac{\sum_{\mathrm{q}s}C_{{\mathrm{q}s}}v_{\mathrm{q}s}^{2}\tau_{\mathrm{q}s}N_{\mathrm{q}s}}{\sum_{\mathrm{q}s}C_{{\mathrm{q}s}}\tau_{\mathrm{q}s}N_{\mathrm{q}s}}},$
(4)

where $N$ is the symmetry degeneracy of corresponding mode $\mathrm{q}s$. The
calculated average group velocity is 0.87 km/s for BaZrS3, comparable to that
of low thermal conductivity material BaSnS2 with 0.8 km/s17. Fig 5(b) shows
the mode resolved phonon lifetime. The lifetime of acoustic phonons is in the
range of 1-100 ps, while is 0.1-10 ps for the optical phonons. Shorter phonon
lifetimes indicate greater scattering and weaker contribution to thermal
transport. Hence, optical phonons dominating thermal transport in BaZrS3 is
mainly due to their large group velocity. Fig S2 also shows the scattering
rates contributed by the interaction of three and four phonons. The four-
phonon scattering rate is more than two orders of magnitude smaller than the
three-phonon scattering rate in a wide frequency range, indicating that it is
reasonable to consider only the three-phonon interaction. The average lifetime
can be written as

$\bar{\tau}=\frac{\sum_{\mathrm{q}s}C_{\mathrm{q}s}v_{\mathrm{q}s}^{2}\tau_{\mathrm{q}s}}{\sum_{\mathrm{q}n}C_{\mathrm{q}s}v_{\mathrm{q}s}^{2}}.$
(5)

The calculated $\bar{\tau}$ is 1.26 ps for BaZrS3, lower than that in BaSnS2
with 2.37 ps17, indicating the potential lattice anharmonicity in BaZrS3.
Besides, we identify the signature of a wave-like tunneling behavior in the
heat carrying phonons. The horizontal lines in Fig 5(b) is the Wigner limit
$\tau_{Wigner}$ reflecting the average phonon interband spacing, i.e.,
$\tau_{Wigner}=N_{s}/\omega_{\max}$, where $N_{s}$ is the total number of
phonon bands and $\omega_{\max}$ is the maximum phonon frequency25. The green
curve gives the phonon lifetime in the Ioffe-Regel limit ($\tau_{Ioffe-
Regel}=1/\omega$)26. The majority of phonon lifetime are in the $\tau_{Ioffe-
Regel}<\tau<\tau_{Wigner}$ range resulting in large $\kappa_{\mathrm{c}}$
contributions. The Grüneisen parameter $\gamma$ reflects the anharmonic
property of the lattice determined by the change of phonon frequencies with
respect to the change of unit-cell volume. The higher the absolute value of
$\gamma$, the stronger the lattice anharmonicity. As seen in Fig 5(c),
acoustic phonons and low-frequency optical phonons have high $\gamma$ range
from -4 to 9. The average Grüneisen parameter is defined as

$\bar{\gamma}=\frac{\sum_{\mathrm{q}s}\left|\gamma_{\mathrm{q}s}\right|C_{{\mathrm{q}s}}}{\sum_{\mathrm{q}s}C_{{\mathrm{q}s}}}$
(6)

The calculated $\bar{\gamma}$ is 1.34 for BaZrS3 comparable to some materials
with strong anharmonicity such as PbTe ($\gamma\sim 1.40$)27 and BiCuSeO
($\gamma\sim 1.50$)28.

To understant where differences in thermal transport originate for TDFC and
FFC calculations, we compare the $\bar{V_{g}}$, $\bar{\tau}$ and
$\bar{\gamma}$ at different temperatures derived from the two methods. As
shown in Fig 5(d), the $\bar{V_{g}}$ using FFC is around 0.9 km/s but
gradually decreases as the temperature rises when using TDFC. This result is
consistent with that the phonon dispersion is softened with the increase of
temperature as shown in Fig 2(a). The change rate of $\bar{V_{g}}$ is -3.7%
after using TDFC at 300 K. This change in $\bar{V_{g}}$ will limit the thermal
transport performance of the material. As the temperature increases, the
$\bar{\tau}$ decreases gradually for both methods, but at each temperature,
the results obtained by TDFC method are larger. The change rate of
$\bar{\tau}$ is 28.4% after using TDFC at 300 K, which will improve the
thermal conductivity. Hence, the differences in thermal transport using TDFC
and FFC calculations lie in the change of phonon lifetime. In fact, the
lifetime of phonons are determined by the strength of anharmonic scatterings
and the scattering possibilities. The former is related to the anharmonicity
and can be characterized by $\gamma$. The latter is determined by the three-
phonon scatterings which include emission and absorption processes and can be
characterized by the phase space $P_{3}$. Fig S3 shows $P_{3}$ at 300 K using
FFC and TDFC methods. It can be seen that the temperature effect only causes a
weak enhancement of the scattering channel in the low frequency region. As
shown in Fig 5(f), the $\bar{\gamma}$ using FFC method remains almost constant
($\sim 1.57$), while generally decreases as the temperature rises using TDFC.
It indicates that the temperature effect can lower the anharmonicity of
BaZrS3, which further increases the phonon lifetime at a specific temperature
and thus delays the trend of thermal conductivity decreasing with increasing
temperature shown in Fig 3(a).

Figure 5: Mode-resolved (a) group velocity (b) lifetime, and (c) Grüneisen
parameter versus phonon frequency in BaZrS3 at 300 K using TDFC. (d-f)
Temperature dependent average group velocity, lifetime, and Grüneisen
parameter using TDFC and FFC.

To precisely evaluate the electronic transport properties of BaZrS3, ionized
impurity (IMP) scattering, acoustic deformation-potential (ADP) scattering and
polar optical phonon (POP) scattering are taken into considerations. Fig S4(a,
b) show the electron and hole mobility decomposed by various scattering
mechanisms with carrier concentration $n=10^{19}\;\rm{cm^{-3}}$. BaZrS3 is a
polar perovskite material, in which there is a large dielectric electron-
phonon interaction, yet the kinetic energy of the electron is sufficient to
keep the polaron from collapsing into a localized state29. In this material
the polar-optical phonon mode scattering will dominate and thus limit
mobility. In BaZrS3, the calculated electron and hole mobility are
$37.9\;\rm{cm^{2}/V\cdot s}$ and $4.3\;\rm{cm^{2}/V\cdot s}$ at 300 K,
respectively, close to that of CH3NH3PbI3 $\sim 30\;\rm{cm^{2}/V\cdot s}$30.
Fig S4(c, d) show the electron and hole scattering rates for IMP, ADP and POP
scattering. The POP scattering rates of most electrons and holes are $\sim
10^{14}\;\rm{s^{-1}}$, which is 1-2 orders of magnitude higher than the
scattering rates of IMP and ADP. Our previous work suggests that Se-alloying
strategy can weaken the electron$-$phonon coupling and prolong the
nonradiative electron$-$hole recombination lifetime in BaZrS3 while the
opposite is true when using the Ti-alloying strategy31. The calculated
electron and hole mobility are $49.5\;\rm{cm^{2}/V\cdot s}$ and
$17.0\;\rm{cm^{2}/V\cdot s}$ at 300 K for BaZrSe3, which confirms the superior
electrical properties of BaZrSe3 over BaZrS3. The electron mobility of
BaZr0.75Ti0.25S3 and BaZr0.5Ti0.5S3 decrease by about half compared with
BaZrS3 as shown in Fig S5(a).

Figure 6: The thermoelectric properties for BaZrS3. (a) seebeck coefficients,
(b) electronic conductivities and (c) $ZT$ values vesus carrier concentration
at 300 K. (d) The maximum $ZT$ at different temperatures.

The thermoelectric performance of a material is quantified by a dimensionless
figure of merit $ZT$ ($ZT=S^{2}\sigma T/(\kappa_{e}+\kappa_{L})$), where $S$
represents the seebeck coefficient, $\sigma$ is the electron conductivity, and
$\kappa_{e/L}$ is the electronic/lattice thermal conductivity. The
thermoelectric parameters of seebeck coefficients $S$ and electronic
conductivities $\sigma$ varying with carrier concentration $n$ at 300 K are
shown in Fig 6. The band- and momentum-resolved relaxation time are considered
in the calculation including ADP and POP scattering. The $S$ is almost
isotropic and the average value is predicted to be $389\;\rm{\mu V/K}$ and
$597\;\rm{\mu V/K}$ at $n=10^{19}\;\rm{cm^{-3}}$ for n-type and p-type system,
respectively. The $S$ for BaZrS3 is comparable to some traditional
thermoelectric materials including Bi2Te3 ($215\;\rm{\mu V/K}$)32, PbTe
($185\;\rm{\mu V/K}$)33, and SnSe ($510\;\rm{\mu V/K}$)34. The large $S$ for
BaZrS3 can be understood by the Mott relation35 that large seebeck
coefficients tend to appear on the edge of the band with large DOS. Fig S6
shows the band structure and DOS of BaZrS3. Both the valence band and the
conduction band of BaZrS3 have multi-valley features, leading large DOS on
both sides of the Fermi level. $\sigma$ increases gradually with the increase
of carrier concentration. At $n=10^{19}\;\rm{cm^{-3}}$, the average
$\sigma=0.20\;\rm{kS/m}$ under p-type doping is significantly lower than that
of competing electronic materials. For n-type system, $\sigma$ has a higher
magnitude of $\sigma=1.92\;\rm{kS/m}$ because of the lower carrier effective
mass and scattering rate for electrons compared with holes. $\kappa_{e}$ can
be calculated according to the Wiedemann-Franz law, i.e. $\kappa_{e}=L\sigma
T$, where $L$ is the Lorenz number. The dimensionless figure of merit $ZT$
vesus carrier concentration at 300 K for BaZrS3 are calculated and shown in
Fig 6(c). A maximum value of $ZT=0.19$ at $n=2.6\times 10^{21}\;\rm{cm^{-3}}$
and $ZT=0.23$ at $n=4.1\times 10^{20}\;\rm{cm^{-3}}$ is realized for p-type
and n-type system, respectively. As the temperature increases, the maximum
$ZT$ increases further up to 0.58 at 500 K mainly due to the decrease of
lattice thermal conductivity as shown in Fig 6(d).

The band structure of BaZrSe3, BaZr0.75Ti0.25S3 and BaZr0.5Ti0.5S3 are shown
in Fig S7. Compared with BaZrS3, their band gap decreases significantly and
the multi-valley features are inherited, which are favorable for
thermoelectric performance. Fig S5(b) shows the seebeck coefficients when $ZT$
reaches its maximum value. It can be seen that, compared with BaZrS3, the
electron seebeck coefficients are significantly increased for both BaZrSe3,
BaZr0.75Ti0.25S3 and BaZr0.5Ti0.5S3. Due to lower lattice thermal conductivity
and better electrical transport performance, after replacing S element with Se
element in BaZrS3, the maximum $ZT$ of n-type BaZrSe3 can reach 0.42 at 300 K,
which is much higher than the inorganic perovskite CsPbI3 (0.08) and
comparable to the hybrid perovskite $\delta$-FAPbI3 (0.45)36. The $ZT$ can be
further increased to 0.91 at 500 K.

## Conclusion

In conclusion, we investigate the lattice thermal transport and thermoelectric
properties of BaZrS3 through first-principles calculations. BaZrS3 have low
lattice thermal conductivity $\kappa_{L}$ with 1.84 W/mK at 300 K due to weak
bonding and strong anharmonicity. The twisted structure and bonding between Ba
atoms and the Zr-S framework significantly soften the lattice and lower the
phonon group velocity. Because of the relatively high group velocity of
optical phonons, they make a major contribution to thermal transport. We
investigate how the temperature effect influence the thermal transport. The
temperature effect mainly reduces the anharmonicity of the material and thus
reduces the decrease of particle-like thermal conductivity with temperature.
Our work also highlights the importance of wave-like thermal transport
contributions to $\kappa_{L}$, especially at high temperatures. The polar
optical phonons dominate the carrier scattering in BaZrS3 and the electron
mobility is $37.9\;\rm{cm^{2}/V\cdot s}$ at 300 K with
$n=10^{19}\;\rm{cm^{-3}}$. BaZrS3 has high seebeck coefficient due to multi-
valley band structure. At 500 K, the maximum $ZT$ value is 0.58 for BaZrS3 and
can be further improved by replacing the S element with Se or Ti-alloying
strategy, indicating excellent thermoelectric performance in the perovskite
family.

## Numerical methods

Calculations relating to electrical properties are implemented using the
Vienna Ab Initio simulation package (VASP) based on density functional theory
(DFT)37. The revised PBE-GGA for solids (PBEsol)38 is evaluated as the
exchange-correlation functional with a kinetic energy cutoff of
$500\;\rm{eV}$. The energy convergence value between two consecutive steps is
set as $10^{-5}\;$eV when optimizing atomic positions and the maximum
Hellmann-Feynman (HF) force acting on each atom is $10^{-3}\;$eV/Å. The
reciprocal space is sampled by a grid of $7\times 5\times 7$ k points in the
Brillouin zone. For the band gap calculations, the hybrid HSE06 functional39
is used with the screening parameter setting as 0.25 Å-1. The scattering rates
and carrier mobilities are calculated using the AMSET package40 considering
ionized impurity (IMP) scattering, acoustic deformation-potential (ADP)
scattering and polar optical (POP) scattering mechanism. The obtained band-
and momentum-relevant scattering rates are then used as inputs of the
BoltzTrap2 package41 to evaluate the thermoelectric properties of the
material.

The calculations of $\kappa$ and other relevent parameters such as phonon
relaxation time are carried out by the ShengBTE software42 which operates
based on the iterative scheme. The code used to calculate wave-like thermal
conductivity has been integrated into ShengBTE software by us. The
$\bf{q}$-mesh in the first irreducible Brillouin Zone is set be
13$\times$11$\times$13, the scale parameter for Gaussian smearing is set to be
0.1, and the supercells of 2$\times$2$\times$2 containing 160 atoms are
chosen. The 2nd- and 3rd-force constants are extracted from the temperature
dependent effective potential (TDEP)43, which represent the phonon-phonon
interactions and electron-phonon interactions at current temperatures as
accurate as possible with a finite order. The results of ab initio molecular
dynamics simulations with time interval of 1 fs and total simulation steps of
4000 are set as the training data of DeePMD-kit44 to train for the machine
learning potential (MTP). Then, the obtained MTP file will be the input file
for the LAMMPS package45. At different temperatures, the canonical (NVT)
ensemble was first employed for running 1 ns with a time step of 1 fs to bring
the system to thermal equilibrium. Then the system was sampled by 0.1 ps
within a period time of 100 ps. These samples were used as the inputs of TDEP
to compute the 2nd(3rd)-order IFCs. The truncation radii for 2nd(3rd)-order
IFCs are set to be 6.94Å and 13.43Å, respectively.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

This work is supported by Natural Science Foundation of China (12304038),
Huzhou Natural Science Foundation (2023YZ50), Shanghai Municipal Natural
Science Foundation under Grant (19ZR1402900), Sichuan Natural Science
Foundation (2022NSFSC1193), the Startup funds of Outstanding Talents of UESTC
(A1098531023601205), National Youth Talents Plan of China (G05QNQR049).

## Reference

## References

  * Pandey et al. 2022 Pandey, T.; Du, M.-H.; Parker, D. S.; Lindsay, L. Origin of ultralow phonon transport and strong anharmonicity in lead-free halide perovskites. _Materials Today Physics_ 2022, _28_ , 100881
  * Tong et al. 2023 Tong, Z.; Zhang, Y.; Pecchia, A.; Yam, C.; Zhou, L.; Dumitrică, T.; Frauenheim, T. Predicting the lattice thermal conductivity in nitride perovskite LaWN3 from ab initio lattice dynamics. _Advanced Science_ 2023, _10_ , 2205934
  * Swarnkar et al. 2019 Swarnkar, A.; Mir, W. J.; Chakraborty, R.; Jagadeeswararao, M.; Sheikh, T.; Nag, A. Are chalcogenide perovskites an emerging class of semiconductors for optoelectronic properties and solar cell? _Chemistry of Materials_ 2019, _31_ , 565–575
  * Babayigit et al. 2016 Babayigit, A.; Ethirajan, A.; Muller, M.; Conings, B. Toxicity of organometal halide perovskite solar cells. _Nature Materials_ 2016, _15_ , 247–251
  * Meng et al. 2016 Meng, W.; Saparov, B.; Hong, F.; Wang, J.; Mitzi, D. B.; Yan, Y. Alloying and defect control within chalcogenide perovskites for optimized photovoltaic application. _Chemistry of Materials_ 2016, _28_ , 821–829
  * Gross et al. 2017 Gross, N.; Sun, Y.-Y.; Perera, S.; Hui, H.; Wei, X.; Zhang, S.; Zeng, H.; Weinstein, B. Stability and Band-Gap tuning of the chalcogenide perovskite BaZrS3 in Raman and optical investigations at high pressures. _Physical Review Applied_ 2017, _8_ , 044014
  * Niu et al. 2016 Niu, S.; Huyan, H.; Liu, Y.; Yeung, M.; Ye, K.; Blankemeier, L.; Orvis, T.; Sarkar, D.; Singh, D. J.; Kapadia, R.; Ravichandran, J. Bandgap control via structural and chemical tuning of transition metal perovskite chalcogenides. _Advanced Materials_ 2016, _29_ , 1604733
  * Perera et al. 2016 Perera, S.; Hui, H.; Zhao, C.; Xue, H.; Sun, F.; Deng, C.; Gross, N.; Milleville, C.; Xu, X.; Watson, D. F.; Weinstein, B.; Sun, Y.-Y.; Zhang, S.; Zeng, H. Chalcogenide perovskites – an emerging class of ionic semiconductors. _Nano Energy_ 2016, _22_ , 129–135
  * Niu et al. 2019 Niu, S.; Zhao, B.; Ye, K.; Bianco, E.; Zhou, J.; McConney, M. E.; Settens, C.; Haiges, R.; Jaramillo, R.; Ravichandran, J. Crystal growth and structural analysis of perovskite chalcogenide BaZrS3 and Ruddlesden–Popper phase Ba3Zr2S7. _Journal of Materials Research_ 2019, _34_ , 3819–3826
  * Liu et al. 2022 Liu, C.; Wu, C.; Song, T.; Zhao, Y.; Yang, J.; Lu, P.; Zhang, G.; Chen, Y. An efficient strategy for searching high lattice thermal conductivity materials. _ACS Applied Energy Materials_ 2022, _5_ , 15356–15364
  * Liu et al. 2022 Liu, C.; Chen, Z.; Wu, C.; Qi, J.; Hao, M.; Lu, P.; Chen, Y. Large thermal conductivity switching in ferroelectrics by electric Field-Triggered crystal symmetry Engineering. _ACS Applied Materials & Interfaces_ 2022, _14_ , 46716–46725
  * Wu et al. 2023 Wu, Y.; Chen, Y.; Peng, L.; Zhang, H.; Zhou, L. Intervalley scattering induced significant reduction in lattice thermal conductivities for phosphorene. _Nanoscale Horizons_ 2023, _8_ , 912–920
  * Osei-Agyemang and Balasubramanian 2019 Osei-Agyemang, E.; Balasubramanian, G. Understanding the extremely poor lattice thermal transport in chalcogenide perovskite BaZrS3. _ACS Applied Energy Materials_ 2019, _3_ , 1139–1144
  * Shao et al. 2023 Shao, H.; Ding, D.; Fang, Y.; Song, W.; Huang, J.; Dong, C.; Zhang, H. Phonon transport in Cu2GeSe 3 : Effects of spin-orbit coupling and higher-order phonon-phonon scattering. _Physical Review B_ 2023, _107_ , 085202
  * Sun et al. 2020 Sun, B. et al. High frequency atomic tunneling yields ultralow and glass-like thermal conductivity in chalcogenide single crystals. _Nature Communications_ 2020, _11_ , 6039
  * Osei-Agyemang et al. 2021 Osei-Agyemang, E.; Koratkar, N.; Balasubramanian, G. Examining the electron transport in chalcogenide perovskite BaZrS3. _Journal of Materials Chemistry C_ 2021, _9_ , 3892–3900
  * Li et al. 2021 Li, Z.; Xie, H.; Hao, S.; Xia, Y.; Su, X.; Kanatzidis, M. G.; Wolverton, C.; Tang, X. Optical phonon dominated heat transport: a first-principles thermal conductivity study of BaSnS2. _Physical Review B_ 2021, _104_ , 245209
  * Wang et al. 2020 Wang, C.; Zheng, C.; Gao, G. Bulk and monolayer ZrS3 as promising anisotropic thermoelectric materials: a comparative study. _The Journal of Physical Chemistry C_ 2020, _124_ , 6536–6543
  * Moyal 1949 Moyal, J. E. Quantum mechanics as a statistical theory. _Mathematical Proceedings of the Cambridge Philosophical Society_ 1949, _45_ , 99–124
  * Simoncelli et al. 2019 Simoncelli, M.; Marzari, N.; Mauri, F. Unified theory of thermal transport in crystals and glasses. _Nature Physics_ 2019, _15_ , 809–813
  * Yang et al. 2022 Yang, X.; Tiwari, J.; Feng, T. Reduced anharmonic phonon scattering cross-section slows the decrease of thermal conductivity with temperature. _Materials Today Physics_ 2022, _24_ , 100689
  * Li et al. 2018 Li, J.; Zhang, X.; Chen, Z.; Lin, S.; Li, W.; Shen, J.; Witting, I. T.; Faghaninia, A.; Chen, Y.; Jain, A.; Chen, L.; Snyder, G. J.; Pei, Y. Low-Symmetry rhombohedral GeTe thermoelectrics. _Joule_ 2018, _2_ , 976–987
  * Zhao et al. 2021 Zhao, Y.; Zeng, S.; Li, G.; Lian, C.; Dai, Z.; Meng, S.; Ni, J. Lattice thermal conductivity including phonon frequency shifts and scattering rates induced by quartic anharmonicity in cubic oxide and fluoride perovskites. _Physical Review B_ 2021, _104_ , 224304
  * Xie et al. 2020 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. _Journal of the American Chemical Society_ 2020, _142_ , 9553–9563
  * Simoncelli et al. 2022 Simoncelli, M.; Marzari, N.; Mauri, F. Wigner formulation of thermal transport in solids. _Physical Review X_ 2022, _12_ , 041011
  * Allen and Feldman 1993 Allen, P. B.; Feldman, J. L. Thermal conductivity of disordered harmonic solids. _Physical Review B_ 1993, _48_ , 12581–12588
  * Morelli et al. 2008 Morelli, D. T.; Jovovic, V.; Heremans, J. P. Intrinsically minimal thermal conductivity in cubic I$-$V$-$VI2 semiconductors. _Physical Review Letters_ 2008, _101_ , 035901
  * Pei et al. 2013 Pei, Y.-L.; He, J.; Li, J.-F.; Li, F.; Liu, Q.; Pan, W.; Barreteau, C.; Berardan, D.; Dragoe, N.; Zhao, L.-D. High thermoelectric performance of oxyselenides: intrinsically low thermal conductivity of Ca-doped BiCuSeO. _NPG Asia Materials_ 2013, _5_ , e47
  * Frost 2017 Frost, J. M. Calculating polaron mobility in halide perovskites. _Physical Review B_ 2017, _96_ , 195202
  * Herz 2017 Herz, L. M. Charge-Carrier mobilities in metal halide perovskites: fundamental mechanisms and limits. _ACS Energy Letters_ 2017, _2_ , 1539–1548
  * Li et al. 2022 Li, Q.; Yan, L.; Chu, W.; He, J.; Luo, H.; Frauenheim, T.; Tretiak, S.; Zhou, L. Control of polaronic behavior and carrier lifetimes via metal and anion alloying in chalcogenide perovskites. _The Journal of Physical Chemistry Letters_ 2022, _13_ , 4955–4962
  * Poudel et al. 2008 Poudel, B.; Hao, Q.; Ma, Y.; Lan, Y.; Minnich, A.; Yu, B.; Yan, X.; Wang, D.; Muto, A.; Vashaee, D.; Chen, X.; Liu, J.; Dresselhaus, M. S.; Chen, G.; Ren, Z. High-Thermoelectric performance of nanostructured bismuth antimony telluride bulk alloys. _Science_ 2008, _320_ , 634–638
  * Nielsen et al. 2012 Nielsen, M. D.; Levin, E. M.; Jaworski, C. M.; Schmidt-Rohr, K.; Heremans, J. P. Chromium as resonant donor impurity in PbTe. _Physical Review B_ 2012, _85_ , 045210
  * Zhang and Zhao 2015 Zhang, X.; Zhao, L.-D. Thermoelectric materials: Energy conversion between heat and electricity. _Journal of Materiomics_ 2015, _1_ , 92–105
  * Sun et al. 2015 Sun, P.; Wei, B.; Zhang, J.; Tomczak, J. M.; Strydom, A.; Søndergaard, M.; Iversen, B. B.; Steglich, F. Large Seebeck effect by charge-mobility engineering. _Nature Communications_ 2015, _6_ , 7475
  * Shi et al. 2023 Shi, W.; Yao, M.; Wu, X.; Zhou, T.; Yong, X.; Deng, T.; Ma, H.; Xi, J. Atomistic Insights into the Origin of High-Performance thermoelectric response in hybrid perovskites. _Advanced Science_ 2023, _10_ , 2300666
  * Kresse and Furthmuller 1996 Kresse, G.; Furthmuller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Phys. Rev. B_ 1996, _54_ , 11169–11186
  * Perdew et al. 2008 Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring the Density-Gradient expansion for exchange in solids and surfaces. _Physical Review Letters_ 2008, _100_ , 136406
  * Heyd et al. 2003 Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. _The Journal of Chemical Physics_ 2003, _118_ , 8207–8215
  * Ganose et al. 2021 Ganose, A. M.; Park, J.; Faghaninia, A.; Woods-Robinson, R.; Persson, K. A.; Jain, A. Efficient calculation of carrier scattering rates from first principles. _Nature Communications_ 2021, _12_ , 2222
  * Madsen et al. 2018 Madsen, G. K.; Carrete, J.; Verstraete, M. J. BoltzTraP2, a program for interpolating band structures and calculating semi-classical transport coefficients. _Computer Physics Communications_ 2018, _231_ , 140–145
  * Li et al. 2014 Li, W.; Carrete, J.; Katcho, N. A.; Mingo, N. ShengBTE: A solver of the Boltzmann transport equation for phonons. _Computer Physics Communications_ 2014, _185_ , 1747–1758
  * Hellman and Abrikosov 2013 Hellman, O.; Abrikosov, I. A. Temperature-dependent effective third-order interatomic force constants from first principles. _Physical Review B_ 2013, _88_ , 144301
  * Wang et al. 2018 Wang, H.; Zhang, L.; Han, J.; E, W. DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics. _Computer Physics Communications_ 2018, _228_ , 178–184
  * Thompson et al. 2022 Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; in 't Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D.; Shan, R.; Stevens, M. J.; Tranchida, J.; Trott, C.; Plimpton, S. J. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. _Computer Physics Communications_ 2022, _271_ , 108171