Unexpectedly high thermoelectric performance of anisotropic Zr2Cl4 monolayer

Recently, the Hf2Cl4-type materials as functional materials have attracted broad interest because of their enormous potential in thermoelectric (TE) applications. However, relevant investigations are still scarce up to now. To explore the Hf2Cl4-type materials with excellent TE properties, we focus on the TE properties of Zr2Cl4 monolayer and calculate the TE parameters based on first-principles calculations and Boltzmann transport equation. Although, as compared to some typical TE materials, it exhibits better heat transport and thus higher lattice thermal conductivity, the figure of merits (ZT) of both p-type and n-type Zr2Cl4 reach an unexpectedly high value of 3.90 and 3.60, respectively, owing to the larger electrical conductivity and higher power factor. Additionally, owing to the prominent difference in electrical conductivity between the x- and y-direction, strong anisotropy in ZT values is observed. Our study reveals that both n-type and p-type Zr2Cl4 monolayers have the potential for future TE applications.


Introduction
With the continuous development of technologies, the demand for energy in human society and various industries, especially renewable and clean energy, is continuously increasing [1][2][3]. The key to solving this issue is improving the efficiency of energy [4,5]. Therefore, numerous thermoelectric (TE) materials, which can efficiently convert waste heat into electricity without air pollution, have been proposed theoretically, paving the way for developing green consumer and industrial products [6,7]. The performance of TE materials is usually described by a figure of merit (ZT), ZT = S 2 σT/κ, where S is the Seebeck coefficient, σ represents the electrical conductivity, κ contains κ e (electronic thermal conductivity) and κ l (lattice thermal conductivity) [8]. Notably, due to the coupling between these TE parameters [9], the ZT value of many TE materials is still below 1 at room temperature [10,11]. Many strategies have been adopted for enhancing the TE performance to date [12][13][14]. For example, Xiao et al incorporated the PbTe phase into the Pb 0.94 Sn 0.06 S matrix and obtained a low lattice thermal conductivity (∼0.61 Wm -1 K -1 ) and an enhanced ZT value (∼1.3) at 923 K in Pb 0.94 Sn 0.06 S-8%PbTe [15]. On the other hand, reducing the dimension had been another method to enhance the ZT value of TE materials since the quantum confinement effect in two-dimensional (2D) materials was found to enhance the power factor (PF) by Hicks and Dresselhaus [16,17] in 1993. From then on, the study on TEs sparked considerable interest in 2D materials. For example, Tang et al [18] investigated the TE transport properties of p-type Bi 2 Te 3 with layered nanostructure and found that there is a significant improvement in the TE ZT (ZT = 1.35 at 300 K). Consequently, it is feasible to explore new monolayer materials as potential candidates for TE applications.
To find excellent TE materials, many methods like experimental synthesis and theoretical prediction have been conducted. In particular, high-throughput calculations for screening promising TE materials have been studied in many groups [19,20]. For example, Chen et al [21] studied the TE properties of more than 48 000 inorganic compounds using highthroughput calculations and compared their calculations with experimental data. They found the computational maximum Seebeck coefficient is consistent with the experimental value. In 2020, Sarikurt et al [22] screened the promising TE materials using high-throughput calculations and proposed that strong anisotropy and prominently high electronic fitness function (EFF) [23] values exist in many 2D layered materials. Based on the investigation by Sarikurt et al, many monolayers were confirmed to be potential TE materials. For instance, the n-type doping monolayer InClSe with a high ZT value (2.82 at 700 K) was found and proved to be a promising candidate for TE application [24]. In addition, Chang et al [25] demonstrate that 2D Ga 2 I 2 S 2 has high anisotropic ZT values and can be regarded as a new potential room-temperature TE material. Qi et al [26,27] investigated and compared the TE performance of Al 2 X 2 Se 2 (X = Cl, Br, I) in recent years. These researches confirmed that high-throughput calculations are reliable for screening materials with high TE performance.
Recently, the Hf 2 Cl 4 monolayer was confirmed by Li et al [28] as an excellent TE material with highly anisotropic ZT value of 3.2 along the x-axis and 1.1 along the y-axis at 700 K. As an isostructural material (prototype WTe 2 ), Zr 2 Cl 4 monolayer (space group: P2_1/m) is considered as a promising TE material because of its high n-type EFF values [22]. Besides, as group VIA elements, Zr and Hf atoms have the similar valence electronic structure and the investigation on the TE performance of Zr 2 Cl 4 has not been reported yet. Therefore, it is meaningful to explore the electrical and thermal transport mechanisms of Zr 2 Cl 4 monolayer.
In this work, we use the density functional theory (DFT) and Boltzmann transport theory to study the TE properties of WTe 2 -type monolayer Zr 2 Cl 4 . In preliminary work, the geometrical parameters and electronic structure of monolayer Zr 2 Cl 4 are calculated. After confirming the thermodynamic stability, we explore the thermal and electronic transport properties of monolayer Zr 2 Cl 4 . Based on the TE parameters (Seebeck coefficient, electrical conductivity, and thermal conductivity) obtained, the ZT thus TE performance of monolayer Zr 2 Cl 4 is estimated. The results show strong anisotropy in the lattice thermal conductivities, electrical conductivities, and ZT values. In addition, an expectedly high ZT value of 3.90 is achieved at 700 K. Our work suggests that Zr 2 Cl 4 monolayer has excellent TE performance and heavy doping effectively enhances its ZT value.

Computational details
In this study, first-principles calculations were performed with the projector augmented wave method as implemented in the Vienna ab-initio simulation package code [29]. The generalized gradient approximation with Perdew-Burke-Ernzerhof (PBE) is selected as the exchange-correlation functional [30]. The Heyd-Scuseria-Ernzerhof (HSE06) hybrid functional [31] is applied to obtain a more accurate electronic band structure and density of states (DOS). The plane-wave cutoff energy and Monkhorst-Pack k-point grid are set to be 500 eV and 15 × 15 × 1, respectively. The convergence energy and force criteria are 10 −5 eV and 10 −2 eV Å −1 , respectively. The supercell we used for calculating phonon dispersion and lattice thermal conductivity is 6 × 5 × 1. A vacuum layer over 20 Å is constructed to avoid the influence of interlayer interactions along the z-direction. The deformation potential (DP) theory is used to approximately evaluate the relaxation time τ and carrier mobility µ according to the formula [32,33]: whereh is Planck constant, T is the temperature, k B is Boltzmann constant, m * , C, and E 1 are the effective mass, elastic constant, and DP constant, respectively. The effective mass is calculated by m * =h 2 / ( . The elastic constant can be obtained by , where E edge is the shift of the band edge, and ∆a is the change of the lattice parameter relative to the equilibrium lattice parameter. Generally, the electron effective mass is divided into band effective mass m * b and DOSs effective mass m * D . Here the effective mass in equation (1) represents m * b . Then the Onsager transport coefficients are calculated by Boltzmann transport equation implemented in the BoltzTraP2 code [34]. Notably, the method we adopt in BoltzTraP2 is adjusting the Fermi level to change the carrier concentration with the band structure unchanged instead of introducing dopant elements. Then we simulate the calculation of Seebeck and electrical conductivity under different carrier concentration. Although the process is different from the experimental case, i.e. the elements doping would give rise to the change of electronic structure, our results can still provide theoretical support in the aspect of doping type and concentration for enhancing the TE performance of monolayer Zr 2 Cl 4 . In addition, DP theory focuses on acoustic phonon dominated scattering process, while optical phonons are not considered. The second-order and third-order interatomic force constants are studied by the Phonopy code [35] and thirdorder.py script [36]. Then, the lattice thermal conductivity κ l is figured out using ShengBTE code [36].

Structural and electronic properties of monolayer Zr 2 Cl 4
The monolayer Zr 2 Cl 4 (space group No. 11, P2_1/m) is obtained by exfoliation from the corresponding bulk (space group No. 31, Pmn2 1 ), as displayed in figure 1(a). To investigate if spin-orbital coupling (SOC) [37] and van der Waals (vdW) [38] interaction is necessary for this work, we also calculate the lattice parameters and band structure of Zr 2 Cl 4 monolayer by standard DFT, DFT + SOC and DFT-D3 methods, as shown in table 1 and figures 2(a)-(d). The optimized lattice parameters of Zr 2 Cl 4 monolayer are a = 3.34 Å and b = 6.23 Å, which agrees with previous results [22]. It can also be seen that there are four different bonds between Zr and Cl atoms. The electronic configuration of Zr 2 Cl 4 monolayer modeled by PBE and HSE06 functional is depicted in figure 1(b). The DOS corresponding to the band structure is calculated by HSE06 functional. Both the conduction band minimum (CBM) and valence band maximum (VBM) are in the Γ-X path, namely, the band gap with HSE06 functional exhibits a direct character of 1.27 eV, which is approximately double than that with PBE functional (0.62 eV). The result is in line with the previous report (0.59 eV) [22], as shown in table 1. Remarkably, the band degeneracy behavior occurs around CBM and VBM in figure 1(b). Specifically, a conduction band valley, namely CB1 in the S-Y path, is numerically close to the CBM, analogous to the case of VBM and VB1, which results in the slope of DOS around CBM and VBM. Such a band degeneracy would enhance the Seebeck coefficient [39,40], as will be further discussed. Similar phenomenon is also observed in GeAs 2 [41] and GeSe [42]. According to the DOS distribution, both CBM and VBM are mainly contributed by the d orbital of Zr atoms. Besides, the results in table 1 and figures 2(a)-(d) indicate the band structure of each method has a small discrepancy, namely, SOC and vdW have negligible influence on the geometrical structure and band structure of Zr 2 Cl 4 monolayer. As a result, we will not consider them in the following calculations.

Structural stability
In order to confirm the stability of Zr 2 Cl 4 monolayer, phonon dispersion, ab-initio molecular dynamics (AIMD), energy hull and elastic constants calculations are carried out. Figure 1(c) shows the phonon spectrum of Zr 2 Cl 4 monolayer. There is no imaginary frequency in the phonon dispersion, indicating Zr 2 Cl 4 monolayer is dynamic stable. Moreover, for the three acoustic modes around the Γ point, both transverse acoustic and longitudinal acoustic modes exhibit a linear trend, while out-of-plane acoustic flexural (ZA) branch shows a quadratic parabola dispersion in the low-frequency range, which is common in 2D materials and ascribed to the rapid weakening of the lateral restoring force [43,44]. The AIMD simulation in figure 3 shows the variation of free energy with time (0-4 ps) and the crystal structures under 300, 500 and 700 K. It can be seen that the average values of the free energy oscillate slightly around a constant after 0.5 ps at each temperature. Furthermore, there is no extra phase forming and no bond breaking, which provides firm evidence that the Zr 2 Cl 4 monolayer is thermally stable under temperatures of 300, 500 and 700 K. Besides, the Zr-Cl energy convex hull in figure 4 reveals that formation energy of Zr 2 Cl 4 is only 0.16 eV/atom above the convex hull, which is comparable to the value of other known 2D materials [45]. In order to confirm the mechanical stability, the elastic constants are calculated, as shown in table 2. It satisfies the Born stability criteria as follows [46]: (1) the second-order elastic stiffness tensor matrix C is positive, (2) all eigenvalues of C are positive, (3) all the leading principal minors of C are positive, and (4) an arbitrary set of minors C are all positive. These results provide the possibility of realizing Zr 2 Cl 4 experimentally.

Thermal transport properties
To explore the TE properties of Zr 2 Cl 4 monolayer, it is mandatory to unambiguously ascertain the lattice thermal conductivity and other correlated thermal parameters. Initially, the convergence test of Q-gird is carried out for obtaining convergent κ l . As shown in figure 5(a), with Q-grid increasing from 7 × 4 × 1 to 63 × 36 × 1, the curves of lattice thermal conductivities along the x-and y-axis become smooth. Therefore, we set Q-grid to 56 × 32 × 1 in this communication. Subsequently, we calculate the lattice thermal conductivities under the temperature from 200 to 700 K. The result in figure 5(b) indicates the κ l decrease with increasing temperature, stemming from the intrinsic phonon-phonon scattering    [47]. It can be found that at 300 K, the κ l along the x-and y-axis are 16.95 and 11.70 W m −1 K −1 , respectively, which is comparable with isostructural Hf 2 Cl 4 monolayer 25 but larger than many 2D TE materials [48,49]. As the temperature rises to 700 K, the κ l along x-axis (yaxis) decreases to 7.23 (4.84) W m −1 K −1 owing to the more frequent phonon scattering at higher temperatures. In addition, an anisotropic ratio (κ lx /κ ly ) is as high as 1.45, which is comparable with ZrSe 2 /HfSe 2 superlattice-monolayer structure (∼1.47) [50] and larger than penta-silicene (1.29) and penta-germanene (1.25) [51]. This behavior can be ascribed to the anisotropic phonon group velocities and mean free path (MFP) according to the formula [16]:

a (Å) b (Å) <Zr-Cl> 1 (Å) <Zr-Cl> 2 (Å) <Zr-Cl> 3 (Å) <Zr-Cl> 4 (Å) Eg
in which C V is the lattice heat capacity, V g is the phonon group velocity and l is the phonon MFP. Figures 5(c) and (d) indicate the relatively low phonon group velocities at 300 K in the y-direction with respect to the x-direction, leading to the reduced lattice heat transport. Figure 5(e) shows that the phonon MFP of ∼1800 nm along the x-axis is larger than that of ∼1500 nm along the y-axis. These account for the larger lattice thermal conductivity along the x-axis. Similar phenomena were also found in many well-known TE materials [50,52]. To further understand the thermal transport properties of Zr 2 Cl 4 monolayer, we calculate the frequency-dependent phonon lifetimes in figure 5(f) and the inset is the mode Grüneisen parameter γ. It can be clearly seen that the phonon lifetimes of acoustic modes are on the order of tens of picoseconds, which is longer than those of optical modes. The longer lifetimes of the acoustic modes suggest their dominant role in the lattice thermal conductivity. The inset shows that except for ZA mode, the averaged values of Grüneisen parameter are nearly zero. Generally, the Grüneisen parameter is used to describe the strength of anharmonic scatterings. As a result, the small γ and long phonon lifetime lead to better phonon transport of Zr 2 Cl 4 monolayer.

Electronic transport properties
In figures 6(a) and (b), we show the calculated Seebeck coefficient of p-and n-type Zr 2 Cl 4 monolayer as a function of  carrier concentration under three temperatures. The negative and positive values of S correspond to n-type and p-type system, respectively, and here we take absolute values |S| in figures 6(a) and (b). Generally, with the decreasing of carrier concentration and the increasing of temperature, the |S| should increase continuously [53]. However, it is remarkable  that, with the falling of carrier concentration to an optimal value, the |S| of n-and p-type Zr 2 Cl 4 monolayer along the xand y-axis show an abnormal behavior named bipolar conduction effect [54] at high temperatures. For instance, we notice that the |S| for the p-type system along the x-axis peaks at a maximum value of 340 µV K −1 at the carrier concentration of 5.07 × 10 12 cm −2 then declines at 700 K. Generally, the bipolar conduction effect naturally arises as a general concern in TE materials with a small band gap, such as Sb 2 Te 3 [55], Bi 2 Te 3 and Bi 2 Te 2 Se [56]. Moreover, it is remarkable that the difference of |S| between p-type and n-type Zr 2 Cl 4 is negligible except for the abnormal part, which can be explained by the formula as follows [53]: where k B , N v , N c , n, p and r are the Boltzmann constant, effective DOS in valence band, effective DOS in conduction band, number of electron carrier, number of hole carrier, and scattering mechanism parameter, respectively. Since figure 1(b) shows that the slope of DOS around VBM is close to that around CBM, under the carrier concentration of 1 × 10 11 cm −2 , the |S| of 635 µV K −1 for p-type system is comparable to that of 636 µV K −1 for n-type system at 300 K in the x-direction. The sharp DOS around CBM and VBM thereby large |S| are ascribed to the band degeneracy behavior as we mentioned in the previous part. We further note that the anisotropy in |S| is very weak. For example, under the carrier concentration of 1 × 10 13 cm −2 , the |S| of p-type Zr 2 Cl 4 monolayer along the x-axis (288 µV K −1 ) is close to the value of 282 µV K −1 along the y-axis. The electrical conductivity obtained from BoltzTrap2 code contains the term of relaxation time τ , which can be figured out according to equation (1) based on DP theory. The results are shown in table 3 and figures 6(c) and (d). As clearly noticed, the electrical conductivity σ is improved with the increasing carrier concentration and decreasing temperature except for the case under the carrier concentration below 1 × 10 12 cm −2 at 700 K, which is due to the bipolar conduction effect as we discussed in the previous part. This behavior is easily understood by the formula [57]: In general, more frequent scattering of carriers occurs at higher temperatures, resulting in shorter relaxation time and lower carrier mobility thereby smaller electrical conductivity [58]. Note also that the σ of p-type Zr 2 Cl 4 is higher than that of n-type Zr 2 Cl 4 in all cases. For example, the σ of 4.39 × 10 4 S m −1 for p-type system is more than twice as high as that of 1.89 × 10 4 S m −1 for n-type system along the x-axis at 300 K at the carrier concentration of 1 × 10 11 cm −2 . Additionally, the electrical conductivity exhibits strong anisotropy, i.e. the σ along the x-axis is larger than that along the y-axis for both n-type and p-type system under constant temperature and carrier concentration. For instance, under the carrier concentration of 1 × 10 14 cm −2 and the temperature of 300 K, the σ for p-type doping along the x-axis (3.62 × 10 7 S m −1 ) is nearly a hundred times of that along the y-axis (8.13 × 10 5 S m −1 ). According to equation (5), these phenomena result from the large discrepancy of carrier mobility between n-and p-type system along the x-and y-axis, as seen in table 3.
To better understand the role of S and σ in the excellent TE performance of layer Zr 2 Cl 4 , we calculate the PF (PF = S 2 σ) shown in figures 6(e) and (f). We can observe that the optimal PF values of p-type Zr 2 Cl 4 are higher than those of n-type Zr 2 Cl 4 in all cases, owing to the larger electrical conductivities. Besides, the strong anisotropy is observed in PF values of n-and p-type system at the three temperatures. In particular, the maximum PF value in the x direction at 300 K reaches 138 mW m −1 K −2 under the hole concentration of 3.59 × 10 13 cm −2 , while that in the y direction is only 2.45 mW m −1 K −2 under the hole concentration of 1.37 × 10 13 cm −2 . The value is significantly larger than bilayer MoS 2 , in which a high PF of 8.5 mW m −1 K −2 was observed by Zhang et al [59] at the room temperature.

TE performance
The electronic thermal conductivity κ e is calculated according to the Wiedemann-Franz law [60,61]: κ e = LσT, where L = 2.45 × 10 −8 WΩ K −2 is the Lorentz constant. For comparison, the electronic thermal at 300 K obtained directly from BoltzTraP2 code: κ e = κ 0 − TσS 2 is calculated [34]. The corresponding comparison between the two formulas in p-and n-type Zr 2 Cl 4 monolayer along x-axis at 300 K is displayed in figures 7(a) and (b), respectively. It is shown that the results obtained by two different methods are generally consistent

Direction
Carrier type well with each other. Then combining the TE parameters obtained before, the dimensionless ZT can be calculated, as depicted in figures 8(a) and (b). In analogy with the PF, the ZT values of Zr 2 Cl 4 also exhibit strong anisotropy. For instance, the ZT value in the x direction peaks at 1.77 with hole doping, which is ten times than that of 0.14 in the y direction at 300 K.
More importantly, the highest ZT value of Zr 2 Cl 4 monolayer reaches 3.90 under the hole concentration of 6.19 × 10 13 cm −2 at 700 K. Such a high ZT value is higher than some other 2D TE materials, such as ZrS 3 (2.44 at 800 K) [62] and isostructural Hf 2 Cl 4 (3.2 at 700 K) [28]. Note also that the highest ZT value of n-type Zr 2 Cl 4 is also as high as 3.60 at  700 K, indicating that doping is effective for enhancing the TE performance of Zr 2 Cl 4 . Given that the highest ZT value can reach 1.77 even at the room temperature, the Zr 2 Cl 4 monolayer has great potential for industrial application in the TE field.

Conclusion
In summary, the TE properties of Zr 2 Cl 4 monolayer are systematically studied by using first-principles calculations combined with semi-classical Boltzmann transport theory. The phonon dispersion and AIMD simulation reveal that the Zr 2 Cl 4 monolayer is thermodynamically stable. The large Seebeck coefficient which results from the band degeneracy and large electrical conductivity lead to the high PF. In addition, the ZT values exhibit significant difference between the x and y directions, owing to the strong anisotropy existing in electrical conductivity and lattice thermal conductivity. Besides, although the great thermal transport property causes unsatisfactory lattice thermal conductivities, an unexpectedly high ZT values of p-type and n-type Zr 2 Cl 4 owing the high PF achieve 3.90 and 3.60 at 700 K, respectively. This study highlights the Zr 2 Cl 4 monolayer as an excellent TE material and its ZT value can be enhanced by doping electronics and holes, which provides a feasible avenue to gain excellent TE materials.

Data availability statement
All data that support the findings of this study are available upon reasonable request to the authors.