Computational prediction of the thermoelectric performance of LaZnOPn (Pn = P, As)

High-efficiency earth-abundant thermoelectric materials present an environmentally friendly route for power generation via waste-heat harvesting. Whilst record thermoelectric efficiencies have been achieved in recent years, the leading materials tend to contain toxic and rare elements, including Pb, Te and Bi, which can be problematic commercially. In this study, the thermoelectric capabilities of two earth-abundant candidate thermoelectric materials LaZnOP and LaZnOAs have been investigated using density functional theory (DFT). The electronic structure and band alignment of both materials have been calculated using hybrid DFT including spin–orbit coupling effects, and indicate that both materials should be p-type dopable with highly mobile charge carriers. We analyse the vibrational properties and calculate the lattice thermal conductivity of LaZnOPn (Pn 1⁄4 P, As) using lattice dynamics and find both materials to exhibit highly anisotropic thermal transport, driven by differences in phonon lifetimes and group velocity. We calculate the electronic transport properties, including predicted ZTs, for LaZnOP and LaZnOAs and conclude that both LaZnOP and LaZnOAs are promising candidate materials for hightemperature thermoelectric applications. Furthermore, we provide insight into the effects of nanostructuring during processing on the thermoelectric potential of both materials.


Introduction
Given the rapidly increasing global demand for clean and renewable energy, any effective decarbonisation strategy must involve maximising the energy efficiency of existing technologies. Thermal energy is an unavoidable product of myriad processes, with 52% of energy input globally estimated to be lost as waste heat, signicantly increasing global energy consumption. 1,2 As such, waste-heat harvesting presents a promising opportunity to generate clean energy, through effective utilisation of some of the massive amount of thermal energy currently wasted to the environment. 3,4 The thermoelectric effect describes the ability of a material to directly and reversibly generate electricity when subject to a thermal temperature gradient, with no moving parts or waste products. 5 Thermoelectric generators present a promising opportunity to improve the efficiency of an abundance of processes in applications spanning power generation, transportation, domestic goods, wearable devices powered by heat from the human body, as well as solid-state refrigeration technologies via the Peltier effect. [6][7][8][9][10][11] The efficacy of a thermoelectric material is measured using the dimensionless thermoelectric gure of merit ZT, given by where a, s, T, k are the Seebeck coefficient, electrical conductivity, temperature and thermal conductivity, respectively, and with the thermal conductivity containing both electronic and lattice contributions (k ¼ k elec + k latt ). Evidently, for a thermoelectric material to obtain a high ZT it must have a high power factor (a 2 s) and low thermal conductivity, a non-trivial requirement which is complicated through the interdependence of the Seebeck coefficient, electrical conductivity and the electronic contribution to thermal conductivity through the concentration of charge carriers in a system. 3,12 To this end, intelligent design principles have been explored to identify novel materials with inherently low lattice thermal conductivity and high electrical conductivity. [12][13][14][15] The "phonon-glass electron-crystal" (PGEC) concept put forth by Slack 16 focussed thermoelectric research towards materials with complex crystal structures that facilitate the transport of charge carriers (electron-crystal), but are also poor conductors of heat (phononglass). [17][18][19][20] In addition, processing methods including doping, nanostructuring and alloying have been successfully employed to improve the inherent thermal and electrical properties of materials, to increase thermoelectric efficiencies. [21][22][23] The widespread adoption of thermoelectric generators for industrial and commercial heat harvesting has been somewhat hindered by poor conversion efficiencies (h), expressed as where T h and T c are the hot-and cold-end temperatures across a thermocouple consisting of n-and p-type thermoelectric materials connected in parallel, with an average ZT given by ZT avg . To achieve adequate efficiencies, a thermoelectric material should therefore have a ZT avg > 1 and be thermally stable across a wide temperature range. 24 Despite signicant growth in terms of record ZT values achieved, many of the champion materials, including Bi 2 Te 3 , Pb(Se, Te), contain toxic and rare elements, and can be prone to stability issues at high temperatures. [25][26][27] These problems have further limited the reach of thermoelectrics to mostly niche applications such as space exploration and Peltier coolers for electronics. [28][29][30] The growing interest in thermoelectric devices for renewable energy applications has driven an inux of research into competitive thermoelectric materials comprised of environmentally benign and earth-abundant constituent elements. [31][32][33] Oxide thermoelectrics are of interest for high-temperature applications in particular, as they tend to possess superior chemical and thermal stability properties across a wide temperature range and under oxidising conditions, as well as being low cost. [34][35][36][37] The performance of oxides thermoelectrics has historically lagged behind the records set by intermetallic and chalcogenide-based materials, with few oxides possessing the low thermal conductivity and high power factors required for industrial deployment. 38,39 A range of layered cobalt oxide materials, including Na x CoO 2 and Ca 3 CO 4 O 9 , were found to possess reasonable p-type thermoelectric properties, beginning in the late 1990's. 40,41 Fujita et al. reported a ZT > 1 at 800 K for Na x Co 2 O 2Àd single crystals prepared using a ux method. 41 Polycrystalline Tb-doped Ca 3 Co 4 O 9 prepared using a solid-state synthesis method were found to posses a ZT of 0.74 at 800 K. 42 While the performance of n-type thermoelectrics generally has not been as promising as the p-type counterparts, the rare earthdoped perovskite-based SrTiO 3 materials are considered the archetypal n-type oxide thermoelectrics, with a record ZT of 0.41 at 973 K reported for Sr 1À3x/2 La x TiO 3Àd . 43 In this work, we present a theoretical study of the of the layered mixed-anion quaternary oxypnictides LaZnOP and LaZnOAs, which are of interest as potential earth-abundant thermoelectrics. [44][45][46] Multinary semiconductors, including the I 2 -II-IV-VI 4 and I-II 2 -III-VI 4 families of materials have been increasingly of interest in recent years across elds spanning photovoltaics, photocatalysis, nonlinear optics and thermoelectric energy generation. [47][48][49][50][51][52] The promising performance of quaternary thermoelectric materials has been found to be broadly driven by their low thermal conductivity. Kuo et al. recently reported on four I 2 -II-IV-VI 4 quaternary selenide compounds (BaAg 2 SnSe 4 , BaCu 2 GeSe 4 , BaCu 2 SnSe 4 , SrCu 2 -GeSe 4 ), which were found to possess ultra-low thermal conductivities. 53 The I-II 2 -III-VI 4 material CuMn 2 InTe 4 was also found to possess a low intrinsic thermal conductivity and strong doping potential, indicating promise for thermoelectric applications. 54 The layered quaternary oxyselenide BiCuOSe, a material which is isostructural with LaZnOPn (Pn ¼ P, As), has been reported to have a p-type ZT > 1 when doped with Pb, originating from a very low intrinsic lattice thermal conductivity (0.4-0.5 W m À1 K À1 ), which has been attributed to weak intralayer bonding and resulting strong bond anharmonicity. 55,56 Herein, we have calculated the electronic band structure, band alignment and electronic transport properties including the Seebeck coefficient, electrical conductivity, electronic component of thermal conductivity and power factor across a range of temperatures and doping concentrations from rst principles using hybrid-density functional theory (DFT) including spin-orbit coupling (SOC) effects. In addition, the lattice dynamics, including the phonon dispersion, Grüneisen parameter and lattice thermal conductivity of LaZnOPn (Pn ¼ P, As) were calculated using DFT within the Generalised Gradient Approximation (GGA). An estimate for ZT is predicted from the calculated lattice thermal conductivity and electronic band structure using Boltzmann transport theory within the constant relaxation time approximation (RTA). Our study indicates that LaZnOP and LaZnOAs have the potential to be promising p-type earth-abundant high-temperature thermoelectrics.

Computational methodology
All density functional theory (DFT) calculations were performed within the Vienna Ab initio Simulation Package, 57-60 using projector augmented-wave (PAW) pseudopotential method to account for interactions between core and valence electrons. 61,62 The kinetic energy cutoff for the plane-wave basis set were checked for convergence, with a 450 eV energy cutoff and 6 Â 6 Â 3 G-centred k-point grid deemed to be sufficient to converge the total energy to within 0.01 eV per atom. A more dense 13 Â 13 Â 6 G-centred k-point grid was used for all density of states (DOS) calculations.
The revised Perdew-Burke-Ernzerhof Generalised Gradient Approximation (GGA) functional (PBEsol) 63 was used for all lattice dynamics calculations, including for the second-and third-order interatomic force constants used to calculate the phonon dispersion, Grüneisen parameter and lattice thermal conductivity. PBEsol has been shown to accurately reproduce lattice parameters in solid systems, and has demonstrated good performance in lattice dynamics calculations, while maintaining a relatively low computational cost. [64][65][66][67][68] The Heyd-Scuseria-Ernzerhof hybrid exchange-correlation density functional (HSE06), 69 which contains 25% screened Hartree-Fock exchange at short ranges, was used for all electronic structure, band alignment and electronic transport property calculations. HSE06 has proven success in calculating the structural properties and electronic structure of materials, including oxide semiconductors compared to experiment, eliminating the wellknown band-gap underestimation error association with localdensity approximation (LDA) and GGA functionals. 67,[70][71][72][73] The structure of LaZnOPn (Pn ¼ P, As) reported by Lincke et al. 74,75 was rst optimised using the above functionals to a force conversion criterion of 0.01 eVÅ À1 , with the basis-set expanded by 1.3 times to 585 eV to eliminate Pulay stress. 76 The equilibrium structures of the compounds in this study are provided in an online repository: https://doi.org/10.5281/ zenodo.3715440.
Spin-orbit coupling (SOC) effects were included for the electronic band structure and electronic transport properties, to account for the relativistic effects originating from heavier elements in the materials. Electronic band structures, density of states and phonon dispersions were generated using the opensource Python package: sumo developed by Ganose et al. 77 The electron (e À ) and hole (h + ) effective masses (m*) at a given k were calculated from parabolic tting of the band edges of the conduction and valence bands with eigenvalue E(k), respectively, according to The temperature dependent polaron charge carrier mobility (m) values were calculated using the PolaronMobility.jl package, 78 using the method described by Davies et al. 79 The static dielectric constant was calculated using density functional perturbation theory (DFPT) using PBEsol. The highfrequency dielectric constant was obtained from optical properties calculated with the HSE06 functional, using the method of Furthmüller et al. 80 Alignment of the core state electron eigenvalues to the vacuum level were done using the core-level alignment approach. A slab-gap model consisting of a non-polar 7 layer thick slab was constructed from the HSE06 bulk optimised unit cell with a (010) surface termination and 30Å vacuum separation between slabs, and calculated using HSE06 + SOC. 81 The planar average of the electrostatic potential was calculated using the MacroDensity package and plotted along the (001) direction, with the value of the electrostatic potential at the plateau taken to be the energy of the vacuum. [82][83][84] Vibrational properties, including phonon dispersion and lattice thermal conductivities were calculated using the supercell and nite displacement approaches, as implemented in Phonopy 85 and Phono3py 86 soware packages. The structures of LaZnOPn (Pn ¼ P, As) were relaxed initially using PBEsol, to a force convergence criterion of 0.0001 eVÅ À1 . The harmonic force constants were calculated from displaced 5 Â 5 Â 4 supercells (800 atoms) using a 0.01Å displacement amplitude and with proportional reduction in Brillouin-zone sampling. The Grüneisen parameter was calculated from displaced supercells at the optimised cell volume, and AE1% volume expansions from this volume. Phonon DOS were calculated using Phonopy on a 36 Â 36 Â 17 G-centred q-mesh. The displaced structures for the third-order interatomic force constants (5168 displaced structures) were generated using Phono3py with a reduced supercell size of 4 Â 4 Â 2 (256 atoms) using a 0.03Å displacement amplitude to reduce computational cost. Lattice thermal conductivity values were calculated by solving the linear Boltzmann Transport Equation (BTE) within the single-mode relaxation time approximation (RTA), as implemented in Phono3py, and were deemed to be converged with a 36 Â 36 Â 17 G-centred q-point mesh, including the nonanalytical term correction and isotopic effects. 86 The lifetimes (s) of a phonon at a frequency u and q-point (q ¼ 2pk) were calculated from the imaginary component of self-energy (L), where The group velocities (v g ) and phonon mean free paths (L) can be calculated from the phonon lifetimes, such that The electronic transport properties of LaZnOPn (Pn ¼ P, As), including the Seebeck coefficient, electrical conductivity and electronic contribution to lattice conductivity were calculated using the BoltzTraP code. 87 The wavefunction of the HSE06 optimised structure was recalculated using a much denser kmesh (22 Â 22 Â 10 G-centred) with SOC included, and the electronic transport properties were calculated by solving the linear Boltzmann Transport Equation (BTA) under the constant RTA and using the rigid band approach. The temperaturedependent electronic contribution to the thermal conductivity (k elec ) is calculated using the Wiedemann-Franz law: where L is the Lorentz number and s is the electrical conductivity.
All crystal structures diagrams were produced using the VESTA soware package. 88 3 Results and discussion

Structure and geometry
Both LaZnOP and LaZnOAs crystallise in the tetragonal P4/nmm space group, and are isostructural with the promising thermoelectric material BiCuOSe. [La 3+ O 2À ] + and [Zn 2+ Pn 3À ] À layers are stacked alternately along the (010) direction, and with La-P coordination between layers (Fig. 1). 46 The [Zn 2+ Pn 3À ] + layer consists of ZnP 4 edge-and corner-sharing tetrahedra, while the [La 3+ O 2À ] + layer contains La 3+ coordinated to four equivalent O 2À . Table 1 shows a comparison of the experimental and calculated lattice parameters of the PBEsol and HSE06 relaxed structures, and corresponding bond distances. 44,46 The optimised lattice parameters are in excellent agreement with experiment for both functionals, with a $1% difference for all values except for those calculated for LaZnOP using PBEsol, which shows an underestimation in the c (010) direction of $3%. PBEsol is known to slightly underestimate the lattice parameters in many systems, which is consistent with the results seen here. 89

Electronic structure
The electronic structure of LaZnOP and LaZnOAs were calculated using HSE06 + SOC from the HSE06 bulk optimised structure. The electronic band structures of LaZnOPn (Pn ¼ P, As) are shown in Fig. 2. The calculated and experimental band gaps and electron and hole effective masses, and corresponding temperature dependent polaron mobilities are summarised in Table 2. Both LaZnOP and LaZnOAs have direct band gaps located at the G point, which were found to be 1.435 eV and 1.242 eV, respectively. The slightly smaller band gap of LaZnOAs can be attributed to the decreased electronegativity of As compared to P in LaZnOP, which results in reduced localisation of electrons in the system.
The calculated band gap of LaZnOP is in excellent agreement with experimental results reported by Kayanuma et al. showing a 1.4 eV optical band gap for LaZnOP single crystals, obtained from diffuse reectance spectra data. 44 A band gap of 0.75 eV was calculated from DFT using the GGA functional PBE96 in the same study, indicating issues with the use of GGA functionals for an accurate description of the electronic structure of LaZnOP. Our results are underestimations compared to larger experimental optical band gaps of $1.7 eV and $1.5 eV taken from diffuse reectance spectra reported by the same group for powdered samples of LaZnOP and LaZnOAs, respectively. 45 Larger band gaps of 2.6 and 2.3 eV, respectively, were reported for LaZnOP and LaZnOAs single crystals from electronic absorption spectra measured using a microcrystal spectrophotometer. 74,75 Another theoretical study reports direct band gaps of 0.62 eV and 0.53 eV, respectively, for LaZnOP and LaZnOAs, again calculated using a GGA functional. 90 The vast range in experimental and theoretical band gap values reported for these systems illustrates the importance of a study using hybrid DFT, which has proven success in reproducing the band gap of semiconducting materials compared to experiment. 71,72 One of the complications in terms of maximising the ZT of a material is the contradiction between the electronic properties required for a high Seebeck coefficient and high electrical conductivity. Previous studies have shown that a large density of states and accompanying minimal dispersion at band edges is favourable for maximising the Seebeck coefficient. 15,26 Conversely, highly dispersive bands with light electrons and holes are required for good electrical conduction. Degeneracy and anisotropy have been successfully exploited in band structure engineering techniques, with heavy bands providing sufficient density of states for a large Seebeck coefficient and the light bands reducing electrical resistivity. [91][92][93] Fig. 2 The electronic band structure of (a) LaZnOP and (c) LaZnOAs calculated using HSE06 + SOC, plotted along a high symmetry path in the Brillouin zone according to the Bradley and Cracknell notation. 94 The conduction and valence bands are coloured by orange and blue. The electronic density of states of (b) LaZnOP and (d) LaZnOAs, calculated using HSE06 + SOC. The VBM has been set to 0 eV in all cases.
Electron and hole effective masses, calculated from parabolic tting of the valence and conduction band edges results in small electron and hole effective masses along the G-X and G-M directions in both compounds, with larger effective masses present in the G-Z direction ( Table 2). In both LaZ-nOP and LaZnOAs, the electron effective masses are smaller than hole effective masses in all directions. Between the two compounds, all electron and hole effective masses are lower in LaZnOAs, with the exception of the electron effective mass in the G-Z direction. The charge carrier effective masses, which are inversely proportional to carrier mobilities indicate that electrons and holes will be mobile in directions within the layers and less mobile moving perpendicular to the layers. The calculated electron and hole mobility values are listed in Table 2. All electron and hole mobilities are higher in LaZnOAs than in LaZnOP, which is reective of the higher frequency optical dielectric constant (u optic ) of LaZnOAs. Signicantly higher charge carrier mobilities are seen in the G-X and G-M directions in both materials compared to the G-Z direction, in agreement with the trend observed in the calculated effective masses. The highly anisotropic mobilities are a result of more effective orbital overlap within the [La 3+ O 2À ] + and [Zn 2+ Pn 3À ] À layers, which is reduced in the through-layer direction due to the weaker intra-layer interactions. Fig. 2 shows the total and partial density of states (DOS), calculated with HSE06 + SOC. There are clear similarities between the compositional make-up of the valence and conduction bands in the two compounds. In both systems, the top of the valence band is predominantly composed of pnictide p states, with O p, Zn p, La d and s orbitals also contributing to the valence band. The bottom of the conduction band is predominantly composed of Zn s states (Fig. S1 †), with La d and f, Zn d and Pn (Pn ¼ P, As) p orbital contributions deeper into the conduction band.

Band alignment
The band alignment of LaZnOP and LaZnOAs (Fig. 3) calculated relative to the O 1s core level within the bulk using HSE06 + SOC reveals that both materials have high-lying conduction band minima and valence band maxima, relative to the vacuum level. These features are characteristic of native p-type semiconducting materials, such as Cu 2 O and the isostructural LaCuOSe, which is a degenerate p-type semiconductor when doped with Mg. 71,95-97 High-lying valence band maxima, and the resulting low ionisation potentials drives the formation of holes in semiconducting materials, while lower electron affinities limit n-type dopability due to charge compensation from the higher concentration p-type defects. From these results we propose that LaZnOP and LaZnOAs show promise as p-type thermoelectric materials, although careful assessment of the defect chemistry must be performed to identify a suitable dopant.
An experimental band alignment for LaZnOP from UPS data found the conduction band minimum to be located at $3 eV below the vacuum, compared to our calculated value of 2.5 eV. 44 In the same study, bipolarity was indicated for undoped and Zndoped LaZnOP samples, with the undoped samples indicating n-type behaviour from small negative Seebeck coefficient values (À2.7 mV K À1 ), and Cu-doped samples showing p-type behaviour and positive Seebeck coefficients. However the samples did not exhibit the behaviour expected of band semiconductors, notably the Seebeck coefficients of the undoped samples did not vary with carrier concentration.
Moving from P 2À to As 2À , the level of the conduction band minimum is largely unchanged, with the valence band maximum shiing up by 0.25 eV. As we can see from the density of states, the top of the valence band is predominantly composed of Pn p states, and the shi in ionisation potential can be rationalised by considering the reduced electronegativity of the pnictide anion moving down the group, meaning it is easier to remove an electron from the valence band. On the other hand, the conduction band minimum in both cases contains mainly Zn s states, which is largely unaffected by the choice of pnictide. Table 2 The experimental and HSE06 + SOC calculated direct band gaps of LaZnOPn (Pn ¼ P, As), the charge carrier effective masses calculated from parabolic fitting of the band edges at the conduction band minimum (CBM) and valence band maximum (VBM), and polaron mobility (m) at T ¼ 300 K. Charge carrier effective masses and polaron mobility are given as electron/hole  Fig. 4 shows the phonon dispersion curves of LaZnOPn (Pn ¼ P, As) plotted along the high symmetry path in the Brillouin zone, alongside the phonon partial density of states (PDOS), calculated using the nite displacement method with a 5 Â 5 Â 2 supercell with the PBEsol functional. In both LaZnOP and LaZnOAs, the calculated phonon spectra contain 24 dispersion relations with no imaginary phonon frequencies present, indicating both structures are dynamically stable. As in the electronic band structure, there is notable anisotropy in the phononic properties, with signicantly lower dispersion along the G-Z (through-layer) direction compared to G-X (in-layer) in both systems. As group velocity is calculated from the derivative of the phonon dispersion, at phonon bands are indicative of small phonon group velocities, and are thus favourable for low lattice thermal conductivity. The phonon partial density of states (PDOS) shows that in both systems the acoustic branches are dominated by the heavier cation species, La 3+ and Zn 2+ . Changing the Pn species from P to the heavier As results in a attening of the lowfrequency phonon modes, and corresponding reduced group velocity. In both compounds, the highest energy optical branches (above 9 THz) originate from phonon modes involving mainly O states which are highly dispersive. The large mass contrast present in LaZnOAs, between the heavier ions and the lighter O 2À results in a phonon band gap between 6.2-9 THz. Phonon band gaps have been shown to reduce the thermal conductivity in materials such as in Zintl compounds, via reduction of the group velocity of modes near the band edges. 98 The mode Grüneisen parameter g p provides an estimation of the degree of anharmonicity of a vibrational phonon mode in a system by considering the change in properties of a particular phonon mode with respect to thermal expansion. Anharmonic phonon scattering, which can be driven by weak chemical bonding in a lattice, is responsible for the Umklapp scattering processes which cause a nite lattice thermal conductivity in the absence of external scattering mechanisms such as grain boundaries and point defects. In systems in which the phonon scattering is dominated by Umklapp scattering processes, a large Grüneisen parameter is indicative of a low lattice thermal conductivity, related to the signicant anharmonic phonon scattering processes. [99][100][101] Indeed, the Grüneisen parameter has been used to rationalise the ultra-low lattice thermal conductivity resulting from strongly anharmonic bonding in materials such as SnSe, PbSe and BiCuOSe. 27,102,103 Fig . 4 shows the Grüneisen parameters of the acoustic modes in LaZnOP and LaZnOAs plotted along the high symmetry paths in Fig. 4 Phonon dispersion and phonon density of states of (a) LaZnOP and (b) LaZnOAs calculated using the PBEsol functional, plotted along a high symmetry pathway. Mode Grüneisen parameters of (c) LaZnOP and (d) LaZnOAs calculated using the PBEsol functional, plotted along a high symmetry pathway. In the plot, the longitudinal acoustic mode is highlighted in blue, and the two transverse acoustic phonon modes are highlighted in orange and yellow. the Brillouin zone. The Grüneisen parameters are positive for all acoustic modes in both LaZnOP and LaZnOAs, apart from at the X point in LaZnOAs. The average acoustic mode Grüneisen parameters in LaZnOP and LaZnOAs are 1.67 and 1.14, respectively, with maximum Grüneisen parameters of 3.67 and 2.83 occuring at the Z point in both systems. The increased anharmonicity in modes corresponding to motion in the (010) direction can be attributed to weaker intra-layer interactions, indicating the layered structure helps to facilitate changes in pressure by easily compressing or expanding in the (010) direction. The relatively large Grüneisen parameters in the acoustic phonon modes in LaZnOP and LaZ-nOAs are indicative that the phonon scattering in these modes are dominated by Umklapp processes, and thus should contribute signicantly to a reduction in thermal transport. The average Grüneisen parameter of BiCuOSe, a material which has ultra-low lattice thermal conductivity stemming from strong anharmonic phonon scattering, has been calculated to be 2.4-2.6 from rst principles. 103 PbSe and PbTe, both excellent thermoelectric materials with strong bond anharmonicity, have been reported to have average Grüneisen parameters of 1.65 and 1.45, respectively. 27 Fig. 5(a) shows the temperature dependent lattice thermal conductivity (k latt ) for LaZnOP and LaZnOAs calculated within the single-mode RTA, along the in-layer (a, b) and through-layer (c) directions. In both compounds, the lattice thermal conductivity decreases with temperature between 0 K and 1000 K in both the in-layer and through-layer directions, following a near 1/T relationship. There is signicant directional anisotropy in the thermal transport properties of LaZnOP and LaZnOAs, resulting in an average k c /k a,b ratio of 4.826 and 3.756 respectively at 300 K. The lattice thermal conductivity of LaZnOAs is smaller than LaZnOP in both directions across the entire temperature range, with an average reduction in the in-layer and through-layer lattice thermal conductivity of 31.53% and 11.14%, respectively. Table 3 outlines the contribution of the acoustic and optical phonon modes to the total lattice thermal conductivity in LaZ-nOP and LaZnOAs at 300 K. It is found that the contribution of optical phonon modes to the total lattice thermal conductivity is greater in the through-layer direction than in the in-layer direction in both materials, by 16.5-17.5% respectively. Consideration of the lattice thermal conductivities with the lifetime of the phonon modes (Fig. 6) shows the acoustic phonon relaxation times are longer than the optical phonons in both systems, in some cases by an order of magnitude or greater.

Vibrational properties
Between the two materials, the lifetime of the low-frequency acoustic and optical phonon modes are shorter in LaZnOAs than in LaZnOP, which indicates increased three-phonon scattering driven by the low-frequency As phonon modes in LaZ-nOAs. The low-frequency optical phonons, which are found to have a signicant contribution to the in-layer lattice thermal conductivity in LaZnOAs, were also found to have lower average group velocities, as shown in Fig. S4. † In LaZnOP, the phonon modes in the frequency range dominated by P vibrational modes were found to have relatively long phonon lifetimes and large group velocities in the in-layer direction. We therefore attribute the reduced thermal transport in the in-layer direction in LaZnOAs to a combination of the greater contribution of optical phonons to the in-layer lattice thermal conductivity, which have shorter phonon lifetimes and group velocities. Fig. 5(b) shows the cumulative contribution to lattice thermal conductivity against phonon mean free path (MFP) at 300 K. In both LaZnOP and LaZnOAs, the maximum phonon mean free path is longer in the in-layer direction than the through-layer direction, with 50% of in-layer vibrational phonon modes possessing mean free paths of less than 95.05 nm and 51.08 nm, respectively. On the other hand, 50% of the through-layer directional phonons have mean free paths of less than 136.93 nm and 106.93 nm, respectively. This indicates nanostructuring to sub 100 nm dimensions could effectively reduce the lattice thermal conductivity of both LaZnOP and LaZnOAs signicantly, including by greater than 50% in the higher lattice thermal conductivity direction in both materials, and should therefore be considered during processing.

Thermoelectric properties
The dimensionless gure of merit ZT describes the efficiency of a thermoelectric material, with a good thermoelectric material exhibiting a ZT > 1. The predicted p-type ZT values of LaZnOP and LaZnOAs calculated using BoltzTraP as a function of temperature and carrier concentration are shown in Fig. 7, calculated using the BoltzTraP code. 87 Fig. 6 Frequency dependent phonon lifetimes of LaZnOP and LaZ-nOAs with the longitudinal acoustic and two transverse acoustic phonon mode lifetimes highlighted in blue, orange and yellow, respectively, and the optical phonon relaxation times in grey. Table 3 The percentage contribution of the transverse acoustic (TA 1 , TA 2 ), longitudinal acoustic (LA), and the optical (O) phonon modes to the total lattice thermal conductivity of LaZnOPn (Pn ¼ P, As) in the through-layer and in-layer directions at 300 K The maximum ZT values for LaZnOP and LaZnOAs in the through-layer direction were found to be 2.10 and 1.92 at 1000 K, respectively, corresponding to optimal charge carrier concentrations of 1.52 Â 10 19 cm À3 and 4.24 Â 10 19 cm À3 . The maximum ZTs in the in-layer directions were 0.88 (1.87 Â 10 20 cm À3 ) and 0.95 (1.18 Â 10 20 cm À3 ) at 1000 K. The directional anisotropy in the thermoelectric performance is due to the much higher lattice thermal conductivity in the in-layer direction in both materials, stemming from higher phonon group velocities and lower anharmonicity in the vibrational modes. The high ZT values predicted indicates that LaZnOP and LaZ-nOAs are very promising p-type high-temperature thermoelectric materials, assuming they are dopable to the optimal carrier concentrations needed. Single-crystal Ca 3 Co 4 O 9 , a known p-type oxide thermoelectric, has been reported to have a ZT of 0.87 at 973 K. 104 A maximum ZT of 1.4 at 923 K was reported for modulation-doped BiCuOSe, which is isostructural with LaZ-nOP and LaZnOAs. 105 In addition, the long mean free paths of the heat-carrying phonons in LaZnOP and LaZnOAs indicate that the ZT values calculated here have the potential to be improved through nanostructuring techniques, although this may have an opposing detrimental affect on charge carrier mobilities.

Mode
Comparing the ZT results in the c direction ( Fig. 7(a) and (c)), the maximum ZT decreases as we move from P to As in LaZnOPn, despite the lattice thermal conductivity of the two materials being very similar at this temperature (0.58 and 0.55, respectively). The disparity in the ZT values must therefore come from the electronic transport properties of the two materials. Fig. S2 and S3 † show the temperature dependence of p-type electrical conductivity, Seebeck coef-cient, electronic lattice thermal conductivity and power factor for LaZnOP and LaZnOAs, for a range of charge carrier concentrations from 10 18 -10 21 cm À3 . In both systems, increasing the carrier concentration results in a steady increase in the electronic lattice thermal conductivity and electronic conductivity values, and accompanying decrease in the Seebeck coefficient. The total thermal conductivity in a material is determined from a sum of the lattice and electronic components to thermal conductivity. The predicted electronic lattice thermal conductivity values are very similar between LaZnOP and LaZnOAs, with an increase in electronic lattice thermal conductivity observed with increasing temperature. At high levels of doping, such as those that are required to achieve the maximum predicted ZTs, the through-layer lattice and electronic thermal conductivities are similar in magnitude, while the throughlayer thermal conductivity is dominated by the lattice component. This further indicates that experimental investigation of the processing conditions of both LaZnOP and LaZnOAs must be considered, as phonon analysis indicates that nanostructuring could be a promising route to improving the ZT of both materials by reducing the in-layer lattice thermal conductivity. Electrical conductivity (s) follows the same increasing trend as electronic thermal conductivity with carrier concentration, but slightly decreases with temperature. The Seebeck coefficient is found to be larger in LaZnOP than in LaZnOAs, which can be attributed to the directly proportional relationship between the Seebeck coefficient (S) and the charge carrier effective mass at the valence band maximum (Table 2).
Importantly, to realise the maximum predicted ZT values detailed doping studies are necessary to understand the defect chemistry of LaZnOP and LaZnOAs. However, according to the doping limit rules, which suggest that the position of the band edges of a material control their dopability, the band alignment of LaZnOP and LaZnOAs clearly support a p-type dopable material. [106][107][108][109][110][111] It is known that a p-type material will have a high-lying valence band maximum and corresponding low ionisation potential, which favours hole formation. Conversely, n-type materials tend to possess high electron affinities and lower-lying conduction band maxima. 97,112 In the case of these systems, the ionisation potentials and electron affinities for both materials are relatively small, indicating p-type dopability, providing there are no compensating low formation energy donor defects present. A study of the defect properties will indicate the doping potential of LaZnOP and LaZnOAs, and also provide guidance on optimised experimental growth conditions.

Conclusion
We have investigated the structural, electronic, lattice dynamical and thermoelectric properties of the layered oxypnictide materials LaZnOP and LaZnOAs from rst principles. Electronic structure and band gap calculations using hybrid density functional theory (DFT) indicate that both materials have direct band gaps of 1.435 eV and 1.242 eV, respectively, and band alignments that show a strong propensity for p-type dopability. Such behaviour must be assessed by detailed investigation of the defect chemistry of both materials. We have evaluated that both materials are promising earth-abundant thermoelectric candidates, with predicted maximum p-type ZT c and ZT a,b values of 2.10 and 0.88 for LaZnOP, and 1.92 and 0.95 for LaZnOAs, respectively. The vibrational properties of both materials indicate that they contain strongly anharmonic bonding and have a low intrinsic lattice thermal conductivity, which are also strongly anisotropic. The presence of lowfrequency vibrational modes attributed to As in LaZnOAs appears to increase the three phonon scattering processes and reduce phonon dispersion, thus reducing phonon lifetimes and group velocities, ultimately resulting in a lower lattice thermal conductivity, when compared to LaZnOP. In addition, we propose that these materials are excellent candidates for nanostructure engineering, which we propose will likely result in a reduction in the lattice thermal conductivity, from analysis of the mean free paths of the heat carrying phonons.

Conflicts of interest
There are no conicts to declare.