Strong anharmonicity and high thermoelectric efficiency in high temperature SnS from first-principles

TL;DR

This study uses first-principles calculations to demonstrate that high-temperature SnS exhibits strong anharmonic vibrational properties similar to SnSe, indicating its potential as an efficient thermoelectric material.

Contribution

The paper reveals that SnS has high thermoelectric efficiency at high temperatures due to its anharmonic vibrational properties, comparable to SnSe, supported by first-principles calculations.

Findings

SnS and SnSe have similar electronic and vibrational properties in the high-temperature phase.

SnS exhibits strong anharmonic effects affecting phonon spectra and thermal conductivity.

Non-perturbative anharmonic effects are crucial for accurate thermal conductivity calculations.

Abstract

SnS and SnSe are isoelectronic materials with a common phase diagram. Recently, SnSe was found to be the most efficient intrinsic thermoelectric material in its high-temperature Cmcm phase above 800 K. Making use of first-principles calculations, here we show that the electronic and vibrational properties of both materials are very similar in this phase and, consequently, SnS is also expected to have a high thermoelectric figure of merit at high temperature in its Cmcm phase. In fact, the electronic power factor and lattice thermal conductivity are comparable for both materials, which ensures a similar figure of merit. As in the case of SnSe, the vibrational properties of SnS in the Cmcm phase are far from trivial and are dominated by huge anharmonic effects. Its phonon spectra are strongly renormalized by anharmonicity and the spectral functions of some particular in-plane modes depict anomalous non-lorentzian profiles. Finally, we show that non-perturbative anharmonic effects in the third-order force-constants are crucial in the calculation of the lattice thermal conductivity. Our results motivate new experiments in the high temperature regime to measure the figure of merit of SnS.

Equations

$$ZT=S^2\sigma T/\kappa ,$$

$$\sigma (T,\mu )$$

$$S(T,\mu )$$

$$\Sigma (\epsilon )=\frac{1}{\Omega N_{\mathbf{k}}}\sum _{n\mathbf{k}}{\tau }_{n\mathbf{k}}^e|{\mathbf{v}}_{n\mathbf{k}}{|}^2\delta (\epsilon -{\epsilon }_{n\mathbf{k}}),$$

$$\mathcal{H}=\sum _a\frac{P_a^2}{2M_a}+\frac{1}{2}\sum _{ab}(\mathbf{R}-\mathcal{R}{)}_a{\Phi }_{ab}(\mathbf{R}-\mathcal{R}{)}_b,$$

$$\frac{{\partial }^{2}F}{\partial \mathcal{R}\partial \mathcal{R}}=\mathbf{\Phi }+\stackrel{(3)}{\mathbf{\Phi }}\mathbf{\Lambda }(0)[1-\stackrel{(4)}{\mathbf{\Phi }}\mathbf{\Lambda }(0){]}^{-1}\stackrel{(3)}{\mathbf{\Phi }},$$

$$\stackrel{(n)}{\mathbf{\Phi }}={⟨\frac{{\partial }^{n}V}{\partial {\mathbf{R}}^n}⟩}_{{\rho }_{\mathcal{H}}}.$$

$$\stackrel{(n)}{\varphi }={\left[\frac{{\partial }^{n}V}{\partial {\mathbf{R}}^n}\right]}_0.$$

$$\begin{gathered} {\Lambda }^{abcd}(z)=-\frac{1}{2}\sum _{\mu \nu }\tilde{F}(z,{\tilde{\Omega }}_{\mu },{\tilde{\Omega }}_{\nu })\times \\ \sqrt{\frac{\hslash }{2M_a{\tilde{\Omega }}_{\mu }}}{e}_{\mu }^a\sqrt{\frac{\hslash }{2M_b{\tilde{\Omega }}_{\nu }}}{e}_{\nu }^b\sqrt{\frac{\hslash }{2M_c{\tilde{\Omega }}_{\mu }}}{e}_{\mu }^c\sqrt{\frac{\hslash }{2M_d{\tilde{\Omega }}_{\nu }}}{e}_{\nu }^d, \end{gathered}$$

$$D_{ab}^{(S)}={\Phi }_{ab}/\sqrt{M_a{M}_b},$$

$$\begin{gathered} \tilde{F}(z,{\tilde{\Omega }}_{\mu },{\tilde{\Omega }}_{\nu })=\frac{2}{\hslash }[\frac{({\tilde{\Omega }}_{\mu }+{\tilde{\Omega }}_{\nu })[1+n_B({\tilde{\Omega }}_{\mu })+n_B({\tilde{\Omega }}_{\nu })]}{({\tilde{\Omega }}_{\mu }+{\tilde{\Omega }}_{\nu }{)}^2-z^2}- \\ \frac{({\tilde{\Omega }}_{\mu }-{\tilde{\Omega }}_{\nu })[n_B({\tilde{\Omega }}_{\mu })-n_B({\tilde{\Omega }}_{\nu })]}{({\tilde{\Omega }}_{\mu }-{\tilde{\Omega }}_{\nu }{)}^2-z^2}], \end{gathered}$$

$${\mathbf{G}}^{-1}(z)=z^{2}1-{\mathbf{M}}^{-\frac{1}{2}}\mathbf{\Phi }{\mathbf{M}}^{-\frac{1}{2}}-\mathbf{\Pi }(z).$$

$$D_{ab}^{(F)}=\frac{1}{\sqrt{M_a{M}_b}}\frac{{\partial }^{2}F}{\partial {\mathcal{R}}_a\partial {\mathcal{R}}_b}.$$

$$\Pi (z)={\mathbf{M}}^{-\frac{1}{2}}\stackrel{(3)}{\mathbf{\Phi }}\mathbf{\Lambda }(z)[1-\stackrel{(4)}{\mathbf{\Phi }}\mathbf{\Lambda }(z){]}^{-1}\stackrel{(3)}{\mathbf{\Phi }}{\mathbf{M}}^{-\frac{1}{2}},$$

$$\mathbf{\Pi }(z)\approx {\mathbf{\Pi }}^{(B)}(z)={\mathbf{M}}^{-\frac{1}{2}}\stackrel{(3)}{\mathbf{\Phi }}\mathbf{\Lambda }(z)\stackrel{(3)}{\mathbf{\Phi }}{\mathbf{M}}^{-\frac{1}{2}}.$$

$${\Pi }_{\mu }(\mathbf{q},\omega )=\sum _{\bar{a}\bar{b}}{e}_{\mu }^{\bar{a}}(-\mathbf{q}){\Pi }_{\bar{a}\bar{b}}(\mathbf{q},\omega +i{0}^{+})e_{\mu }^{\bar{b}}(\mathbf{q}).$$

$$\begin{gathered} \sigma (\mathbf{q},\omega )= \\ \frac{1}{\pi }\sum _{\mu }\frac{-\omega Im{\Pi }_{\mu }(\mathbf{q},\omega )}{({\omega }^2-{\tilde{\Omega }}_{\mu }^2(\mathbf{q})-Re{\Pi }_{\mu }(\mathbf{q},\omega ){)}^2+(Im{\Pi }_{\mu }(\mathbf{q},\omega ){)}^2}. \end{gathered}$$

$${\Omega }_{\mu }^2(\mathbf{q})={\tilde{\Omega }}_{\mu }^2(\mathbf{q})+Re{\Pi }_{\mu }(\mathbf{q},{\tilde{\Omega }}_{\mu }(\mathbf{q})),$$

$$\begin{gathered} {\kappa }_l^{\alpha \beta }=\frac{{\hslash }^2}{\Omega N_{\mathbf{q}}{k}_B{T}^2}\times \\ \sum _{\mathbf{q}\mu }{c}_{\mu }^{\alpha }(\mathbf{q})c_{\mu }^{\beta }(\mathbf{q}){\tilde{\Omega }}_{\mu }^2(\mathbf{q})n_B({\tilde{\Omega }}_{\mu }(\mathbf{q}))\left[n_B({\tilde{\Omega }}_{\mu }(\mathbf{q}))+1\right]{\tau }_{\mu }(\mathbf{q}), \end{gathered}$$

$$\begin{gathered} \frac{1}{{\tau }_{\mu }(\mathbf{q})}=2{\gamma }_{\mu }(\mathbf{q})=\frac{\pi }{{\hslash }^2{N}_{\mathbf{q}}}\sum _{{\mathbf{q}}^{\prime }\nu \eta }|{\stackrel{(3)}{\Phi }}_{\mu \nu \eta }(\mathbf{q},{\mathbf{q}}^{\prime },{\mathbf{q}}^{\prime \prime }){|}^2 \\ \times [(1+n_B({\tilde{\Omega }}_{\nu }({\mathbf{q}}^{\prime }))+n_B({\tilde{\Omega }}_{\eta }({\mathbf{q}}^{\prime \prime })))\delta ({\tilde{\Omega }}_{\mu }(\mathbf{q})-{\tilde{\Omega }}_{\nu }({\mathbf{q}}^{\prime })-{\tilde{\Omega }}_{\eta }({\mathbf{q}}^{\prime \prime })) \\ +2(n_B({\tilde{\Omega }}_{\nu }({\mathbf{q}}^{\prime }))-n_B({\tilde{\Omega }}_{\eta }({\mathbf{q}}^{\prime \prime })))\delta ({\tilde{\Omega }}_{\mu }(\mathbf{q})+{\tilde{\Omega }}_{\nu }({\mathbf{q}}^{\prime })-{\tilde{\Omega }}_{\eta }({\mathbf{q}}^{\prime \prime }))], \end{gathered}$$

$$\begin{gathered} {\stackrel{(3)}{\Phi }}_{\mu \nu \eta }(\mathbf{q},{\mathbf{q}}^{\prime },{\mathbf{q}}^{\prime \prime })=\sum _{\bar{a}\bar{b}\bar{c}}\sqrt{\frac{{\hslash }^3}{8M_{\bar{a}}{M}_{\bar{b}}{M}_{\bar{c}}{\tilde{\Omega }}_{\mu }(\mathbf{q}){\tilde{\Omega }}_{\nu }({\mathbf{q}}^{\prime }){\tilde{\Omega }}_{\eta }({\mathbf{q}}^{\prime \prime })}} \\ \times e_{\mu }^{\bar{a}}(\mathbf{q})e_{\nu }^{\bar{b}}({\mathbf{q}}^{\prime })e_{\eta }^{\bar{c}}({\mathbf{q}}^{\prime \prime }){\stackrel{(3)}{\Phi }}_{\bar{a}\bar{b}\bar{c}}(\mathbf{q},{\mathbf{q}}^{\prime },{\mathbf{q}}^{\prime \prime }), \end{gathered}$$

Full text

††preprint: APS/123-QED

I Introduction

Thermoelectricity is a technologically interesting material property that allows to transform residual heat into useful electricityGoldsmid (2010); Behnia (2015). The efficiency of this energy transformation is controlled by the dimensionless figure of merit

[TABLE]

where $S$ is the Seebeck coefficient, $\sigma$ the electrical conductivity, $T$ the temperature, and $\kappa=\kappa_{e}+\kappa_{l}$ the sum of electronic $\kappa_{e}$ and lattice $\kappa_{l}$ thermal conductivities. Therefore, a good thermoelectric performance requires a high power factor $P_{F}=S^{2}\sigma$ together with a low thermal conductivity.

Monochalcogenides have proven to be efficient thermoelectric materialsHeremans et al. (2008); Zhang et al. (2013); Yang et al. (2008); Cho et al. (2011) mainly due to their strongly anharmonic lattice that implies a low lattice thermal conductivityDelaire et al. (2011); Li et al. (2014); Iizumi et al. (1975); O’Neill et al. (2017); Ribeiro et al. (2018). PbTe is an appropriate example of the potential technological relevance of thermoelectric monochalcogenides: it shows a high $ZT$ in the $600-800$ K temperature rangeRavich (2013), as high as 2.2 when nanostructuredHsu et al. (2004), and has been successfully applied in spacecraftsRowe (2018). In the last years SnSe has attracted a great deal of attention since it was measured to be the most efficient intrinsic thermoelectric materialZhao et al. (2014). Its figure of merit soars to $2.6$ after a structural phase transitionZhao et al. (2014); Adouby (1998); Chattopadhyay et al. (1986); Von Schnering and Wiedemeier (1981) at around $800$ K from the low-symmetry $Pnma$ phase to the high-symmetry $Cmcm$. In the high-symmetry phase the electronic band gap is reduced without affecting its ultralow thermal conductivity, providing the record $ZT$. A recent theoretical work shows that the phase transitionAseginolaza et al. (2019) is second order and non-perturbative anharmonicity is very important to get a thermal conductivity in agreement with experiments.

SnS is isoelectronic to SnSe and shows very similar electronic and vibrational propertiesGuo et al. (2015); Chattopadhyay et al. (1986); Von Schnering and Wiedemeier (1981) at low temperatures. Experimentally it also shows a phase transitionChattopadhyay et al. (1986); Von Schnering and Wiedemeier (1981) from the $Pnma$ to the $Cmcm$ structure and a very low thermal conductivity in the formerHe et al. (2018); Tan et al. (2014) phase. Therefore, it is expected to be a very efficient thermoelectric material in the high temperature phase, which together with the fact that S is a much more earth abundant element than Se, makes it a very interesting candidate for technological applications. Actually, in Refs. He et al., 2018; Tan et al., 2014 it is shown how the $ZT$ of undoped SnS increases very fast before the phase transition as in the case of SnSe. However, as far as we are aware, there are no experimental transport measurements for the high temperature phase of SnS. First-principles calculations of its thermoelectric properties are also absent in the literature, hindered by the unstable modes obtained within the harmonic approximation as in the case of SnSeSkelton et al. (2016); Dewandre et al. (2016).

In this work, by performing ab initio calculations we propose that $Cmcm$ SnS is expected to be a very efficient intrinsic thermoelectric material, as good as SnSe in this phase. We show that the $P_{F}$ of SnSe and SnS are expected to be very similar in this phase, as long as the electronic relaxation time is similar in both materials. By including anharmonicity in the phonon calculation at a non-perturbative level within the Stochastic Self-Consistent Harmonic ApproximationErrea et al. (2014); Bianco et al. (2017); Monacelli et al. (2018) (SSCHA), we show that the phonon spectrum of SnS suffers a strong anharmonic renormalization. The phase transition is driven by the collapse of a zone-border phonon. Anharmonicity is so large that the spectral function of some vibrational modes deviates from the Lorentzian-like shape as it happens in other monochalcogenidesLi et al. (2014); Ribeiro et al. (2018). Finally, we calculate the lattice thermal conductivity of $Cmcm$ SnS obtaining ultralow values below $\approx 1.0$ Wm*-1K-1*. Non-perturbative anharmonic corrections to the third-order force-constants are important in its calculation as it happens in SnSeAseginolaza et al. (2019). There is a clear anisotropy between in-plane and out-of-plane thermal conductivities. The similarity of the power factors and the lattice thermal conductivities of SnSe and SnS suggest that the latter may be an earth abundant efficient thermoelectric material and motivate more experimental effort to measure its $ZT$ in the high-temperature phase.

This article is organized as follows. In section II we briefly review the theoretical background for the calculation of anharmonic phonons, thermal conductivity, and electronic transport properties. In section III we specify the computational details. In section IV we discuss the results of our work. Conclusions are given in section V.

II Theoretical Background

II.1 Electronic transport properties

Within the semiclassical Boltzmann transport theoryScheidemantel et al. (2003) the electrical conductivity and the Seebeck coefficient can be calculated respectively as

[TABLE]

where $e$ is the electron charge, $\mu$ the chemical potential, $f(T,\mu,\varepsilon)$ the Fermi-Dirac distribution function, and $\Sigma(\varepsilon)$ the transport distribution function. The latter is defined as

[TABLE]

where $\Omega$ is the unit cell volume, $N_{\mathbf{k}}$ the number of $\mathbf{k}$ points in the sum, and $\varepsilon_{n\mathbf{k}}$, $\mathbf{v}_{n\mathbf{k}}$ and $\tau^{e}_{n\mathbf{k}}$ are, respectively, the energy, Fermi velocity and relaxation time of the electronic state with band index $n$ and crystal momentum $\mathbf{k}$. Our goal here is to compare the power factors $P_{F}(T,\mu)=\sigma(T,\mu)S^{2}(T,\mu)$ of SnSe and SnS coming from their different band structure without explicitly calculating the electronic relaxation times. We thus assume that $\tau^{e}_{n\mathbf{k}}=\tau^{e}$ is just the same constant for both compounds. In these conditions it is easy to see from Eqs. (2)-(4) that the power factor is proportional to $\tau^{e}$. In the following we will limit ourselves to the analysis of $P_{F}(T,\mu)/\tau^{e}$, which only depends on the band structure of the compounds.

II.2 Free energy of strongly anharmonic crystals

We study the vibrational properties of SnS within the Born-Oppenheimer (BO) approximation, in which the Hamiltonian $H$ that determines the dynamics of the ions consists of the ionic kinetic energy and the BO potential $V(\mathbf{R})$. $\mathbf{R}$ denotes $R^{\alpha s}(\mathbf{l})$ in component free notation, which specifies the atomic configuration of the crystal. $\alpha$ is a Cartesian direction, $s$ labels an atom within the unit cell, and $\mathbf{l}$ a lattice vector. From now on, we will use a single composite index $a=(\alpha,s,\mathbf{l})$ and bold letters to indicate quantities in component-free notation. We will keep this composite index for Fourier transformed components adding a bar, $\bar{a}$, to distinguish that in this case $\bar{a}$ just denotes a Cartesian index and an atom in the unit cell.

As it will be shown below and as it happens in $Cmcm$ SnSeAseginolaza et al. (2019); Skelton et al. (2016); Dewandre et al. (2016), the harmonic approximation collapses for $Cmcm$ SnS. Truncating the Taylor expansion of $V(\mathbf{R})$ for this phase at second order and diagonalizing the resulting harmonic force-constants $\bm{\phi}$ large imaginary frequencies are obtained. This makes the calculation of any thermodynamic and transport property involving phonons impossible at the harmonic level. We overcome this problem by solving the ionic Hamiltonian within the SSCHA, a variational method that includes anharmonic effects at a non-perturbative level in the calculation of the vibrational free energyErrea et al. (2014); Bianco et al. (2017); Monacelli et al. (2018).

The SSCHA performs a variational minimization of the free energy with respect to a trial density matrix $\rho_{\mathcal{H}}$ that solves an auxiliary harmonic Hamiltonian

[TABLE]

where $P$ is the kinetic energy and $M_{a}$ the atomic mass of atom $a$. The variational parameters in the minimization are the $\mathbf{\Phi}$ force-constants, which should not be confused with the harmonic force-constants $\bm{\phi}$, and the $\bm{\mathcal{R}}$ positions. The $\bm{\mathcal{R}}$ positions are referred as the centroid positions, i.e., the most probable atomic positions. The $\mathbf{\Phi}$ force-constants are related to the broadening of the ionic wave functions around the centroid positions. At the minimum, the SSCHA yields a free energy $F$ that takes into account anharmonicity without approximating the BO potential. The minimization can be performed by calculating atomic forces and energies in stochastic atomic configurations in supercells using importance sampling and reweighting techniquesErrea et al. (2014); Bianco et al. (2017); Monacelli et al. (2018). The supercell atomic configurations are created according to the probability distribution function related to $\rho_{\mathcal{H}}$. Since the BO energy landscape is sampled stochastically, the SSCHA method does not use any fit or approximation on the $V(\bm{R})$. It is, therefore, unbiased by the starting point.

II.3 Free energy Hessian and second-order phase transition

In a displacive second-order phase transition, at high temperature the free energy $F$ has a minimum in a high-symmetry configuration ($\bm{\mathcal{R}}_{hs}$), but, on lowering the temperature, $\bm{\mathcal{R}}_{hs}$ becomes a saddle point at the transition temperature $T_{c}$. Therefore, the free energy Hessian evaluated at $\bm{\mathcal{R}}_{hs}$, $\partial^{2}F/\partial\bm{\mathcal{R}}\partial\bm{\mathcal{R}}|_{\bm{\mathcal{R}}_{hs}}$, at high temperature is positive definite but it develops one or multiple negative eigendirections at $T_{c}$. The SSCHA free energy Hessian can be computed by using the analytic formulaBianco et al. (2017)

[TABLE]

with

[TABLE]

Here $\left\langle\right\rangle_{\rho_{\mathcal{H}}}$ denotes the quantum statistical average taken with the density matrix $\rho_{\mathcal{H}}$. All these averages are evaluated here stochastically as described in Ref. Bianco et al., 2017. The $\overset{(n)}{\mathbf{\Phi}}$ non-perturbative $n$-th order force-constants should not be confused with the $n$-th order perturbative force-constants $\overset{(n)}{\bm{\phi}}$, which are calculated as derivatives of the BO potential at a reference position [math] and not as quantum statistical averages:

[TABLE]

In Eq. (6) the value $z=0$ of the fourth-order tensor $\mathbf{\Lambda}(z)$ is used. For a generic complex number $z$ it is defined, in components, by

[TABLE]

with $\tilde{\Omega}_{\mu}^{2}$ and $e_{\mu}^{a}$ the eigenvalues and corresponding eigenvectors of

[TABLE]

respectively. In Eq. (9)

[TABLE]

where $n_{B}(\omega)=1/(e^{\beta\hbar\omega}-1)$ is the bosonic occupation number. Evaluating through Eq. (6) the free energy Hessian at $\bm{\mathcal{R}}_{hs}$ and studying its spectrum as a function of temperature, we can predict the occurrence of a displacive phase transition and estimate $T_{c}$. This technique has been successful to study phase-transition temperatures in high-pressure hydrides, monochalcogenides, and transition metal dichalcogenides undergoing charge-density wave transitionsRibeiro et al. (2018); Bianco et al. (2018, 2019).

II.4 Dynamical properties of solids and phonon frequencies

As shown in Ref. Bianco et al., 2017, even if the SSCHA is a ground-state theory, it is possible to formulate a valid ansatz in order to calculate dynamical properties of crystals such as phonon spectral functions. The one-phonon Green function $\mathbf{G}(z)$ for the variable $\sqrt{M_{a}}(R_{a}-\mathcal{R}_{a})$ can be calculated as

[TABLE]

With this definition, in the static limit the Green function becomes the dynamical matrix obtained with the free energy Hessian: $\mathbf{G}^{-1}(0)=-\mathbf{D}^{(F)}$, with

[TABLE]

We will label with $\omega_{\mu}$ the eigenvalues of $\mathbf{D}^{(F)}$. The SSCHA self-energy is given by

[TABLE]

where $M_{ab}=\delta_{ab}M_{a}$ is the mass matrix. We have explicitly verified that neglecting $\overset{(4)}{\bm{\Phi}}$ in Eq. (6) has a completely negligible impact on the eigenvalues of $\mathbf{D}^{(F)}$. We consistently neglect $\overset{(4)}{\bm{\Phi}}$ in Eq. (12) as well. This reduces the SSCHA self energy to the so-called bubble self energy, namely

[TABLE]

The cross section in an inelastic, e.g. neutron, experiment is proportional to the spectral function $\sigma(\omega)=-\omega\mathrm{T}rIm\mathbf{G}(\omega+i0^{+})/\pi$Cowley (1968). Its peaks signal the presence of collective vibrational excitations (phonons) having certain energies and linewidth. In order to recognize the contribution of each phonon mode to this spectral function, we first take advantage of the lattice periodicity and Fourier transform the spectral function and the self energy, and second we neglect the mixing between phonon modes and assume that $\mathbf{\Pi}(z)$ is diagonal in the basis of the eigenvectors:

[TABLE]

The cross section is then given by

[TABLE]

In Eqs. (16) and (17) $\tilde{\Omega}_{\mu}^{2}(\mathbf{q})$ and $e_{\mu}^{\bar{a}}(\mathbf{q})$ are, respectively, the eigenvalues and eigenvectors of $\mathbf{D}^{(S)}(\mathbf{q})$, the Fourier transform of Eq. (10).

The cross section calculated as in Eq. (17) does not have any given lineshape. However, when $\Pi_{\mu}(\mathbf{q},\omega)$ is small compared to $\tilde{\Omega}_{\mu}^{2}(\mathbf{q})$, it is justified to approximate $\Pi_{\mu}(\mathbf{q},\omega)\sim\Pi_{\mu}(\mathbf{q},\tilde{\Omega}_{\mu}(\mathbf{q}))$, which turns $\sigma(\mathbf{q},\omega)$ into a sum of Lorentzian functions. In this Lorentzian approximation the peaks appear at the $\Omega_{\mu}(\mathbf{q})$ phonon frequencies, with

[TABLE]

and the linewidths are proportional to $Im[\Pi_{\mu}(\mathbf{q},\tilde{\Omega}_{\mu}(\mathbf{q}))]$.

II.5 Thermal conductivity

We calculate the lattice thermal conductivity within the single mode relaxation time approximation (SMA)Khitun and Wang (2001) making use of the eigenvalues and eigenvectors of $\mathbf{D}^{(S)}(\mathbf{q})$ (as it will be shown below it is not possible at the harmonic level due to the instabilities obtained) as well as the non-preturbative third-order force-constants $\overset{(3)}{\mathbf{\Phi}}$. In the SMA the lattice thermal conductivity is written as followsPaulatto et al. (2013):

[TABLE]

where, for the phonon mode $\mu$ with momentum $\mathbf{q}$, $c_{\mu}^{\alpha}(\mathbf{q})$ is the Cartesian component $\alpha$ of its lattice group velocity and $\tau_{\mu}(\mathbf{q})$ its lifetime. $N_{\mathbf{q}}$ is the number of $\mathbf{q}$ points used in the sum. The Bose-Einstein occupation of each mode is given by the Boltzmann Transport Equation (BTE) and the inverse phonon lifetime (with $\gamma_{\mu}(\mathbf{q})$ the half width at half maximum) is calculated asPaulatto et al. (2013)

[TABLE]

with $\mathbf{q}+\mathbf{q}^{\prime}+\mathbf{q}^{\prime\prime}=\mathbf{G}$, $\mathbf{G}$ being a reciprocal lattice vector. Here $\overset{(3)}{\Phi}_{\mu\nu\eta}(\mathbf{q},\mathbf{q}^{\prime},\mathbf{q}^{\prime\prime})$ is the third order force-constants matrix written in the space of the normal modes:

[TABLE]

where $\overset{(3)}{\Phi}_{\bar{a}\bar{b}\bar{c}}(\mathbf{q},\mathbf{q}^{\prime},\mathbf{q}^{\prime\prime})$ are the Fourier transformed non-perturbative third-order force-constants. We also calculate the thermal conductivity with the perturbative third-order force-constants by substituting the non-perturbative $\overset{(3)}{\mathbf{\Phi}}$ by the perturbative $\overset{(3)}{\bm{\phi}}$ in Eqs. (20) and (21).

III Computational Details

We calculate the electronic bands using ab initio Density Functional Theory (DFT) calculations within the local density approximation (LDA)Perdew and Zunger (1981) and the generalized gradient approximation in the Perdew Burke Ernzerhof (PBE) parametrizationPerdew et al. (1996) as implemented in the Quantum ESPRESSOGiannozzi et al. (2009, 2017) software package. Harmonic phonons and perturbative third-order force-constants $\overset{(3)}{\bm{\phi}}$ are calculated using Density Functional Perturbation TheoryBaroni et al. (2001); Paulatto et al. (2013). We use projector augmented waveBlöchl (1994) (PAW) pseudopotentials that include $5s^{2}$ $5p^{2}$ $4d^{10}$ electrons in the case of Sn and $3s^{2}$ $3p^{4}$ in the case of S or Se. For the perturbative third-order force-constants we use norm-conserving pseudopotentials which were shownAseginolaza et al. (2019) to provide very similar third-order force-constants compared to the PAW result. A $16\times 16\times 16$ sampling of the first Brillouin zone of the primitive cell and an energy cutoff of $70$ Ry are employed for the DFT self-consistent calculation. For the electronic transport calculations we use the Boltztrap software packageMadsen and Singh (2006). For the sum in Eq. 4 we perform a non self-consistent DFT calculation in a $30\times 30\times 30$ sampling of the first Brillouin zone. We use experimental lattice parameters at the transition temperature as we got better agreement with experiments for SnSe in a previous workAseginolaza et al. (2019). The experimental lattice parameters taken from Refs.Adouby (1998); Chattopadhyay et al. (1986) are $a=22.13$ a0, $b=8.13$ a0, $c=8.13$ a0 for SnSe and $a=21.69$ a0, $b=7.84$ a0, $c=7.84$ a0 (a0 is the Bohr length) for SnS. The structures of the high temperature $Cmcm$ and low temperature $Pnma$ phases are shown in Figure 1.

Anharmonic phonons and non-perturbative third-order force-constants are calculated within the SSCHA using a $2\times 2\times 2$ supercell. For the SSCHA calculation we use forces calculated within DFT. Once we get the anharmonic force-constants, we substract the harmonic ones and interpolate the difference to a $6\times 6\times 6$ grid. Then, we add this interpolated difference to the harmonic dynamical matrices in a $6\times 6\times 6$ grid, which yields anharmonic force-constants in a fine grid. By Fourier interpolating the latter we can calculate phonon frequencies at any point in the Brillouin zone. We impose the acoustic sum rule to the third-order force-constants with an iterative method prior to their Fourier interpolationPaulatto et al. (2013); Aseginolaza et al. (2019). The lattice thermal conductivity is calculated with Eq. (19) using a $10\times 10\times 10$ grid. For the calculation of the phonon linewidths we use a $20\times 20\times 20$ mesh in Eq. (20) with a Gaussian smearing of 1 cm*-1* for the Dirac deltas.

IV Results and Discussion

IV.1 Electronic transport

Figure 2 (a) shows the electronic band structures of SnS and SnSe in the high symmetry phase.

It shows that the electronic properties of these materials are very similar because their electronic band structures are basically the same as expected for isoelectronic compounds with the same atomic structure. The major difference is that the indirect (the conduction and valence bands that constitute the gap are denoted with an arrow in Figure 2 (a)) energy gap ($0.45$ eV for SnSe and $0.7$ eV for SnS) is bigger in the case of SnS, in agreement with experimentsVidal et al. (2012); Zhao et al. (2014) and previous calculationsGuo et al. (2015). As expected, the calculated electronic gaps within LDA underestimate the experimental values ($0.86-0.948$ eV for SnSe and $0.9-1.142$ eV for SnS). Using these band structures we have calculated the Seebeck coefficient, which within the approximation of a constant electronic relaxation time it is independent of it, and the electrical conductivity over the electronic relaxation time $\sigma/\tau^{e}$. The Seebeck coefficient is very similar for both materials, but $\sigma/\tau^{e}$ is slightly larger in the case of SnSe due to the smaller electronic gap. Using these two quantities we have calculated $P_{F}/\tau^{e}$, shown in Figure 2 (b). As we can see, $P_{F}/\tau^{e}$ is very similar for both materials, but slightly higher in the case of SnSe. As we can see, $P_{F}/\tau^{e}$ increases as temperature increases and the difference between SnSe and SnS is less than $5\%$ at $1000$ K. These results make clear that regarding the electronic transport properties these two materials are very similar in the high temperature phase provided that the relaxation time for the electrons is similar for both materials, which is expected for isoelectronic and isostructural compounds.

IV.2 $Pnma$ to $Cmcm$ phase transition

As it was already pointed outChattopadhyay et al. (1986); Aseginolaza et al. (2019), symmetryOrobengoa et al. (2009); Perez-Mato et al. (2010) dictates that it is possible to have a second-order phase transition between the $Cmcm$ and $Pnma$ phases. The transition is dominated by the distortion pattern associated to a non-degenerate mode ($Y_{1}$) at the zone border $Y$ point. This meansAseginolaza et al. (2019) that, in a second-order displacive phase transition scenario, the transition temperature $T_{c}$ is defined as $\partial^{2}F/\partial Q^{2}(T=T_{c})=0$ where $Q$ is the order parameter that transforms the system continuously from the $Pnma$ ($Q\neq 0$) to the $Cmcm$ ($Q=0$) phase. As the distortion is dominated by the $Y_{1}$ phonon, $\partial^{2}F/\partial Q^{2}(T)$ is proportional to $\omega_{Y_{1}}^{2}(T)$, which we can calculate diagonalizing $\mathbf{D}^{(F)}$.

Figure 3 shows $\omega_{Y_{1}}^{2}(T)$ within the LDA and PBE approximations.

As in the case of SnSeAseginolaza et al. (2019), the second derivative of the free energy is positive at high temperatures and decreases lowering the temperature. For both approximations, it becomes negative at the critical temperature $T_{c}$, which means that the $Pnma$ phase is not any longer a minimum of the free energy and the structure distorts adopting the $Pnma$ phase. $T_{c}$ strongly depends on the approximation of the exchange-correlation functional: it is $600$ K for LDA and $465$ K for PBE. Our LDA calculation agrees better with the experimental value, around $900$ KChattopadhyay et al. (1986). We associate the discrepancy between LDA and PBE to the different pressures obtained in the transition direction, $P_{zz}$. In fact, as shown in the case of SnSeAseginolaza et al. (2019), $T_{c}$ depends strongly on the pressure in this $z$ direction. The pressure in Figure 3 includes anharmonic vibrational effects on the energy following the procedure outlined in Ref. Monacelli et al., 2018. For the same lattice parameter LDA displays a much smaller pressure, as generally LDA predicts smaller lattice volumes than PBE. The underestimation with respect to experiments may be attributed to the small supercell size used for the SSCHA calculations ($2\times 2\times 2$). Even if experimentally $T_{c}$ is around $100$ K higher in SnS than in SnSe, our LDA calculations give basically the same transition temperature for both materials as $T_{c}=616$ K in SnSe according to our previous calculationsAseginolaza et al. (2019). However, within PBE SnSe does show a lower transition temperature since $T_{c}=299$ K for SnSe KAseginolaza et al. (2019).

IV.3 Anharmonic phonons

Figure 4 (a) compares the harmonic phonon spectrum with the anharmonic one calculated within the Lorentzian approximation at $800$ K within the LDA. In the anharmonic spectrum shown the phonon energies correspond to the $\Omega_{\mu}(\mathbf{q})$ values of Eq. (18). The linewidth obtained in the Lorentzian approximation is also shown.

The phonon spectrum suffers from a huge anharmonic renormalization. The harmonic spectrum shows broad instabilities, which are stabilized by anharmonicity. The $Y_{1}$ mode is unstable below the transition temperature, but it is stabilized after the transition. By having a look at the the phonon linewidths, we can see that two modes at the $\Gamma$ point ($\Gamma_{1}$ and $\Gamma_{2}$) not only suffer a strong anharmonic renormalization, but they also have a large linewidth compared to the rest of the modes in the first Brillouin zone. These modes describe optical in-plane atomic displacements (see Figure 1, $\Gamma_{2}$ has the same atomic displacements as $\Gamma_{1}$ but in the other in-plane direction), which are the same atomic displacements of $Y_{2}$ and $Y_{3}$ at the point $Y$ with a different periodicity due to the different momentum. The $Y_{2}$ and $Y_{3}$ in-plane modes also show a very large linewidth. On the contrary, the linewidth of mode $Y_{1}$ is not so large even if it is responsible for the phase transition (see Figure 1).

In strongly anharmonic materialsDelaire et al. (2011); Li et al. (2014); Ribeiro et al. (2018, 2018); Li et al. (2014); Paulatto et al. (2015); Bianco et al. (2018); Aseginolaza et al. (2019), the phonon spectral functions $\sigma(\mathbf{q},\omega)$ show broad peaks, shoulders, and satellite peaks that cannot be captured by the simple Lorentzian picture. In Figure 4 (b) and (c) we show the spectral function keeping the full frequency dependence on the self-energy (see Eq. (17)). The calculation is done for the $\Gamma$ and $Y$ points. The great majority of the modes describe a Lorentzian shape. However, the modes with a large linewidth within the Lorentzian approximation (see Figure 4 (a)) are those that clearly deviate from the Lorentzian profile ($\Gamma_{1}$, $\Gamma_{2}$, $Y_{2}$, $Y_{3}$). This non-Lorentzian shape makes clear that these modes are strongly anharmonic and the frequency dependence of the self-energy is crucial to account for their spectral function. In this case, as we can see in in Figures 4 (b) and (c), the non-Lorentzian shapes of the strongly anharmonic modes do not create appreciable shoulders or satellite peak in the total spectral function, however, their contribution is far from trivial.

IV.4 Lattice thermal transport

In Figure 5 (a) we show the lattice thermal conductivity of $Cmcm$ SnS as a function of temperature calculated using $\overset{(3)}{\bm{\Phi}}$ and $\overset{(3)}{\bm{\phi}}$ for solving the BTE within the SMA. We recall that $\overset{(3)}{\bm{\Phi}}$ are non-perturbative third-order force-constants calculated using Eq. 7 and $\overset{(3)}{\bm{\phi}}$ are perturbative third-order force-constants calculated using Eq. 8. In Figure 5 (b) we show the lattice thermal conductivities of $Cmcm$ SnS and SnSe using $\overset{(3)}{\bm{\Phi}}$.

We can see that the non-perturbative calculation using $\overset{(3)}{\bm{\Phi}}$ is lower than the perturbative one using $\overset{(3)}{\bm{\phi}}$ for the three Cartesian directions. This result makes clear that the non-perturbative anharmonicity is very important to calculate the thermal conductivity in this kind of thermoelectric materials. By looking at the values of the lattice thermal conductivity we can see that both materials show very similar ultralow values, below $\approx 1.0$ Wm*-1K-1*. In-plane results are slightly higher for SnSe and out-of-plane calculations higher for SnS, in agreements with another calculationGuo et al. (2015) where the thermal conductivities of SnS and SnSe for the low-temperature $Pnma$ phase are calculated with harmonic phonons and perturbative third-order force-constants. Theoretical calculations following the same procedure also show that the thermal conductivities of $Pnma$ SnSe and SnS are very similarSkelton et al. (2017, 2016), in agreement with experimentsZhao et al. (2014); Tan et al. (2014). Our calculations confirm that in the high-temperature $Cmcm$ phase the thermal conductivity of these two compounds is also very similar. Both materials show a clear anisotropy between in-plane and out-of-plane calculations in agreement with experimental resultsIbrahim et al. (2017) for the low-temperature phase close to the phase transition.

V Conclusions

In conclusion, we have calculated the electronic and vibrational transport properties of $Cmcm$ SnS using first-principles calculations. We have seen that the electronic transport properties of SnS and SnSe are comparable and that a similar power factor is expected for these isoelectronic and isostructural compounds. As in the case of SnSe, SnS suffers a second-order phase transition from the $Cmcm$ to the $Pnma$ phase driven by the collapse of a zone border phonon. We have also seen that SnS shows a strongly anharmonic phonon spectrum. Many phonon modes have a very large linewidth and show non-Lorentzian profiles in the spectral function. Finally, we have calculated the lattice thermal conductivity of $Cmcm$ SnS and we have seen that nonperturbative anharmonicity substantially corrects the third order force-constants. The thermal conductivity of both materials is very similar and ultralow. Therefore, by comparing the electronic and vibrational transport properties of SnS and SnSe in the $Cmcm$ high-temperature phase, we conclude both should be good thermoelectrics. Thus, we suggest that SnS may be an earth-abundant very efficient high-temperature thermoelectric material. This work motivates more experimental effort in this regime for its characterization.

VI ACKNOWLEDGMENTS

Financial support was provided by the Spanish Ministry of Economy and Competitiveness (FIS2016-76617-P); and the Department of Education, Universities and Research of the Basque Government and the University of the Basque Country (IT756-13). U.A. is also thankful to the Material Physics Center for a predoctoral fellowship. Computer facilities were provided by the Donostia International Physics Center (DIPC), PRACE (2017174186) and DARI (A0050901202).