Convergence of separate orbits for enhanced thermoelectric performance of layered ZrS2

Minimizing the band splitting energy to approach orbital degeneracy has been shown as a route to improved thermoelectric performance. This represents an open opportunity in some promising layered materials where there is a separation of p orbitals at the valence band edge due to the crystal field splitting. In this work, using ab initio calculations and semiclassical Boltzmann transport theory, we try to figure out how orbital degeneracy influences the thermoelectric properties of layered transition-metal dichalcogenide ZrS2. We tune the splitting energy by applying compressive biaxial strain, and find out that near-degeneration at the Γ point can be achieved for around 3% strain. As expected, the enhanced density-of-states effective mass results in an increased power factor. Interestingly, we also find a marked decline in the lattice thermal conductivity due to the effect of strain on phonon velocities and scattering. The two effects synergetically enhance the figure of merit. Our results highlight the convenience of exploring this optimization route in layered thermoelectric materials with band structures similar to that of ZrS2.


Introduction
Searching for improved thermoelectric materials is a crucial requirement for thermoelectric energy recovery. The efficiency of thermoelectric conversion is governed by the dimensionless figure of merit s k = ZT S T , 2 where S is the Seebeck coefficient, s its electrical conductivity, k its thermal conductivity including both electronic k ( ) e and lattice k ( ) l contributions [1,2], and T the absolute temperature. Enhancing ZT involves working towards a high power factor s ( ) S 2 and a low thermal conductivity, specifically a low lattice thermal conductivity. However, due to correlations between these coefficients, it is a great challenge to optimize both at the same time [1]. A number of strategies have been tried with different degrees of success [3][4][5][6][7][8]. A majority of efforts have focused on reducing lattice thermal conductivity by enhancing phonon scattering through nanostructuring or alloying [3][4][5]. Band engineering, which can markedly increase the Seebeck coefficient while maintain a high electrical conductivity, has also been proposed as a way to enhance the power factor [6][7][8].
Trying to achieve band degeneracy is one of the widely used approaches to band engineering. The idea is to place the extrema of multiple conduction or valence bands within a few k B T of each other [8][9][10]. The degenerate bands cause a localized peak in the density of states, thereby enhancing the Seebeck coefficient. This has little influence on the carrier mobility, meaning that enhanced power factors can be obtained. Tang et al achieved a high ZT>1 in CoSb 3 skutterudites by converging the secondary conduction band with the primary light conduction bands at high temperature (800 K) [11]. The creation of solid solutions has also been demonstrated as an effective approach towards band degeneracy. Liu et al obtained an enhanced ZT∼1.3 in n-type Mg 2 Si 1-x Sn x solid solutions, and found that the Seebeck coefficient and power factor reach their maximum when x=0.7, since the two light and the heavy conduction bands become effectively degenerate at X point as confirmed by ab initio calculations [12].
Likewise, orbital degeneracy has been recently proposed as a way to optimize the thermoelectric performance of some layered materials, e.g. CaAl 2 Si 2 -type Zintl compounds and transition-metal dichalcogenides MX 2 [13,14]. For compounds with cubic symmetry, the p orbitals at the valence band edge are threefold degenerate at the G point. On the other hand, in layered materials p z orbital is separated from p x and p y due to the crystal field splitting. Accordingly, a non-degenerate G( ) p z band and a doubly degenerate G( ) p x y , band are formed. The splitting energy is defined as x y z , the energy difference between the two bands at the G point. Zhang et al demonstrated that excellent electronic transport performance can be realized by manipulating the bands to bring this splitting energy closer to zero, which was done by means of solid solution and biaxial strain [13,14]. As for biaxial strain, tensile strain is more effective for a positive D while negative D requires compressive strain. In spite of this successful demonstration, the effects of orbital degeneracy on electron and phonon transport are still not fully understood. We note that a recent work reported on the potential thermoelectric performance of layered ZrS 2 , where there is a separation of p orbitals at the valence band edge [15]. However, there have been no studies on the effect of strain on thermoelectric properties for bulk layered ZrS 2 . In this paper, we try to use a biaxial strain to control the orbital degeneracy of layered ZrS 2 . We then explore the effect on its thermoelectric transport coefficients, with an eye to improving the thermoelectric performance of this and related materials.

Methodology
The energy band structures are calculated within the framework of density functional theory (DFT) using the projector-augmented-wave formalism [16] and the Perdew-Burke-Ernzerhof exchange and correlation functional [17], as implemented in VASP [18]. During the geometry optimization, both position and lattice constant are fully relaxed until the forces on all atoms are lower than 0.01 eV Å -1 . The plane-wave energy cutoff is set at 400 eV, and the irreducible Brillouin zone is sampled with a 15×15×9 Monkhorst-Pack mesh. Considering the van der Waals interactions between atomic layers, we also include DFT-D2 function in the calculations [19]. The transition-metal dichalcogenide ZrS 2 has a 1T-CdI 2 type structure with the space group P-3m1, and our optimized lattice constants for unstrained ZrS 2 are a=3.65 Å and c=5.82 Å, which are in good agreement with previous results [15,20]. The biaxial strain ε is defined as (a 0a)/a 0 ×100%, where a and a 0 are the in-plane strained and unstrained lattice parameters, respectively.
We calculate the thermoelectric transport coefficients in the framework of the Boltzmann transport equation (BTE) within the relaxation time approximation (RTA) [21][22][23], which describe the deviation from equilibrium of the carrier occupancies under the action of an external field. In the steady state, they can be written as where  F represents the external forces. The electronic transport coefficients can be derived from the calculated electronic structure, and the key point is to find so-called transport distribution by solving BTE, which is expressed as once it is defined, the Seebeck coefficient S and electrical conductivity s are given by, where m and f 0 are the chemical potential and equilibrium Fermi function, respectively. Next, it is worthwhile to note that we only need to consider the temperature gradient as external field when discussing the phonon BTE.
Here, we adopt the method provided by Li et al to calculate the lattice thermal conductivity [24][25][26]. To fully solve phonon BTE, they define where l  v is the group velocity of mode l, and D l in the dimension of velocity is a correction of deviation to RTA prediction. t l 0 is the relaxation time of mode l and computed as [26].
are the three-phonon scattering rates and G ll¢ is the scattering probabilities from isotopic disorder. Finally, the lattice thermal conductivity can be obtained in term of l  F as [26], where W is the volume of the unit cell, α and β denote the Cartesian components.

Results and discussion
In this section, we would like to first present the calculated whole thermoelectric coefficients under different strains. It is crucial to find out how the thermoelectric coefficients change with the increased strain, and whether an enhanced power factor and a depressed lattice thermal conductivity will be defined just as our prediction. As stated above, the strain is in fact an approach used to manipulate the level of the orbital degeneracy, thus, the connection between orbital degeneracy and thermoelectric coefficients will be explained later from band structures and related mechanisms. The layered transition-metal dichalcogenide ZrS 2 has a negative crystal field splitting energy D < 0, which means the light hole band G( ) p z is above the heavy hole band G( ) p x y , at G point. The splitting energy is about 0.2 eV in ZrS 2 . Some compounds instead have a positive splitting energy such as TiS 2 [13]. It is desirable for optimizing the electronic transport performance that the two separated bands are nearly degenerated as proposed by Zhang et al [13], and compressive biaxial strain is tested more effective for a negative D. We thereby try a series of biaxial strain ranging from 0% to 3% in order to control the magnitude of the splitting energy. In the electronic transport calculations, we consider an experimental carrier concentration (1.4×10 18 cm −3 ) of intrinsic ZrS 2 [27], which allows us to calculate the thermoelectric coefficients with respect to temperatures. Figure 1 shows the calculated dependence of p-type electronic transport coefficients on temperatures under different strains. One can find there is a critical temperature around 550 K, where the Seebeck coefficient and the electrical conductivity exhibit a nearly opposite tendency. This is because of their inversely proportional relation with carrier density [1]. In higher temperatures, the excited carriers result in the rapidly decrease of Seebeck coefficient and increase of electrical conductivity. Clearly, the strain has a significant effect on the magnitude of electronic transport coefficients. The most important phenomenon confirmed is that the Seebeck coefficient and electrical conductivity increase simultaneously as the increase of strain, which lead to the enhanced power factors as can be seen in figure 1(c). Compared with the same phenomenon in literatures [6,[8][9][10], this is because of the possibly increased level of band degeneracy. Here, the orbital splitting energy in ZrS 2 would be expected to decrease as the increasing strain. Thereby, the electronic transport performance can be optimized as the increased orbital degeneracy, and further explanations would be addressed later. In figure 1(d), the electrical thermal conductivity is calculated from the Wiedemann-Franz law [28]: k s = L T, e which shows similar trend to electrical conductivity due to their proportional relation.
It is well known that the thermoelectric performance can be optimized by turning the carrier concentrations. We therefore calculate the power factors with respect to carrier concentration at 525 K as shown in figure 2. The black line in inserted figure corresponds to the results in figure 1(c) with the fixed carrier concentration. It is found the power factors can be enhanced tenfold by optimizing the carrier concentration through doping. On the other hand, the strain also plays a positive role in enhancing the power factor, and significant effect can be achieved as the increase of carrier concentrations. Clearly, the strain has a strong effect on the magnitude of electronic transport coefficients, which is supported by an increased density of states effective mass arising from the orbital degeneracy, as we will discuss later. It is worthwhile to note that our calculated electronic transport coefficients for unstrained ZrS 2 are well consistent with previous results [15].
For solving phonon BTE, the harmonic and anharmonic interatomic force constants are needed, which allow us to define the transport features of phonons. These force constants are computed on a well converged supercell of 5×5×3 through ab initio calculations. The harmonic forces are obtained in Phonopy package [29]. For the anharmonic force constants, we impose a tested cutoff so as to take the fourth nearest neighbors into consideration. We are now able to explore the strain effect on phonon transport.
As phonon spectrum can reflect most of the phonon properties, we thus calculate the phonon spectrum at different strains as presented in figure 3. The result of unstrained ZrS 2 is very consistent with previous calculations [15], characterized by a marked dispersion of acoustic modes at G point and small hybridization of longitudinal acoustic phonons with low-lying optical phonons. Moreover, ZrS 2 exhibits a high phonon frequency with the allowed maximum frequency larger than 350 cm −1 , as compared to the low values in compound with heavy atoms [30][31][32]. The calculated velocities of the two transverses and one longitudinal acoustic modes along G-M are 1771/3812 m s −1 and 5547 m s −1 , respectively. Such high velocities will no doubt lead to the high lattice thermal conductivity. It is obvious to notice that there are some significant changes in phonon spectrum when the compressive biaxial strains are applied. Firstly, the strain markedly distorts the first transverse acoustic branch and lowers its dispersive level, which results in a shift of its velocity from 1771 m s −1 to merely 457 m s −1 at 3% strain. As acoustic phonons contribute nearly ∼80% to the lattice thermal conductivity [33][34][35][36][37], we expect to obtain a depressed phonon transport due to the decrease of the velocities of acoustic branches under strain. On the other hand, although the optical branches are slightly pushed to higher frequencies as the increased strain, the coupling of local optical branches becomes a little stronger, which may intensify the scattering on optical phonons. These strain-induced changes in phonon spectrum are expected to impede the phonon transport and lower the lattice thermal conductivity.
Calculated dependence of lattice thermal conductivity on temperatures is shown in figure 4(a). As expected, the lattice thermal conductivity nearly halved as the strain increases from 0% to 3%. Thus, the biaxial strain is also an effective approach to lower the lattice thermal conductivity due to its distortion on phonon bands and possible enhancement of phonon scattering by strain effect. It is worthwhile to note that we just consider the strain up to 3% for the structural stability, and further increasing strain may shift the acoustic modes into negative frequency. To further address the scattering probability from strain-induced disorder, we calculate the room temperature isotopic scattering rates as plotted in figure 4(b). We are able to find much higher scattering rates of optical phonons than acoustic phonons, which is mainly attribute to the intrinsic high frequency and low velocity of optical phonons. Considering the strain effect, the scattering rates for optical phonons has a significant increase, characterized by the frequency between 7 and 9 THz. It is also found that the scattering center gradually shifts to higher frequency in line with the move and locally increased coupling of optical branches. Although the isotopic scattering is negligible as compared to anharmonic one, the anharmonic scattering nearly maintain its magnitude and just merely shift, i.e. a shift of the average rate from 2.8 ps −1 (0%) to 2.5 ps −1 (3%) under room temperature.
It is necessary to compare the phonon properties with available experimental and theoretical results, as shown in table 1. The calculated frequencies of 254 (cm −1 ) and 323 (cm −1 ) for the A 1g and E g vibration modes at G, respectively, are in good agreement with the experimental values of 251 (cm −1 ) and 333 (cm −1 ) by Raman measurement [38] and also the previous theoretical values of 255.46 (cm −1 ) and 334.78 (cm −1 ) [15] for intrinsic ZrS 2 . However, our calculated lattice thermal conductivity at room temperature of 11.1 W mK −1 is somewhat lower than the previous theoretical value of 17.8 W mK −1 [15], and the reason is possible that the quite small supercell of 3×3×2 they chosen for anharmonic calculation [15]. We test this supercell size and find the phonon spectrum does not thermally stable, thus, the lattice thermal conductivity may be not converged. We  also compare unstrained ZrS 2 with recently reported single layer (SL) ZrS 2 . Although the frequencies are much the same, the lattice thermal conductivity for SL-ZrS 2 [39] are much smaller, this decrease mostly originates from the enhanced phonon scattering rates in low-dimensional structure.
With all transport coefficients available, we can now evaluate the figure of merit ZT. However, due to the RTA, the electrical conductivity and electrical thermal conductivity can only be calculated with respect to the carrier relaxation time, i.e. s t, k t. e It is feasible to use a constant relaxation time by fitting with experimental results [40][41][42]. Here, the value of t is estimated at the experimental carrier concentration by fitting the electrical resistivity at 300 K to the experimental value of 0.3 W cm [27]. This yields a relaxation time of 0.64×10 −14 s. The carrier relaxation time is reasonably assumed to be independent of the strain as the lattice parameters for each step tuned by the biaxial strain is much small [13]. Figure 5(a) shows the evaluated ZT as a function of carrier density at 525 K, enhanced ZT can be achieved as the increase of strain. Specifically, this is related to the enhanced power factor and decreased lattice thermal conductivity discussed above. However, due to the high lattice thermal conductivity, ZrS 2 did not show a promising thermoelectric performance, but it does not matter when we discuss the strain effect on the thermoelectric performance. In figure 5(b), we also plot the ZT values as regards to temperature at optimal carrier concentrations. We see that the ZT values increase with temperature and may possibly exceed 1 at higher temperature, which indicates that ZrS 2 could be potentially applied as thermoelectrics working under temperatures around 1000 K.
To confirm the enhanced electronic transport coefficients as the strain-induced orbital degeneracy at the valence band edge, we calculate the band structures under different strains as shown in figure 6. Relative motion of the light G( ) p z band and the heavy G( ) p x y , band at G point can be obtained as the strength of biaxial strain increases from 0% to 3%. While the position of G( ) p z does not change much as the increasing strain, the G( ) p x y , gradually moves up in energy. Thus, the two separated bands with different orbital composition converge as the increasing strain and almost degenerate at ε=3%. This is precisely consistent with the enhanced electronic transport performance. In a semiconductor consisting of several conduction or valence bands for electronic   where subscripts 1 and 2 refer to the transport contribution of carriers from individual band. In principle, this kind of effect will be greatly enhanced when the energy of the band extrema within ∼2k B T of each other. Thereby, the decreasing splitting energy D leads to the enhanced electronic transport performance. From another point of view, the Seebeck coefficient usually increases with the increasing density of states effective mass, as given by [1], where N v includes the orbital degeneracy. When the energy difference between G( ) p z and G( ) p x y , is less than 2k B T, the degenerate orbits will effectively increase the number of nonequivalent pockets of carriers N v , and therefore the Seebeck coefficient. Since it has no adverse effect on the carrier mobility, exceptionally high power factors can be achieved m µ ( ( ) ) ⁎ m m PF e d 3 2 [9]. Along with the decreased lattice thermal conductivity, an enhanced thermoelectric performance is realized through this kind of strain-induced orbital degeneracy. Consequently, orbital degeneracy is an effective approach to optimize the thermoelectric performance in this kind of layered materials.
Finally, considering the better thermoelectric performance at high temperatures we studied, it is necessary to discuss the thermal stability of ZrS 2 at high temperatures. We perform a Born-Oppenheimer molecular dynamics simulation based on DFT, also implemented in the VASP [18], for ZrS 2 with a 2×2×2 supercell at 1000 and 800 K with a time step of 1 fs. After 3 ps, no structural destruction is found (see figure 7) except for some thermal fluctuations (see figure 8) at 1000 and 800 K. So, layered ZrS 2 is thermal stable below 1000 K. In fact, the experimental synthesized temperature 1023 K [27] also indicated the thermal stability of ZrS 2 at high temperatures.

Summary
We have used the ab initio calculations and semiclassical Boltzmann transport theory to suggest that it is feasible to optimize the thermoelectric performance through orbital degeneracy for the layered transition-metal dichalcogenide ZrS 2 , which can be extend to other systems with similar band features. We manipulate the splitting energy by means of compressive biaxial strain. As a result, the degenerate orbits lead to the enhanced electronic transport coefficients due to the increased density of states effective mass and transport contribution of carriers from individual band. Meanwhile, the lattice thermal conductivity is reduced because of the straininduced distortion of phonon band and the scattering probabilities from isotopic disorder. Combining these, an Figure 6. Calculated band structures with orbital weights, the red squares denote the band with p z states and blue squares denote the p x y , states. The splitting level is characterized by the energy difference between the two bands at G point. enhanced thermoelectric performance is thereby obtained. The present work provides a comprehensive understanding on the role of orbital degeneracy in enhancing the thermoelectric performance, and will stimulate further studies on designing high-performance layered thermoelectric materials.