Thermoelectric properties of the bismuth oxychalcogenides Bi2SO2, Bi2SeO2 and Bi2TeO2

We present a detailed theoretical study of the thermoelectric properties of the bismuth oxychalcogenides Bi2ChO2 (Ch = S, Se, Te). The electrical transport is modelled using semi-classical Boltzmann transport theory with electronic structures from hybrid density-functional theory, including an approximate model for the electron lifetimes. The lattice thermal conductivity is calculated using first-principles phonon calculations with an explicit treatment of anharmonicity, yielding microscopic insight into how partial replacement of the chalcogen in the bismuth chalcogenides impacts the phonon transport. We find very good agreement between the predicted transport properties and a favourable cancellation of errors that allows for near-quantitative predictions of the thermoelectric figure of merit ZT. Our calculations suggest recent experiments on n-doped Bi2SeO2 have achieved close to the largest ZT possible in bulk materials, whereas the largest reported ZT for Bi2TeO2 could be improved sixfold by optimising the carrier concentration. We also predict that much larger ZT > 2.5, competitive with the benchmark thermoelectric SnSe, could be obtained for Bi2SO2 and Bi2SeO2 with heavy p-type doping. This study demonstrates the predictive power of this modelling approach for studying thermoelectrics and highlights several avenues for improving the performance of the Bi2ChO2.


Introduction
Reducing anthroprogenic greenhouse gas emissions and mitigating climate change while meeting the ever-increasing demand for energy are among the most important challenges in contemporary science.Achieving these goals requires a range of technologies including improved clean energy sources and supporting energy-storage systems as well as methods for improving the efficiency of energy-intensive processes.The latter is particularly important given that 60% of global energy is currently wasted as heat, and this is projected to remain above 50% into the next decade [1,2].
One of the leading contenders for improving energy efficiency is thermometric (TE) power, which harnesses the Seebeck effect in a TE material to recover waste heat as electrical energy [3].Thermoelectric generators (TEGs) are solid-state devices with no moving parts and have diverse applications ranging from established uses in the aerospace industry to prospective uses for powering remote-sensing devices, recovering heat from combustion engine exhaust gases, and repurposing decommissioned oil rigs as geothermal power plants [3,4].
The performance of a TE material is typically expressed by the dimensionless figure of merit ZT: [1,3] ZT = S 2 σT κ el + κ latt (1) where S is the Seebeck coefficient, σ is the electrical conductivity, and κ el and κ latt are the electronic and lattice (phonon) components of the thermal conductivity κ. S, σ and κ el are properties of the electronic structure and are usually tightly coupled through the carrier concentration such that the optimum balance is .Atoms colours: Bi-purple; Ch-gold; O-red.This image was created using the VESTA software [36].
achieved in heavily-doped semiconductors [1].On the other hand, the κ latt depends on the structure and chemical bonding and is generally minimised in materials formed from heavy elements, with weak chemical bonding [5], and/or with structural features such as phase transitions and active lone pairs that lead to strongly-anharmonic lattice dynamics [6][7][8][9].These requirements have historically led to a focus on 'heavy' chalcogenide materials such as the current industry standards for low-and high-temperature applications, viz.Bi 2 Te 3 (ZT ≃ 1 from 350-450 K) [3] and PbTe (ZT ≃ 2.2 with endotaxial nanostructuring [10]).The high performance of PbTe has been linked to 'band convergence' at high temperature [11] and to strong intrinsic phonon anharmonicity [12][13][14], although the extent of the latter is disputed [15,16].However, the environmental toxicity of Pb and the rarity of Te means both materials are unsuitable for mass-produced TEGs, and considerable effort is being devoted to alternatives including alloys of Bi 2 Te 3 and PbTe with other chalcogens [17][18][19][20] and other chalcogenides such as the flagship TE material SnSe [21][22][23][24][25][26].
While chalcogenides have the advantage of favourable electrical properties and low lattice thermal conductivity, a major drawback is that they are prone to oxidation and can degrade or decompose at high temperature [27].On the other hand, oxide materials, while sought after for their chemical stability and large Seebeck coefficients, typically show undesirable high lattice thermal conductivity [28].Oxychalcogenides are an emerging class of TE materials that in principle combine the favorable properties of oxides and chalcogenides.Among the most widely-studied oxychalcogenides are the Bi 2 ChO 2 (Ch = S, Se, Te), which have been examined as potential n-type TEs to complement p-type chalcogenides such as SnSe in TEGs [27].
Given the inherent complexity of balancing the electrical and thermal properties to maximise the ZT, and the relatively small number of studies on Bi 2 SO 2 and Bi 2 TeO 2 to date, more research is required to obtain a complete picture of the thermoelectric properties of the Bi 2 ChO 2 .
We have recently developed a fully ab inito modelling protocol to determine the figure of merit of thermoelectric materials as a function of doping level and temperature by combining accurate models for the electrical-transport properties and the lattice thermal conductivity [37,38].This procedure has been shown to yield accurate results for a growing number of thermoelectric materials including the chalcogenides SnS/SnSe and GeS/GeSe [26,39] and the oxide perovskites CaTiO 3 , SrTiO 3 and BaTiO 3 [40].In this work we apply the same approach to Bi 2 SO 2 , Bi 2 SeO 2 and Bi 2 TeO 2 with both n-and p-type doping to determine the optimal conditions to maximise the ZT.We observe excellent agreement between our calculated properties and recent experimental measurements, and in particular a favourable cancellation of errors between the electrical and thermal transport properties that allows the calculations to yield near-quantitative predictions of the ZT.We find that recent experiments on bulk n-type Bi 2 SeO 2 are likely to have achieved close to optimum performance, whereas the reported ZT of Bi 2 TeO 2 can be improved sixfold by optimising the carrier concentration.Finally, we also predict that heavy p-type doping of Bi 2 SO 2 and Bi 2 SeO 2 could produce much larger ZT ≈ 2.5 competitive with the benchmark thermoelectric SnSe.

Computational methodology
Periodic density-functional theory (DFT) calculations were performed using the Vienna Ab initio Simulation Package (VASP) code [41].
The PBEsol functional [42] was used to describe electron exchange and correlation, as we have previously found that this gives good results for lattice-dynamics calculations [43].While we would normally expect a dispersion correction to be required to accurately describe the structures of layered materials, we found that 'bare' PBEsol could accurately reproduce the crystal structures of the three bismuth oxychalcogenides, which may be due to the interlayer interactions being more electrostatic in nature than in comparable layered chalcogenides.The ion cores were modelled using projector augmented-wave (PAW) pseudopotentials [44,45] with the valence configurations: Bi-6s 2 6p 3 ; O-2s 2 2p 4 ; S-3s 2 3p 4 ; Se-4s 2 4p 4 ; and Te-5s 2 5p 4 .The Kohn-Sham wavefunctions of the valence electrons were represented in a plane-wave basis with a kinetic-energy cutoff of 600 eV, and Brillouin-zone integrations were performed using Γ-centered Monkhorst-Pack k-point meshes [46] with 4 × 4 × 6 subdivisions.Both parameters were chosen to converge the absolute total energy and cell pressure to < 1 meV atom −1 and < 1 kbar (0.1 GPa) respectively.The size of the charge-density grids was determined automatically to avoid aliasing errors (PREC = Accurate in VASP), and non-spherical contributions to the gradient corrections inside the PAW spheres were accounted for (LASPH = .TRUE.).
An initial structure of bulk Bi 2 TeO 2 was obtained from the Inorganic Crystal Structure Database (ICSD: 194 944) [29] and symmetrised using spglib [47].Structures of Bi 2 SO 2 and Bi 2 SeO 2 were derived from this by replacing the chalcogen.Each structure was fully optimised to tolerances of 10 −8 eV and 10 −3 eV Å −1 on the electronic total energy and forces.
Lattice-dynamics and thermal-conductivity calculations were set up and post-processed using the Phonopy and Phono3py software [38,48].The second-and third-order force constants (FC2/FC3) were determined using the supercell finite-displacement method, with 6 × 6 × 2 expansions of the conventional cell (720 atoms) and a displacement step of 10 −2 Å for the FC2 and 3 × 3 × 1 expansions (90 atoms) and a step size of 3 × 10 −2 Å for the FC3.During the single-point force calculations additional support grids with 8 × the number of points as the standard charge-density grids were employed to reduce numerical noise (ADDGRID = .TRUE.).Phonon density of states (DoS) curves g(ν) were computed by interpolating the phonon frequencies onto regular Γ-centered 24 × 24 × 24 q-point grids and applying a Gaussian smearing with a width σ = 0.064 THz, corresponding to a full-width at half-maximim (FWHM) of 5 cm −1 .Phonon dispersion curves ν(q) were computed by interpolating frequencies along strings of q-points passing between the high-symmetry wavevectors in the I4/mmm Brillouin zone.The lattice thermal conductivities were computed using the single-mode relaxation-time approximation (SM-RTA) [38] from modal properties evaluated on regular Γ-centered 16 × 16 × 16 q-point meshes.This mesh was chosen based on explicit testing to converge the κ latt to within ∼3% of the values obtained with larger 24 × 24 × 24 sampling meshes.For the phonon dispersion and κ latt calculations, a conversion to the primitive cell was performed using the transformation matrix: The impact of non-analytical corrections to the phonon frequencies at q → Γ on the κ latt was assessed using the high-frequency dielectric constant ε ∞ and Born effective-charge tensors Z * calculated using the density-functional perturbation theory (DFPT) routines in VASP [49].Electronic-transport calculations were performed using the AMSET code [37].Uniform band-structure calculations were performed using the HSE06 hybrid functional [50] and expanded 8 × 8 × 12 k-point meshes to obtain accurate bandgaps E g and sets of Kohn-Sham wavefunction coefficients.(We note that PBEsol predicts a metallic electronic structure for Bi 2 TeO 2 .)These meshes were increased by a factor of 20 using interpolation for the transport calculations.Electron relaxation times were calculated by summing scattering rates from three different processes, viz.acoustic deformation potential (ADP), polar optical phonon (POP) and ionised impurity (IMP) scattering.Deformation potentials were computed by performing a series of single-point electronic-structure calculations on deformed structures generated using AMSET.For these calculations, the HSE06 functional was used with the 'base' k-point meshes used for the geometry optimisations.Elastic constants were computed using PBEsol and the finite-differences routines in VASP.High-frequency and static dielectric constants ε ∞ and ε s = ε ∞ + ε ionic were determined using the linear-optics and finite-difference routines in VASP, respectively [49].The ε ∞ was determined using the 'base' 4 × 4 × 6 k-point mesh with the HSE06 functional, while the ionic contribution was determined using an expanded 8 × 8 × 12 mesh with PBEsol.The POP frequencies ω po were determined from the phonon frequencies at q = Γ, obtained using Phonopy, and the infrared (IR) activities determined using the Phonopy-Spectoscopy package [51] and Born effective-charge tensors Z * from density-functional perturbation theory (DFPT) calculations with PBEsol [49].(We note that the fourth scattering mechanism modelled by AMSET, piezoelectric (PIE) scattering, is not relevant to Bi 2 ChO 2 since the piezoelectric moduli vanish by symmetry in the centrosymmetric I4/mmm spacegroup.) We also performed additional calculations to determine the formation energies E F of the Bi 2 ChO 2 and to compare them to the related Bi 2 O 3 and Bi 2 Ch 3 .Details of these calculations are provided in section S1 of the Supporting Information.

Structure and lattice dynamics
The optimised lattice parameters of the three Bi 2 ChO 2 are collected in table 1.The calculations on all three structures agree very well with experiments [29,52,53], with a slight 0.5%-2.3%underestimation of the lattice parameters that can be accounted for by the 'athermal' DFT calculations being performed at 0 K without zero-point corrections.(We note that the Bi 2 SO 2 structure reported in [52] is in the Pnma spacegroup, but has a very similar structure to Bi 2 SeO 2 and Bi 2 TeO 2 , and the tetragonal I4/mmm spacegroup can be recovered with a loose symmetry tolerance of 1 Å.)The predicted lattice constants and cell volumes increase in the order Bi 2 SO 2 < Bi 2 SeO 2 < Bi 2 TeO 2 , as expected given the increase in atomic radii from 100 to 140 pm [54], and this trend is also evident in the experimental structures.
The calculated phonon dispersion and atom-projected density of states (PDoS) curves for the three oxy-chalgoenides are shown in figure 2. All three dispersions are free from imaginary harmonic modes, indicating that the structures are dynamically stable.We note that this indicates the I4/mmm structure of Bi 2 SO 2 represents a minimum on the structural potential-energy surface.All three phonon dispersion curves are split into two regions, viz. a low-frequency region composed of modes dominated by motion of the heavier Bi and Ch atoms, and a higher-frequency region dominated by motion of the O atoms.The two regions are separated by a 'phonon bandgap' of ∼1-2.5 THz.As the mass of the chalcogen increases, the Ch peak in the PDoS occupies a narrower frequency range, which has the effect of flattening the associated low-frequency bands in the dispersion.Since the Bi PDoS extends over the whole of the lower-frequency region, this has a smaller effect on the overall bandwidth, although the bandwidth does fall from Bi 2 SO 2 to Bi 2 TeO 2 .More interestingly, whereas the dispersion of the O-based higher-frequency bands is similar in Bi 2 SO 2 and Bi 2 SeO 2 , the dispersion in Bi 2 TeO 2 has a markedly smaller bandwidth, which implies that the chalcogen layer exerts a significant influence on the oxide layer.

Electronic structure and transport properties
The calculated electronic band structure and DoS of Bi 2 TeO 2 obtained using the HSE06 hybrid functional is shown in figure 3, and corresponding plots for Bi 2 SO 2 and Bi 2 SeO 2 are provided in figures S8/S10 in section S3.1 of the SI.The electronic-structure calculations predict all three oxychalgoenides to be indirect-gap semiconductors with bandgaps E g of 1.46, 1.1 and 0.33 eV for Bi 2 SO 2 , Bi 2 SeO 2 and Bi 2 TeO 2 respectively.In Bi 2 SO 2 and Bi 2 SeO 2 the valence-band maximum and conduction-band minimum (VBM/CBM) lie at the k = N and Γ wavevectors, respectively, with a secondary VBM along the M → S path.Increasing the chalcogen mass results in a 'convergence' of a secondary CBM at k = X, such that the VBM and CBM of Bi 2 TeO 2 are located along k = M → S and at k = X.
To investigate the composition of the band edges, we also calculated atom-projected (partial) electronic density of states (DoS) curves using PBEsol (figures S9, S11 and S12).While GGA functionals such as PBEsol typically underestimate the size of the band gap, as in the present case, this is usually due to a rigid shift of the band levels [55] and GGA calculations are therefore typically sufficient to determine the band composition.These calculations show that the valance band edges are composed predominately of O 2p states, with small Table 1.Optimised lattice parameters and crystallographic spacegroups of the three Bi2ChO2 (Ch = S, Se and Te) structures compared to experimental data [29,52,53].The percentage differences to experiments in the optimised unit-cell volumes are shown in parentheses.a The structure in [52] is reported in the Pnma spacegroup, and the independent a and c lattice constants are given for comparison.Our predicted band gaps agree well with the reported E g = 1.5 eV for Bi 2 SO 2 [56] and 0.23 eV for Bi 2 TeO 2 [29].Our predicted bandgap for Bi 2 SeO 2 is however significantly lower than the E g = 1.77 eV obtained from experiments [57].We suspect this may be an issue with the data fitting in the optical measurements, since the value drops significantly to ≈0.9 eV with a very small 2 at.%Te doping, which is much closer to our predicted value.In keeping with this, our calculated value is also in good agreement with the E g = 1.3 eV obtained in a previous theoretical study [58].

System
Having confirmed the accuracy of our calculated bandgaps, we next proceeded to calculate the electrical transport properties using semi-classical Boltzmann transport theory and approximate models for the electron relaxation times [37].
The transport tensors S, σ and κ el are determined for a series of specified extrinsic carrier concentrations ('doping levels') n and temperatures.The intrinsic charge-carrier concentration of Bi 2 SeO 2 is reported to be n ≈ 10 15 cm −3 at room temperature [29], whereas Te-doped Bi 2 SeO 2 and pristine Bi 2 TeO 2 attain values on the order of 10 18 cm −3 in experiments [29,57], and even higher n ≈ 10 19 and 10 20 cm −3 have been reported for Bi 1.9 Ta 0.1 SeO 2 and Bi 2 Se 1−x O 2 Cl x respectively [32,33].On the other hand, n as high as 10 21 cm −3 have been reported for p-type LaCuOSe and BiCuOSe in experiments [59,60].We therefore consider a range of n = 10 16 − 10 21 cm −3 , given that studies with the lower carrier concentrations report poor electrical properties [57].While most experiments focus on n-type Bi 2 ChO 2 , since other oxychalcogenides are often p-type doped we consider both types of doping here [27].We consider a temperature range of 300-900 K, up to the reported decomposition temperature of Bi 2 TeO 2 [29,61].
The diagonal xx = yy and zz components correspond to the transport along the independent a = b and c directions in the tetragonal crystal structure.Since most experiments are performed on films or pressed pellets composed of randomly-oriented crystal grains, as opposed to single crystals, we focus mainly on the scalar averages (2xx + zz)/3 = (xx + yy + zz)/3 and provide a brief discussion of the anisotropy.
Figure 4 shows the averaged electrical properties σ, S, S 2 σ (power factor, PF) and κ el as a function of carrier concentration at a fixed T = 800 K, typical of the higher temperatures probed in measurements [32,33,57], and as a function of temperature for a fixed n = 10 20 cm −3 corresponding to the largest n reported in experiments [32].
At 800 K, the σ of n-doped Bi 2 SO 2 and Bi 2 SeO 2 are predicted to be comparable and very low for n < 10 19 cm −3 , setting a practical minimum doping level.The narrow bandgap of Bi 2 TeO 2 results in a small intrinsic n-type conductivity on the order of 100 S cm −1 below n ≈ 10 19 cm −3 , which is comparable to the measured σ = 75 S cm −1 at n = 1.06 × 10 18 cm −3 and T = 700 K [29].At larger n, the conductivity is predicted to increase sharply, as for the other oxychalcogenides, and Bi 2 TeO 2 is predicted to support the largest σ of the three systems.
With n-type doping, the absolute Seebeck coefficients of Bi  The steep increase in σ at larger n, together with the shallow variation in S, means that the PFs of Bi 2 SO 2 and Bi 2 SeO 2 are dominated by the electrical conductivity and reach a plateau around the largest n = 10 21 examined in the calculations.For Bi 2 TeO 2 , on the other hand, the larger |S| at larger n means the PF shows a broad peak at a lower n ≈ 10 20 cm −3 .We predict Bi 2 SeO 2 and Bi 2 TeO 2 can achieve comparable PFs, with a lower n required for the telluride, whereas we predict Bi 2 SO 2 to attain a maximum n-type PF around of around 60% of that of the heavier oxychalcogenides.
The σ obtained with p-type doping is predicted to be some 2-8 × smaller than with n-type doping.However, for Bi 2 SO 2 and Bi 2 SeO 2 the calculations predict much larger S = 250 − 600 µV K −1 with n = 10 19 − 10 21 cm −3 and ∼5× larger PFs.The calculated S of p-doped Bi 2 SeO 2 changes sign around 7.5 × 10 16 cm −3 , indicating intrinsic n-type behaviour that requires a small extrinsic hole concentration to displace.This results in a pronounced dip in the |S| in figure 4(c).The sharp drop in the |S| of p-doped Bi 2 SO 2 at lower doping levels suggests a similar phenomenon may occur for the sulphide at a slightly lower n than the minimum 10 16 cm −3 we tested.For Bi 2 TeO 2 , a doping level of approx.10 19 cm −3 is needed to offset the intrinsic electron carriers, again evident from a dip in the |S|.Above this n, an enhanced S compared to n-type doping results in a 1.5 × larger peak S 2 σ ≈ 3 mW m −1 K −2 compared to the ∼2 mW m −1 K −2 obtained with n-type doping, albeit at a larger n ≈ 5 × 10 20 cm −3 .The predicted S of Bi 2 SO 2 and Bi 2 SeO 2 at the largest n = 10 21 cm −3 are around ∼50% higher than the peak S of Bi 2 TeO 2 , which offsets the lower σ and leads to considerably larger PFs.
With a fixed n = 10 20 cm −3 the σ of all three systems, with both and n-and p-type doping, shows a metallic-like fall with temperature characteristic of a degenerate semiconductor.For n-and p-doped Bi 2 SO 2 and Bi 2 SeO 2 and n-doped Bi 2 TeO 2 the |S| increase with temperature, whereas the |S| of p-doped Bi 2 TeO 2 shows a shallow maximum around 600 K. Since the n-doped systems have relatively large σ and small |S|, the S 2 term dominates the temperature dependence of the PF and the largest S 2 σ are obtained at higher T. On the other hand, the smaller σ and smaller relative change in S with temperature in the p-doped systems means that the PFs peak at low-to-mid temperatures of < 300 K for Bi 2 SO 2 and approx.350 and 450 K for Bi 2 SeO 2 and Bi 2 TeO 2 respectively.
As a result of the degenerate behaviour at large n, the κ el follows the Wiedemann-Franz law: where L is the Lorentz number.As a result, the κ el rises with n, and the fall in σ with temperature and compensating T term in equation ( 3) leads to a relatively flat variation with temperature.This results in relatively large κ el for p-doped Bi 2 TeO 2 and all three n-doped systems, which will offset the higher PFs in the thermoelectric figure of merit (cf equation ( 1)).
To determine the dominant electronic scattering processes, we analysed the contribution to the overall scattering rates as a function of energy for the three processes considered in the calculations, viz.acoustic deformation potential (ADP), polar optic phonon (POP) and ionised impurity (IMP) scattering (figures S13-S18 in section S3.2 of the SI).In most cases, we find that POP scattering appears to make the largest contribution to the rates, with the exception of p-doped Bi 2 SO 2 for which the rates of ADP and POP scattering are comparable.This is similar to the orthorhombic Pnma phases of SnS and SnSe [26].
We also examined the anisotropy in the electrical transport properties (figures S19-S21).In the tetragonal I4/mmm spacegroup, a = b ̸ = c by symmetry, and transport along the a = b and c directions corresponds to transport within and across the layers respectively (cf figure 1).In all three chalcogenides, the σ is largest along the a = b direction and smallest along the c direction for both n-and p-type doping, with the largest anisotropy in the p-doped systems.For n-type doping the |S| are predicted to be relatively isotropic, whereas for p-type doping the calculations predict significant anisotropy at n comparable to the intrinsic n-type carrier concentrations but isotropic S at larger n.For p-doped Bi 2 TeO 2 this 'crossover' leads to a significant drop in the averaged |S| around n = 10 19 cm −3 .Aside from this the anisotropy in the PFs is dominated by the σ and the S 2 σ are 1.5 − 4× larger along the a = b direction with n-type doping and at least an order of magnitude larger with p-type doping.
Finally, to establish the accuracy of our predicted transport properties we performed a detailed comparison to selected measurements from experimental studies on Bi 2 SeO 2 [30][31][32][33][34] and Bi 2 TeO 2 [29] (section S2.1 of the SI).The calculations predict the measured transport properties of single-crystal Bi 2 SeO 2 [30] to within a factor of two.However, comparison to measurements on more typical consolidated powder samples of Bi 2 SeO 2 and Bi 2 TeO 2 show larger variability, particularly in the predicted power factors, which we ascribed to 'energy filtering' and 'carrier pocketing' effects [62].We therefore suggest that these calculations may serve as a useful reference point for assessing the quality of consolidated powder and/or thin-film samples, and for determining the impact of sample morphology on the electrical transport.

Lattice thermal conductivity
We now consider the lattice (phonon) contribution to the thermal conductivity, κ latt .As for the electrical properties, we focus mainly on the average (scalar) values representative of films or pressed pellets, which we denote κ latt .The calculated κ latt of the three oxychalcogenides is shown in figure 5(a), and table 2 compares values at T = 800 K to similar calculations on the bismuth chalcogenides Bi 2 Ch 3 , [20] SnS and SnSe, [63] and the oxide perovskite SrTiO 3 [40].(This temperature was chosen as a common comparison point based on the fact that most of the materials, including the Bi 2 ChO 2 perform best at high temperature, but SnS and SnSe undergo thermal phase transitions around 800 K [64] and Bi 2 Te 3 melts at around 900 K).
Our calculations predict a κ latt ordering of Bi 2 SO 2 ≈ Bi 2 SeO 2 > Bi 2 TeO 2 .The predicted κ latt at T = 800 K are larger than the equivalent Bi 2 Ch 3 , and the κ latt of Bi 2 SeO 2 is higher than that of both the lower-symmetry Pnma and higher-symmetry R3m phases of Bi 2 Se 3 .The κ latt of Bi 2 SO 2 and Bi 2 SeO 2 are also larger than those of SnS and SnSe.Partial substitution of the chalcogen with oxygen therefore results in a larger κ latt compared to similar chalcogenides.On the other hand, all three Bi 2 ChO 2 are predicted to have significantly lower κ latt Table 2. Comparison of the averaged lattice thermal conductivity κ latt of the three oxychalcogenides examined in this work at T = 800 K to the bismuth chalcogenides Bi2S3, Bi2Se3 and Bi2Te3, [20] SnS and SnSe, [63] and SrTiO3 [40].We also show the decomposition of the κ latt into harmonic and lifetime components κ/τ CRTA and τ CRTA as defined in equation ( 4), together with the average mode Grüneisen parameters (γ).[20] 0.18 0.20 0.91 -Bi2Se3 (R 3m) [20] 0.68 0.30 2.28 -Bi2Te3 (R 3m) [20] 0.33 0.20 than the oxide perovskite SrTiO 3 , indicating the oxychalcogenides partially retain the favourable low lattice thermal conductivity associated with the chalcogenides.This may be due in part to the layered structure.As the Bi 2 Ch 3 are polar solids, we tested the effect of applying non-analytical correction (NAC) the the phonon frequencies at q → Γ on the calculated κ latt and found this to be negligible (figures S22, S24 and S26 in section S3.3 of the SI).

System
We also assessed the cumulative contribution of the phonon modes as a function of frequency to the thermal conductivity at T = 800 K (figures S23, S25 and S26).This analysis shows that 73%-82% of the κ latt is through modes below the phonon bandgaps at ∼6 THz (cf 2), indicating that the majority of the thermal transport is through the acoustic and low-frequency optic modes.
Within the single-mode relaxation-time approximation (SM-RTA), the κ latt are determined as a sum of microscopic contributions κ λ from individual phonon modes λ [38].To better understand the origin of the differences in the κ latt we use the analysis technique developed in our previous work [63,65,66] in which the κ latt is written as the product of a sum over harmonic modal properties and an averaged phonon lifetime according to: where the C λ , ν λ and τ λ are the modal heat capacities, group velocities and lifetimes, respectively, and the sum is normalised by the unit-cell volume V and the number of wavevectors N included in the summation.We note that the harmonic sum κ/τ CRTA is a tensor, and we consider here the average of the diagonal elements as for the κ latt .
The calculated κ latt are compared to the κ/τ CRTA and τ CRTA in figures 5(b) and (c), respectively, and values at 800 K are included in table 2. The κ/τ CRTA saturate at higher temperatures where the C λ reach the Dulong-Petit limit, and differences in the harmonic function at higher T therefore reflect differences in the phonon group velocities.On the other hand, the τ CRTA decrease with temperature, as the increased thermal population of modes at higher T leads to stronger phonon scattering.
The κ/τ CRTA decrease in the order Bi 2 SO 2 > Bi 2 SeO 2 > Bi 2 TeO 2 , which reflects the impact of differences in the chalcogen mass and chemical bond strength on the phonon dispersions (cf figure 2).On the other hand, Bi 2 SeO 2 is predicted to have longer τ CRTA than Bi 2 SO 2 , which balances the lower group velocities to yield similar κ latt .This phenomenon has also been observed in calculations on the Sn chalcogenides (cf table 2), for which the balance is such that a different ordering of the κ latt can be obtained with different calculation parameters [63].Taking the averaged lifetimes as indicative of the strength of the intrinsic phonon anharmonicity, this implies that the (heavier) selenide is less anharmonic than the (lighter) sulphide.On the other hand, Bi 2 TeO 2 has both lower group velocities and shorter τ CRTA than Bi 2 SO 2 and Bi 2 SeO 2 and hence the lowest κ latt of the three systems.
The κ/τ CRTA of the oxychalcogenides are larger than for the Bi 2 Ch 3 and SnS/SnSe.In addition to the chemical bonding and atomic masses, the harmonic function also depends sensitively on the number of atoms n a in the primitive unit cell and the crystal symmetry (spacegroup) [66].The I4/mmm oxychalcogenide structure has the same n a = 5 as the R 3m structure adopted by Bi 2 Se 3 and Bi 2 Te 3 and is lower symmetry, yet has 4-5 × larger κ/τ CRTA .We therefore attribute the larger group velocities in the oxychalcogenides primarily to the lighter O atoms and stronger chemical bonding.To test the latter, we calculated the formation energies E F of the three Bi 2 ChO 2 and compared them to the comparable Bi 2 O 3 and Bi 2 Ch 3 (table S1 in section S1 of the SI).We obtained values between −3.56 and −4.02 eV per F.U. for the Bi 2 ChO 2 , which lie between our calculated E F = −6.52 and −1.15 to −2.04 eV for Bi 2 O 3 and the Bi 2 Ch 3 , respectively, indicating that the oxychalcogenides likely do have bond strengths intermediate between the oxides and chalcogenides.
Conversely, however, the τ CRTA are shorter than comparable chalcogenides, which may be indicative of stronger anharmonicity.In our recent work on the titanate perovskites we found that the strong bonding led to large group velocities (large κ/τ CRTA ) but also to strong phonon-phonon interactions and hence short lifetimes (short τ CRTA ) [40].The analysis here therefore suggests the the κ latt of the oxychalcogenides may be determined by a similar trade-off.This is partially supported by the calculated mode and averaged Grüneisen parameters γ λ /γ (table 2 and figures S28-S33).The averaged values in table 2 show that Bi 2 SO 2 and Bi 2 SeO 2 have comparable γ to SrTiO 3 and larger γ than SnS and SnSe.(The data from our previous study [20] is incompatible with the newer versions of the software we use, so we were not able to calculate the γ for the Bi 2 Ch 3 ).Interestingly, however, the γ of Bi 2 TeO 2 is a much smaller 0.76, comparable to SnS and SnSe, which suggests weaker intrinsic anharmonicity and that the shorter τ CRTA may arise predominantly from a larger 'phase space' for phonon scattering [67].
As for the electrical properties, the κ latt is anisotopic and larger along the a = b direction than the c direction (figures S22, S24 and S26).For Bi 2 SO 2 and Bi 2 SeO 2 the anisotropy is relatively small, with the in-layer κ latt being ∼50% larger than the cross-layer transport.For Bi 2 TeO 2 , on the other hand, the anisotropy is much larger and the κ latt along the a = b direction is around 4 × larger than along the c direction.With reference to the anisotropy in the electrical properties discussed in the previous section, the a = b direction therefore corresponds to the 'easy axis' for both electrical and thermal transport in all three systems.
For Bi 2 SeO 2 the comparisons suggest the calculations may overestimate the measured κ latt at low temperature but that the agreement between the predictions and measurements improve at higher temperature.The measured room-temperature κ latt ≈ 2.5 W m −1 K −1 of single-crystal Bi 2 SeO 2 along the a = b direction [30] is ∼20% smaller than our predicted 3.1 W m −1 K −1 .The measured room-temperature κ latt of consolidated powders range from 1.2-1.6W m −1 K −1 , [31][32][33][34] which are 40%-60% smaller than our predicted average 2.6 W m −1 K −1 but comparable to the smaller predicted κ zz = 1.5 W m −1 K −1 along the c axis.The agreement improves at higher temperature, such that the measurements of κ latt = 0.7 − 0.8 W m −1 K −1 at ~825 K are around 15%-25% smaller than our average of 0.94 W m −1 K −1 .
The calculations predict the room-temperature κ latt of Bi 2 TeO 2 measured in [29] with near quantitative accuracy, but due to the very shallow temperature dependence in the experiments the calculations underestimate the thermal conductivity at 570 K by around 40%.This is likely due to an incomplete separation of the electronic and lattice components of the thermal conductivity using the Wiedemann-Franz law.

Thermoelectric figure of merit
Combining the calculated electrical conductivity, Seebeck coefficient and lattice thermal conductivity allows us to predict the thermelectric figure of merit ZT as a function of doping level and temperature (figure 6).From these plots, we determine a predicted maximum averaged ZT max for each system, for both n-and p-type doping, and the values and associated physical properties are tabulated in table 3.
For the n-type doping typically targeted in experiments, we predict ZT max = 0.33 and 0.45 for Bi 2 SO 2 and Bi 2 SeO 2 , both at n = 2.5 × 10 19 cm −3 and T = 900 K.The predicted ZT max of Bi 2 SeO 2 is obtained at a smaller n but higher temperature to the best experimental ZT max = 0.38 for Bi 1.90 Ta 0.1 SeO 2 [32].The dependence of the electrical and thermal transport properties on the doping level and temperature is such that even modest ZT values require moderate n > 10 19 cm −3 and high temperatures.
For Bi 2 TeO 2 the calculations predict a larger ZT = 0.21 than the reported value of 0.13 under the experimental conditions of n = 1.1 × 10 18 cm −3 and T = 573 K, which is due to a 25% larger predicted absolute Seebeck coefficient.Moreover, the calculations predict that an even higher ZT = 0.81 should be achievable at this temperature with a modest increase in the doping level to n = 10 19 cm −3 .The largest ZT max = 1.05 is obtained at n = 5 × 10 19 cm −3 and T = 900 K, and we predict a ZT > 1 should be achievable over two narrow windows with n between 10 19 -10 20 cm −3 and T from 700-900 K. Bi 2 TeO 2 is therefore the only one of the three oxychalcogenides the calculations predict can attain a benchmark ZT > 1 with n-type doping.
A comparison of the predicted ZT to the ZT max obtained in selected experimental studies [30][31][32][33][34] (section S2.3 of the SI) suggests that the errors in the electrical and thermal transport balance such that the figures of merit are predicted with near-quantitative accuracy.With this in mind, our calculations suggest the ZT max obtained for Bi 2 SeO 2 in the experiments in [31] and [32] may be close to the optimum that can be obtained.
Under the conditions where the predicted ZT max are obtained, Bi 2 SO 2 and Bi 2 SeO 2 have very similar Seebeck coefficients, and the similar ZT max are a balance of a lower conductivity but also a lower thermal conductivity in the selenide.Bi 2 TeO 2 is predicted to have a similar S to the sulphide and selenide but a much larger σ which, despite its larger κ el , results in an overall higher ZT max .
Interestingly, the larger S obtained with p-type doping results in much higher predicted ZT max for Bi 2 SO 2 and Bi 2 SeO 2 .At a comparable n = 5 × 10 19 cm −3 to the doping levels achieved with n-type doping in experiments, we predict ZT max = 0.72 and 1.12 should be achievable for the sulphide and selenide, respectively, which are 2 − 2.5× larger than the ZT max predicted for n-type doping.With even larger n > 10 20 cm −3 to compensate for the lower σ obtained with p-type doping we predict ZT max = 2.53 and 2.62.A benchmark ZT > 1 requires requires heavy doping, with n on the order of 10 20 cm −3 or above, for both Bi 2 SO 2 and Bi 2 SeO 2 , but, depending on n, can be achieved at moderate T > 500 K.While larger than typical n-type doping levels, the heavy p-type doping n are comparable to those achieved in experiments on p-type LaCuOSe and BiCuOSe [59,60].
For Bi 2 TeO 2 , we predict higher p-type ZT max = 1.36 and 1.51 can be achieved with n = 5 × 10 19 and 10 20 cm −3 , and that ZT > 1 can be obtained with a more modest range of doping levels and a slightly wider temperature range of T > 400 K compared to Bi 2 SO 2 and Bi 2 SeO 2 .As for the sulphide and selenide, part of the larger ZT max is due to a larger S, but the higher doping level also leads to a σ comparable to the n-doped material but with a much lower κ el .Similar calculations on p-type SnS and SnSe predicted averaged ZT max = 1.75 and 2.81 respectively.Compared to n-type Bi 2 SO 2 and Bi 2 SeO 2 , the Sn chalcogenides show larger σ and |S|, resulting in ∼4× larger PFs, and comparable or lower κ.With large n, however, p-doped Bi 2 SO 2 is predicted to have a ∼2× larger σ and comparable S to SnS, although the resulting 2 × larger PF is partially mitigated by the higher κ.The conductivity of p-doped Bi 2 SeO 2 is again 1.6 × larger than predicted for SnSe, with a comparable S, but the 50% larger κ results in a similar predicted ZT max .Compared to Sn chalcogenides, the p-doped Bi 2 ChO 2 can therefore sustain a large S at doping levels that lead to a considerably larger σ and PF, but the ZT max are limited by the proportionally larger κ el (cf equation ( 3)) and larger intrinsic κ latt .
Finally, as for the electrical properties the ZT is anisotropic (section S3.4 of the SI).For n-type Bi 2 SO 2 and Bi 2 SeO 2 , the largest and smallest ZT are predicted to occur respectively along the a = b and c directions, but the ZT along the two independent directions are similar and close to the average due to a balance of a larger σ along the a/b direction and a smaller κ latt along the c direction.For n-type Bi 2 TeO 2 , on the other hand, the combination of a reasonable σ and a much lower κ latt along the c direction, together with a smaller anisotropy in the Seebeck coefficient, results in a much higher ZT > 2.5 at the intrinsic n ≈ 10 18 cm −3 and high T.This suggests that Bi 2 TeO 2 could achieve a still larger ZT max if the growth orientation could be controlled.With p-type doping, on the other hand, in all three systems the much larger electrical Table 3. Predicted maximum averaged thermometric figure of merit ZTmax of n-and p-type doped Bi2SO2, Bi2SeO2 and BiTeO2 with the associated doping levels n and temperatures together with the properties in equation (1), viz. the electrical conductivity σ, Seebeck coefficient S, power factor S 2 σ (PF) and electronic, lattice and total thermal conductivity κ el , κ latt and κtot.For each of the p-type doped materials we also list the maximum ZTmax obtained with the doping level limited to a maximum of n = 5 × 10 19 cm −3 on the same order as typical doping levels achieved in experiments on n-type Bi2SeO2 [30][31][32][33] .
ZTmax n (cm −3 ) T(K) σ (S cm −1 ) conductivity more than offsets the larger κ latt along the a/b direction, resulting in ZT up to 20 × larger than along the c axis.In this case, the tetragonal symmetry means that, despite the large anisotropy, these values remain similar to the averages.(We note that in this case we calculate the average ZT from the averaged S, σ and κ el /κ latt and not by averaging the ZT along the three crystal axes, because this is more representative of how the power factor and ZT are calculated from experimental measurements.)

Conclusions
In summary, we have applied a fully ab initio modelling approach to investigate the electrical and thermal transport properties and thermoelectric figure of merit ZT of the n-and p-doped bismuth oxychalcogenides as a function of doping level and temperature.Electrical transport calculations show that the power factors of Bi 2 SO 2 and Bi 2 SeO 2 are limited by the degredation of the Seebeck coefficient due to the heavy doping required to obtain a reasonable electrical conductivity, whereas the low intrinsic |S| of Bi 2 TeO 2 is offset by a larger σ.On the other hand, the alternative p-type doping is predicted to yield a smaller σ but considerably larger S for Bi 2 SO 2 and Bi 2 SeO 2 , potentially resulting in larger PFs at high n, A detailed microscopic analysis of the lattice (phonon) thermal conductivity suggests the stronger chemical bonding and lighter O atoms in the oxychalcogenides lead to larger group velocities but stronger anharmonic phonon-phonon interactions than in the bismuth chalcogenides Bi 2 Ch 3 and SnS/SnSe.
Comparison of the transport properties to experiments suggest these calculations can accurately reproduce measurements on single crystals and provide a reference point for consolidated powders and thin films, and that a favourable cancellation of errors results in near-quantitative predictions of the figure or merit.With this in mind, we predict that recent experiments on Bi 2 SeO 2 may have achieved close to the optimum bulk ZT possible, whereas a sixfold improvement in the reported ZT max of Bi 2 TeO 2 , and an industrially-viable ZT > 1 at high temperature, should be possible by optimising the carrier concentration.More interestingly, the calculations predict that, if the heavy p-type doping demonstrated for other oxychalcogenides is possible for Bi 2 SO 2 and Bi 2 SeO 2 , much larger ZT max ∼ 2.5 competitive with SnSe may be achievable.
Our results therefore motivate further research on optimising the electrical transport of Bi 2 TeO 2 and exploring the possibility of p-doped Bi 2 SO 2 and Bi 2 SeO 2 .A related avenue may be exploring the impact of defects, such as the possibility of stacking faults from offsets between layers, which could occur with some preparation methods such as deposition in thin-film form.Finally, this study further demonstrates the utility of this ab initio approach as a complement to experiments in designing and optimising high-performance thermoelectric materials for waste-heat recovery.

Figure 2 .
Figure 2. Calculated phonon dispersion and density of states (DoS) curves of Bi2SO2 (a), Bi2SeO2 (b) and Bi2TeO2 (c).On each DoS plot the total DoS is shown in black and projections onto the Bi, chalcogen and O atoms are shown as purple, gold and red shaded areas respectively.

Figure 3 .
Figure 3. Calculated electronic band structure and density of states (DoS) of Bi2TeO2.The valence and conduction bands are shown in blue and yellow, respectively, the valence-band maximum (VBM) is set to E = 0 eV, and the VBM and conduction-band minimum (CBM) are marked by green and red circles.
2 SO 2 and Bi 2 SeO 2 are predicted to decrease monotonically with n between n = 10 17 − 10 20 cm −3 , from 600 to 100 µV K −1 , above which we predict a slower decrease to |S| ≈ 50 µV K −1 .On the other hand, the large intrinsic n of Bi 2 TeO 2 results in |S| < 200 µV K −1 at low n with a peak S = −226 µV K −1 at n ≈ 10 19 cm −3 and a fall to comparable |S| to Bi 2 SO 2 and Bi 2 SeO 2 at n = 10 21 cm −3 .

Figure 5 .
Figure 5. Analysis of the lattice thermal conductivity of Bi2SO2 (blue), Bi2SeO2 (red) and Bi2TeO2 (orange).Plot (a) shows the averaged κ latt as a function of temperature.Plots (b) and (c) show the decomposition into harmonic and lifetime components κ/τ CRTA and τ CRTA , respectively, as defined in equation (4).