Giant thermoelectric figure of merit in multivalley high-complexity-factor LaSO

The ingredients of $ZT$ are the electronic transport coefficients (electrical and electronic-thermal conductivity, Seebeck thermopower) obtained from the electronic structure, and the lattice thermal conductivity, which we discuss specifically in Sec. II.3.

The electronic transport coefficients are computed from the ab initio density-functional band structure as a function of temperature and doping in the relaxation-time approximation to the linearized Boltzmann transport equation, an approach known as Bloch-Boltzmann theory allen ; bt2 . We use a model of temperature- and energy-dependent relaxation time described at length in previous papers noi ; sno , as implemented in a publicly-available custom code sno ; thermocasu . We concentrate on $n$-type doping, although $p$-doping also gives interesting values.

Ab initio density-functional structure optimization and band-structure calculations are performed within the generalized-gradient approximation (GGA) using the PBE functional pbe (hereafter simply labelled as GGA) and the projector augmented wave method paw using the VASP code vasp , with the La, S, and O_s datasets at the maximum suggested cutoff. LaSO has an orthorhombic structure (Fig. 1) with space group $Cmca$. It is optimized following quantum forces and stress in a conventional cell containing 8 formula units. The material is made of strongly buckled La-O bonded planes intercalated in the $a$ direction by molecular-crystal-like layers of S dimers. The resulting computed lattice constants are $a$=13.323, $b$=5.950, and $c$=5.945 Å. As expected, transport is found to be less efficient in the $a$ direction than in the $b$$c$ plane. The electronic states are calculated on a (12$\times$24$\times$24) k-points grid. The minimum gap in GGA is 1.3 eV, so the $n$-type thermoelectric coefficients are essentially unaffected by valence states (we checked this using a scissor operator which adjusts the gap to the value obtained by hybrid functional calculations, see Sec. II.7). The band structure will be discussed further in Sec. II.5.

We estimate the lattice thermal conductivity $\kappa_{L}$ using the model of Ref. antimonides . In the present Subsection, the equation numbers refer to that paper. We compute the phonon spectrum and Grüneisen parameters from first principles using density-functional perturbation theory with the Quantum Espresso QE code, using ultrasoft pseudopotentials and GGA. The phonon-phonon relaxation time is computed from Eqs.3 and 7, where for the parameter $\gamma$ we use the Grüneisen density of states obtained from the calculated spectrum, instead of the average Grüneisen value (which we can calculate and find to be 0.78). Beside phonon-phonon scattering, Casimir and disorder scattering can be easily introduced via Eqs. 21 and 22 (assuming for example that disorder originates from O vacancies). The sound velocity is approximated à la Debye (constant below the Debye energy $k_{B}$$\Theta$, zero above); its value is the harmonic average $s$=3134 m/s of the computed small-wavevector acoustic-branch velocities. The thermal conductivity is finally obtained from Eq.2, with energy dependence only.

Fig. 2 shows the lattice thermal conductivity for two Debye temperatures, $\Theta_{S}$=215 K as obtained from our data according to the definition of Eq. 10, and $\Theta_{P}$=335 K with the definition by Pässler passler . The choice of $\Theta$ changes considerably the conductivity, because the phonon-phonon scattering time depends exponentially on $\Theta$. In the following we will compare results obtained with $\kappa_{L}$ for the perfect crystal with these two Debye tempeatures (the higher curves for each of the two $\Theta$’s). Fig. 2 also shows $\kappa_{L}$ for a few combinations of size and disorder parameters; weak disorder or micron-sized-crystalline Casimir scattering clearly have comparatively minor effects on the scale set by the choice of $\Theta$, so we will not expound upon them further.

We compute the transport coefficients with a code thermocasu ; sno employing parts of the BoltzTrap2 bt2 transport code as libraries, and including energy- and temperature-dependent relaxation time. The ab initio bands (assumed rigid, i.e., not changing with doping or temperature) are interpolated by a Fourier-Wannier technique bt2 over a k-point grid with 64 times more points than the ab initio one, i.e. approximately equivalent to a (48$\times$96$\times$96) grid.

As discussed for instance in Ref. noi , a constant relaxation time $\tau$=$\tau_{0}$ neglects relevant physics in the electron scattering (thermal distributions, local masses, energies of phonons, and more). So we adopt a temperature- and energy-dependent relaxation time $\tau(T,E,\mu)=1/P_{\rm imp}+1/P_{\rm ac}+1/P_{\rm polar}$ which enters the kernel $\sigma$ of the integral in Eq. 9 of Ref. bt2 (additional details can be found in Ref. sno ). The $P$’s are the scattering rates of charged impurities, acoustic phonons, and polar optical phonons given in Refs. noi ; ridley1 ; ridley2 . Piezoelectric scattering is zero by symmetry ridley2 ; bilbao . We point out that there is no effective mass approximation in the transport coefficient calculations, except those implicit in the relaxation time.

The parameters entering the model for $\tau$ can be computed directly. The conduction-band deformation potential is $D$=2.1 eV, the density is 5270 kg$\cdot$m^-3 and the effective conduction mass is set to a conservatively large $m_{c}^{*}$=0.5 $m_{e}$. The phonon energies, dielectric tensor, and sound velocity are obtained from linear-response calculations QE as in Sec. II.3. The sound velocity is 3134 m/s as discussed in the previous Section. For the dielectric constants we use the harmonic averages of the (very similar) $b$ and $c$ components, $\varepsilon_{\infty}$=6.23 and $\varepsilon_{\rm static}$=18.73 (as it turns out, the $a$ direction is quite irrelevant for thermoelectricity).

A key point is the scattering by polar (i.e. optical) phonons, which dominates the total scattering rate. The LO phonon energies are calculated directly using the linear-response approach (Sec. II.3) with non-analiticity along the crystallographic directions. Since LO phonons scatter electrons with momentum parallel to the LO polarization vector, a direction-dependent polar-phonon $\tau$ would be needed. To simplify the procedure, we use a set of effective LO energies (20, 27, and 30 meV) obtained as averages of the three LO groups of frequencies over the $a$, $b$, and $c$ directions. We deem this simplification to be quite acceptable in view of the other significant uncertainties in the calculations, such as the choice of $\Theta_{D}$ in the lattice thermal conductivity, or our choice of using a single phonon replica in the polar-phonon scattering rate.

The relaxation time $\tau$($E$,$T$) is sketched vs. $E$ in Fig. 3. As typical of polar insulators, the polar-phonon scattering dominates, and its downward jump across the LO phonon energies is in the low energy region relevant for transport. Here, we use this relaxation time accounting self-consistently for chemical potential changes, and do not employ any average-time or constant-time approximations.

We now discuss the band structure of LaSO, and two versions of the complexity factor. The first, $C_{b}$, is a geometrical value derived from the number of valleys and the masses from the band structure. The second, $C_{t}$, is a transport value obtained a posteriori from the calculated transport properties, with the procedure of Ref. gibbs . The data are summarized in Fig. 4, reporting the conduction-band Fermi surface of $n$-doped LaSO at several chemical potentials, and Fig. 5, displaying the bands, the carrier density (the conduction density of states multiplied by the Fermi distribution), and the complexity factor $C_{t}$.

The low-energy valleys of the LaSO conduction band occur at four internal points of the Brillouin zone (BZ) and at four zone-border points (which count as two internal points). By inspection, the first four conduction states (marked in different colors in Fig. 4) provide 16 available valleys within about 150 meV of the band edge. These are visible in the Fermi surface in Fig. 4. Four valleys provided by the first band are internal to the BZ. Eight valleys, from the second band, are located on the square faces and therefore shared with adjacent BZs (hence their wight is 1/2). Eight valleys come from the third band (specifically on the segments Z-Z₁ and S-Y/S-Y₁, Fig. 5), and eight more from the fourth band on the same segment.

The same result is found by counting the relevant bands in Fig. 5. The stationary point on the P–U segment and its three symmetry partners contribute a total of 8 bands; points S and Z with their two symmetry partners contribute 2 bands each (accounting for their being at zone border); that amounts again to 16 available band minima.

The carrier density vs. energy (in the central panel of Fig. 5) shows that all these bands are occupied at the chemical potentials and temperatures of interest (i.e., $\mu$ near the lowest conduction edge, and temperatures of the order of 1 to 3 times room temperature). Accordingly, all the bands just mentioned may be considered active (i.e., contributing to the trasport), so that the effective total multiplicity of the occupied valleys is $N_{v}$=16. From Fig. 5, we can also infer that the optimal doping level, namely the doping at which $ZT$ is maximal for a given $T$, will probably fall in the low 10²⁰ cm^-3, and that the Seebeck coefficient may be interesting due to the fast rise of the density of states near the band edge.

As it turns out, the relevant directions for conduction are in the basal $b$-$c$ plane of the material. The masses in the in-plane directions are quite isotropic, so $K_{v}$$\simeq$1, and the geometrical complexity factor is $C_{b}$$\simeq$16. If the $a$ component of the mass tensor were to be considered, a larger anisotropy would arise, resulting in $K$$\simeq$1.2-1.3, i.e., a $C_{b}$ of about 20. (The anisotropy can be appreciated, for example, from the different curvature of the bands at point S, respectively along the Y-S-Y₁ segment and along the U-S-R segment.)

The transport complexity factor $C_{t}$ is in the rightmost panel of Fig. 5. It is computed from the calculated transport coefficients in the basal plane as outlined in Ref. gibbs using a constant relaxation time $\tau_{0}$=10 fs, purely for consistency with Ref. gibbs . We set the scattering parameter $\lambda$=1 in the Seebeck coefficient model of Ref. gibbs , to take into account the dominant polar scattering we have in the present case. Similarly to most quantities in thermoelectricity, $C_{t}$ is a temperature- and chemical potential-dependent tensor. To compare it with the geometric $C_{b}$ (a scalar) we pick a low T=300 K and $\mu$=0, and average over in-plane directions. $C_{t}$ at 300 K is relatively flat at low energy, and its energy average is roughly 15, in quite decent agreement with $C_{b}$=16 obtained above. This is in line with our having considered just the in-plane, largely isotropic transport. We discuss in Sec. II.6 (especially with reference to Figs. 6 and 7) the connection of our $C$ values with the rule of thumb of Ref. gibbs .

Fig. 6 reports the Seebeck coefficient, the electrical conductivity, the electronic thermal conductivity, and the power factor vs. doping and parametrized by $T$. Fig. 7 displays the same quantities vs. $T$ at optimal doping (i.e. the doping at which $ZT$ is a maximum at the given temperature). For simplicity, only the $b$ component is plotted in Fig. 6. As seen in Fig. 7, the $c$ and $b$ components are very close, and the $a$ component ends up producing small power factor and $ZT$, so it can be ignored.

The large conductivity and Seebeck coefficient result in a large power factor for the in-plane transport, with a maximum of 15 mW/(K²m) at 400 K. We can now make contact between the large complexity factors $C_{t}$$\sim$$C_{b}$$\simeq$16 discussed in the previous Section and the rule of thumb of Gibbs et al. gibbs . For this complexity factor, Fig. 3 of Ref. gibbs suggests that the expected maximum power factor should be between 6 and 22 mW/(K²m). Using the setting of Ref. gibbs ($\tau$=$\tau_{0}$=10 fs and $T$=600 K), we indeed obtain a power factor of 21 mW/(K² m); with the full relaxation time treatment, the power factor at 600 K (Fig. 7) is about 12 mW/(K²m). In both cases, our results are consistent with the general prediction gibbs .

Finally, in Fig. 8, we show the $ZT$ tensor for the two instances of lattice thermal conductivity discussed in Sec. II.3. The left panels report the diagonal components of $ZT$ as a function of $T$ at optimal doping, and the right panels the $b$ component as a function of doping for different temperatures. The $b$ and $c$ components differ by only a fraction of a percent. The main results in this Figure are i) $ZT$ is very large, reaching values well over 3 and, respectively, 6 at high $T$ for the two versions of the lattice thermal conductivity; ii) the optimal doping is in the low to mid $10^{20}$ cm^-3, depending on $\kappa_{L}$ and $T$. $ZT$ may still be rather interesting even at lower doping: for example, it is already above 2 at 800 K and 2$\times$10¹⁹ cm^-3 (Fig. 8, upper panel). We recall again that all coefficients refer to electrons in perfect-crystal bands, subject to scattering from phonons and charged impurities, so that no scattering is accounted for from disorder, dislocations, neutral impurities, extended defects, etc. which could affect transport in ways we cannot quantify.

Given the relatively high carrier density required to obtain interesting $ZT$s, we now look into the possibility of $n$-type doping of LaSO. Screening a number of options, we find that Hf seems to be a reasonable candidate donor.

Dopability is difficult to assess in any material, and LaSO is no exception. In particular, here we do not address possible compensation by native defects or other contaminants, but only the solubility and ionization of donors. Also, our discussion is based on equilibrium thermodynamics, so the possibility remains that epitaxial growth (which often occurs out of equilibrium), ion implantation, and certain kinds of diffusion doping processes may do better than we predict here.

Solubility and carrier concentrations are estimated from the formation energies and thermal ionization levels obtained via ab initio calculation. We use VASP, on 144-atom supercells with a 4$\times$4$\times$4 k-point grid and the GGA, the details being the same as in the calculations in the previous Sections. For charged states we use the simple monopole correction by Leslie and Gillan LG , with the dielectric constant calculated in Sec. II.4.

We concentrate on potential donors substituting on the La site (namely Zr, Hf, Ce, Sb, Bi, Sn, and Si) as we find that potential substitutions for O or S, such as F or Cl, tend to go interstitial or have large formation energy. Si can be discarded offhand because of its huge formation energy. Sb, Bi, and Sn can be neglected as well since they have no state (in particular no donor state) in the vicinity of the gap. This effective trivalent behavior is presumably due to the hybridization cost of the $s$ orbital, and the pyramidal bonding geometry of the La site which does not favor $sp$ hybrids.

For the remaining dopants Zr, Hf, and Ce, we find the thermal levels reported in Table 1. Neutral Ce has an electron in an $f$ orbital when substituting for trivalent La, and its donor level is deep in GGA. This will not improve when using non-local-density methods (such as hybrid functionals) which remove self-interaction and tend to lower the energy of very localized occupied states. Luckily, instead, Zr and Hf have shallow thermal levels, in fact lying just above the conduction edge in GGA (although there as usual are large uncertainties in these estimates, at the very least $\pm$0.1-0.2 eV). A python notebook with the formation energies and their post processing is at https://gitlab.com/vfiore/laso-eform.

To check for the effects of more advanced exchange-correlation functionals, we calculated the thermal level for Zr using the hybrid HSE functional hse , and found a value similar to GGA. This is further confirmed by the electronic band structure, which exhibits impurity-related resonant states within the lower portion of the conduction band, at the same position both in GGA and HSE. In passing, the gaps are 1.3 eV in GGA and 2.9 eV and HSE, the difference being within 15% of that predicted by the dielectric correction rule FB .

Another interesting result of these calculations is that the upper valence band and the bottom conduction band are predominantly sulphur-like, so that conduction effectively occurs in-plane in the S layers. Carriers living in the S layer may thus be partially decoupled from charged impurities in the La-O layers, leading to a reduced impurity scattering.

To present synthetically the results as a function of an operational quantity, we report in Fig. 9 the density of carriers vs. the oxygen partial pressure at a typical device operating temperature $T_{\rm op}$=700 K, and with Hf incorporated in LaSO in thermal equilibrium during growth at temperature $T_{\rm gr}$. The rationale for this is that, as we now discuss in more detail, the partial pressure is related to the formation energy of the donor via the oxygen chemical potential, hence determines the doping level.

The formation energy is related to the oxygen chemical potential, hence to the oxygen partial pressure. The dominant solubility limit for Zr and Hf is the formation of dioxides. Since La forms a sesquioxide, the Hf substitution causes an excess of oxygen. Therefore the formation energy $E_{f}$ increases with the oxygen chemical potential, and it is at its largest in oxygen-rich conditions, namely

with $\mu_{\rm O}$$\leq$$\mu_{\rm O_{2}}$ and

where $E_{\rm def}$ and $E_{\rm bulk}$ are the defected and pristine supercell energies, the $\mu$’s the chemical potentials of the bulk elements, $\mu_{\rm O_{2}}$ the chemical potential of O in a O₂ molecule, and $\Delta H$ the oxide formation enthalpies. $E_{\rm def}$ and $E_{\rm bulk}$ are calculated directly, the bulk chemical potentials are taken from experiment, and formation enthalpies are from the Materials Project mp . Finally, the O partial pressure at growth temperature $T_{\rm gr}$ is related to the oxygen chemical potential $\mu_{\rm O}$ by the standard relation $p$=$p_{0}\,\exp{((\mu_{\rm O}-\mu_{\rm O_{2}})/k_{B}T_{gr})}$ and $p_{0}$=19 kPa the normal-conditions oxygen partial pressure.

The dopant density increases as the chemical potential goes more and more negative, i.e. the partial pressure is reduced, and this offers some leeway to increase the density by adopting O-lean growth conditions. As Fig. 9 shows, Hf can indeed produce useful carrier densities in the 10²⁰ cm^-3 range for low but not unreasonable O partial pressures. Since the O vacancy is most likely a deep donor as it is in most oxides, oxygen deficiency should not cause notable counterdoping. Zr, in turn, is unfortunately ruled out by its comparatively larger formation energy.