Thermoelectric Properties of Hexagonal M2C3 (M = As, Sb, and Bi) Monolayers from First-Principles Calculations

Hexagonal M2C3 compound is a new predicted functional material with desirable band gaps, a large optical absorption coefficient, and ultrahigh carrier mobility, implying its potential applications in photoelectricity and thermoelectric (TE) devices. Based on density-functional theory and Boltzmann transport equation, we systematically research the TE properties of M2C3. Results indicate that the Bi2C3 possesses low phonon group velocity (~2.07 km/s), low optical modes (~2.12 THz), large Grüneisen parameters (~4.46), and short phonon relaxation time. Based on these intrinsic properties, heat transport ability will be immensely restrained and therefore lead to a low thermal conductivity (~4.31 W/mK) for the Bi2C3 at 300 K. A twofold degeneracy is observed at conduction bands along Γ-M direction, which gives a high n-type electrical conductivity. Its low thermal conductivity and high Seebeck coefficient lead to an excellent TE response. The maximum thermoelectric figure of merit (ZT) of n-type can approach 1.41 for Bi2C3. This work shows a perspective for applications of TE and stimulate further experimental synthesis.


Introduction
Thermoelectric (TE) technology can directly convert heat energy into electrical power, playing an important role in solving current energy and environmental crises [1][2][3]. However, low conversion efficiency and high cost are currently facing two crucial bottlenecks [4]. Generally, the conversion efficiency of a TE material is evaluated in terms of a dimensionless thermoelectric figure of merit (ZT) [5][6][7], where S, σ, and T are the Seebeck coefficient, electrical conductivity, and absolute temperature. κ is the thermal conductivity, which is composed of the lattice thermal conductivity κ l and electronic thermal conductivity κ e . However, optimizing one parameter without affecting another is difficult due to the complex competition [8]. The common measures to enhance TE performance mainly concentrate on regulating electrical transport coefficients by band-structure engineering [9] and/or suppressing

Methods
In this paper, the M 2 C 3 is calculated by utilizing the Vienna ab initio simulation package (VASP) [30] within the framework of the Perdew-Burke-Ernzerhof (PBE) [31] generalized gradient approximation [32]. A 9 × 9 × 1 k-mesh and kinetic energy cutoff of 500 eV were used for structure optimization in the Brillouin zone (BZ). A vacuum layer of 20 Å thickness along the z direction was employed to eliminate interlayer interactions. All crystal structures were fully optimized until the total energy variation was less than 10 −6 eV/Å, and the residual forces atoms were less than 0.01 eV/Å. In order to accurately evaluate electronic structure, the Heyd-Scuseria-Ernzerhof (HES06) hybrid density functional was adopted [33]. The electrical transport properties were obtained by semiclassical Boltzmann transport theory as implemented in the BoltzTraP code [34]. This method has successfully predicted many TE materials [35,36]. A dense 45 × 45 × 1 k-mesh was used in the BZ.
Using the ShengBTE code [37], the phonon transport properties were calculated from the Boltzmann transport equation, in which the harmonic second-order interaction force constants (2nd IFCs) and the anharmonic third-order IFCs (3rd IFCs) were used as input. The phonon dispersions and 2nd IFCs were calculated by utilizing PHONONPY packages [38]. A 2 × 2 × 1 supercell with 3 × 3 × 1 k-mesh was used. The anharmonic 3rd IFCs were obtained by using the 2 × 2 × 1 supercell with the finite-difference method [39]. The interactions including the sixth-nearest-neighbor atoms were taken into account for the 3rd IFCs. Based on the test of k-mesh, a dense 35 × 35 × 1 k-mesh was used to calculate k l .

Atomic and Electronic Structures
As illustrated in Figure 1a,b the monolayer M 2 C 3 is a hexagonal crystal system with high space group P6/mmm (No. 191). The optimized lattice constant is 5.86, 6.39, and 6.70 Å for As 2 C 3 , Sb 2 C 3 , and Bi 2 C 3 , which is excellently consistent with a previous theoretical prediction [25]. There are four M (M = As, Sb, and Bi) atoms and six C atoms in the primitive cell. Interestingly, it looks like an enlarging of arsenene [40] with the insertion into a 2D mesh of C atoms from a top view. From the side, the monolayer M 2 C 3 possesses puckered configuration. As shown in Table 1, the interatomic distance of M-C shows a lengthening trend, resulting in the enhancement of atomic vibration frequency. More details are summarized in Table 1. are four M (M = As, Sb, and Bi) atoms and six C atoms in the primitive cell. Interestingly, it looks like an enlarging of arsenene [40] with the insertion into a 2D mesh of C atoms from a top view. From the side, the monolayer M2C3 possesses puckered configuration. As shown in Table 1, the interatomic distance of M-C shows a lengthening trend, resulting in the enhancement of atomic vibration frequency. More details are summarized in Table 1.
The band structures and corresponding projected density of states (PDOS) of monolayer M2C3 are shown in Figure 1. Obviously, As2C3 and Sb2C3 are indirect band gap semiconductors with the valence band maximum (VBM) and the conduction band minimum (CBM) located at the Γ and K points, respectively. Monolayer Bi2C3 exhibits a direct band gap of 0.81 eV, which is smaller than that of As2C3 (1.42 eV) and Sb2C3 (0.92 eV). Near the Fermi level, one can see that the bands primarily stem from the M-p orbitals. Interestingly, the two lowest conduction bands (CB) display overlap along Γ-M direction. A twofold degeneracy is observed in Sb2C3 and Bi2C3, which gives rise to a high n-type electrical conductivity. The PDOS fully shows that the valence band (VB) and CB are mainly occupied by the p orbitals. Remaining orbitals almost have no contribution around the Fermi level. Meanwhile, we also find that it exhibits stair-like PDOS, which can increase the Seebeck coefficient [41]. Therefore, an intrinsic CB degeneracy and a stair-like PDOS are obtained in monolayer M2C3, which are considered to be the electronic transport characteristics of high performance TE devices [42].    The band structures and corresponding projected density of states (PDOS) of monolayer M 2 C 3 are shown in Figure 1. Obviously, As 2 C 3 and Sb 2 C 3 are indirect band gap semiconductors with the valence band maximum (VBM) and the conduction band minimum (CBM) located at the Γ and K points, respectively. Monolayer Bi 2 C 3 exhibits a direct band gap of 0.81 eV, which is smaller than that of As 2 C 3 (1.42 eV) and Sb 2 C 3 (0.92 eV). Near the Fermi level, one can see that the bands primarily stem from the M-p orbitals. Interestingly, the two lowest conduction bands (CB) display overlap along Γ-M direction. A twofold degeneracy is observed in Sb 2 C 3 and Bi 2 C 3 , which gives rise to a high n-type electrical conductivity. The PDOS fully shows that the valence band (VB) and CB are mainly Nanomaterials 2019, 9, 597 4 of 10 occupied by the p orbitals. Remaining orbitals almost have no contribution around the Fermi level. Meanwhile, we also find that it exhibits stair-like PDOS, which can increase the Seebeck coefficient [41]. Therefore, an intrinsic CB degeneracy and a stair-like PDOS are obtained in monolayer M 2 C 3 , which are considered to be the electronic transport characteristics of high performance TE devices [42].

Electrical Transport Properties
The S, σ, and κ e are indispensable for analyzing the TE performance of monolayer M 2 C 3 . The electrical transport properties are obtained by solving semiclassic Boltzmann transport equation [33,43] combing with a constant relaxation time approximation. Here, we further imitate the doping effects of electronic transport by using the rigid band approximation. Based on these methods, the shape of electronic band structure is assumed to be invariant under light doping, and only moves up or down at the Fermi level for n-and p-type doping, respectively [44,45]. The negative and positive µ equate to n-and p-type, respectively. Based on Boltzmann transport equation, the electrical transport coefficients as functions of chemical potential µ and temperature can be derived [33], where αβ and V are Cartesian index and the volume of the primitive cell, and the electrical transport distribution function αβ ( ) is given by where N 0 , i, τ, and ν are the sum of q points, the band index, the electron relaxation time, and the electron group velocity. Figure 2a-c shows the Seebeck coefficients for monolayer M 2 C 3 . Obviously, the temperaturedependent decreasing behavior of the Seebeck coefficients is slowing down along with increasing the temperature. Surprisingly, the monolayer As 2 C 3 has a very high Seebeck coefficient of 2.27 mV/K, which is visibly higher than that of Sb 2 C 3 (1.37 mV/K) and Bi 2 C 3 (1.31 mV/K) at room temperature. Compared with some high-performance TE materials (PbTe [19], SnSe [23], Bi 2 O 2 Se [46]), the monolayer M 2 C 3 exhibits an ultrahigh Seebeck coefficient. These high Seebeck coefficients mainly originate from energy-dependent PDOS as shown in Figure 1. For a doped semiconductor, the Seebeck coefficient can be given by [47], where k B is the Boltzmann constant. Equation (5) implies that the stair-like PDOS contains several sharp peaks, which can enhance carrier concentration n ( ) and give a high Seebeck coefficient. The electrical conductivity with respect to scattering time σ/τ is presented in Figure 2d-f. Unlike the Seebeck coefficient, the maximum σ/τ for As 2 C 3 is lower than that of Sb 2 C 3 and Bi 2 C 3 . Meanwhile, it is found that no matter what kind of materials, the σ/τ of n-type is always higher than that of the p-type, which is mainly attributed to the PDOS. We can see clearly that the slope of σ/τ will be flattened in low µ region (−0.5-0.5 eV) with the temperature increasing. The κ e with respect to the relaxation time is shown in Figure 2g-i via the Wiedemann-Franz law [48] κ e = LσT where L = π 2 κ B 2 /3e 2 is the Lorenz number. Similar to the σ/τ, the κ e /τ displays analogous curves.
The κ e /τ gradually increases along with varying the absolute value of the µ from the Fermi energy level (µ = 0). These transport coefficients indicate that the monolayer M 2 C 3 has a good performance of electronic transport ability with n-type.
Nanomaterials 2019, 9, x FOR PEER REVIEW 5 of 11 temperature increasing. The κe with respect to the relaxation time is shown in Figure 2g-i via the Wiedemann-Franz law [48] = (6) where L = π 2 κB 2 /3e 2 is the Lorenz number. Similar to the σ/τ, the κe/τ displays analogous curves. The κe/τ gradually increases along with varying the absolute value of the μ from the Fermi energy level (μ = 0). These transport coefficients indicate that the monolayer M2C3 has a good performance of electronic transport ability with n-type.

Phonon Transport Properties
The phonon scattering curves of monolayer M2C3 is shown in Figure 3. No imaginary phonon frequencies appear in the phonon spectra, which indicates that the monolayer M2C3 is thermodynamically stable at ambient pressure. The acoustic branches of M2C3 exhibit a common phenomenon in 2D systems with a parabolic dispersion of out-of-plane acoustic mode and two linear dispersions of in-plane modes at the Γ point [49]. As shown in Figure 3, one can see that the M atoms dominate the low frequency region (below~10 THz), while the remaining area is from the C atomic contributions. Meanwhile, we also find that the acoustic modes for As2C3, Sb2C3, and Bi2C3 exhibit a downward moving trend, which can be attributed to the larger of atomic mass. It is noted that the low frequency optical modes are alternating and softening with the three acoustic branches for the monolayer M2C3, resulting in strong acoustic-optical interactions. This is similar to the PbSe [50], which has strong anharmonic effects. The boundary frequency of lowest optical branch displays a decreasing trend with the following order: As2C3 (4.18 THz) > Sb2C3 (3.08 THz) > Bi2C3 (2.12 THz). To further analyze phonon scattering properties, the corresponding phonon density of states (PhDOS) is presented. From the PhDOS, we can see clearly that the acoustic phonon branches mainly contain the M vibrations, while the contributions from C (xy) vibrations are mainly limited to 10 to 25 THz. The high frequency region from 45 to 50 THz is fully occupied by C (z) vibrations. In addition, the PhDOS also shows several peaks

Phonon Transport Properties
The phonon scattering curves of monolayer M 2 C 3 is shown in Figure 3. No imaginary phonon frequencies appear in the phonon spectra, which indicates that the monolayer M 2 C 3 is thermodynamically stable at ambient pressure. The acoustic branches of M 2 C 3 exhibit a common phenomenon in 2D systems with a parabolic dispersion of out-of-plane acoustic mode and two linear dispersions of in-plane modes at the Γ point [49]. As shown in Figure 3, one can see that the M atoms dominate the low frequency region (below~10 THz), while the remaining area is from the C atomic contributions. Meanwhile, we also find that the acoustic modes for As 2 C 3 , Sb 2 C 3 , and Bi 2 C 3 exhibit a downward moving trend, which can be attributed to the larger of atomic mass. It is noted that the low frequency optical modes are alternating and softening with the three acoustic branches for the monolayer M 2 C 3 , resulting in strong acoustic-optical interactions. This is similar to the PbSe [50], which has strong anharmonic effects. The boundary frequency of lowest optical branch displays a decreasing trend with the following order: As 2 C 3 (4.18 THz) > Sb 2 C 3 (3.08 THz) > Bi 2 C 3 (2.12 THz). To further analyze phonon scattering properties, the corresponding phonon density of states (PhDOS) is presented. From the PhDOS, we can see clearly that the acoustic phonon branches mainly contain the M vibrations, while the contributions from C (xy) vibrations are mainly limited to 10 to 25 THz. The high frequency region from 45 to 50 THz is fully occupied by C (z) vibrations. In addition, the PhDOS also shows several peaks especially in optical branches region, which can give rise to the small phonon group velocity. Nanomaterials 2019, 9, x FOR PEER REVIEW 6 of 11 especially in optical branches region, which can give rise to the small phonon group velocity. Based upon the phonon kinetic theory, the κl can be calculated as below: where CV, ν α is the phonon specific heat and the phonon group velocity along α direction, and τqλ is the phonon relaxation time for the phonon mode λ at the wave vector q. The κl of monolayer M2C3, at temperatures from 300 K to 800 K, are plotted in Figure 4a. Obviously, the intrinsic κl shows evident temperature dependence, proportional to the inverse of the temperature 1/T. In the high temperature, this dependence behavior is deemed to a common phenomenon of crystals, which attribute to the intrinsic enhancement of phonon-phonon scattering. At room temperature, we can see clearly that the Bi2C3 exhibits an intrinsic low κl of 4.31 W/mK, which is lower than those of the Sb2C3 (9.53 W/mK) and As2C3 (20.82 W/mK). In general, the κl is mainly dominated by the acoustical modes for the monolayer systems [51]. Because the mass of Bi is larger than that of As and Sb, acoustical phonons have lower frequencies for Bi2C3 than for As2C3 and Sb2C3. These acoustic modes with low frequency might cause the small phonon group velocities, leading to the low κl.
The phonon group velocities are the important indicator for the assessment of heat transport ability. Using the phonon dispersion, the phonon group velocities can be obtained by The corresponding group velocities are plotted in Figure 4b. Much large values of the group velocities can be observed in low frequency region, while the high frequency modes exhibit relatively small group velocities. Meanwhile, it can be seen clearly that the group velocities from large to small is monolayer As2C3 > Sb2C3 > Bi2C3, which is consistent with the above analysis. This phenomenon similar to the GeS and SnS [52], which may attribute to large atomic mass.
To identify the underlying mechanism of low κl, we introduce the Grüneisen parameters and phonon relaxation time as shown in Figure 4c,d. The Grüneisen parameter can fully reflect the anharmonic interactions of a crystal, which is essential for analyzing the intrinsic characteristics of κl. The Grüneisen parameter can be described by Based upon the phonon kinetic theory, the κ l can be calculated as below: where C V , ν α is the phonon specific heat and the phonon group velocity along α direction, and τ qλ is the phonon relaxation time for the phonon mode λ at the wave vector q. The κ l of monolayer M 2 C 3 , at temperatures from 300 K to 800 K, are plotted in Figure 4a. Obviously, the intrinsic κ l shows evident temperature dependence, proportional to the inverse of the temperature 1/T. In the high temperature, this dependence behavior is deemed to a common phenomenon of crystals, which attribute to the intrinsic enhancement of phonon-phonon scattering. At room temperature, we can see clearly that the Bi 2 C 3 exhibits an intrinsic low κ l of 4.31 W/mK, which is lower than those of the Sb 2 C 3 (9.53 W/mK) and As 2 C 3 (20.82 W/mK). In general, the κ l is mainly dominated by the acoustical modes for the monolayer systems [51]. Because the mass of Bi is larger than that of As and Sb, acoustical phonons have lower frequencies for Bi 2 C 3 than for As 2 C 3 and Sb 2 C 3 . These acoustic modes with low frequency might cause the small phonon group velocities, leading to the low κ l . The phonon group velocities are the important indicator for the assessment of heat transport ability. Using the phonon dispersion, the phonon group velocities can be obtained by The corresponding group velocities are plotted in Figure 4b. Much large values of the group velocities can be observed in low frequency region, while the high frequency modes exhibit relatively small group velocities. Meanwhile, it can be seen clearly that the group velocities from large to small is monolayer As 2 C 3 > Sb 2 C 3 > Bi 2 C 3 , which is consistent with the above analysis. This phenomenon similar to the GeS and SnS [52], which may attribute to large atomic mass.
To identify the underlying mechanism of low κ l , we introduce the Grüneisen parameters and phonon relaxation time as shown in Figure 4c,d. The Grüneisen parameter can fully reflect the anharmonic interactions of a crystal, which is essential for analyzing the intrinsic characteristics of κ l . The Grüneisen parameter can be described by It is noted that the large Grüneisen parameters can be observed at a low frequency region. Usually, large Grüneisen parameters (absolute value) indicate strong anharmonicity, which can give rise to low Nanomaterials 2019, 9, 597 7 of 10 κ l . The averages of the acoustic Grüneisen parameters (absolute value) are calculated to be 2.61, 4.25, and 4.46 for As 2 C 3 , Sb 2 C 3 , and Bi 2 C 3 , respectively. Obviously, the Bi 2 C 3 exhibits larger Grüneisen parameter than that of As 2 C 3 and Sb 2 C 3 , indicating that Bi 2 C 3 possesses stronger anharmonicity. To further explore the thermal transport properties, we show the phonon relaxation time in Figure 4d. The phonon relaxation time of Bi 2 C 3 is much shorter than that of the As 2 C 3 and Sb 2 C 3 at 300 K, which is a significant factor for the low κ l of Bi 2 C 3 . More phonon transport details are summarized in Table 2. . (9) It is noted that the large Grüneisen parameters can be observed at a low frequency region. Usually, large Grüneisen parameters (absolute value) indicate strong anharmonicity, which can give rise to low κl. The averages of the acoustic Grüneisen parameters (absolute value) are calculated to be 2.61, 4.25, and 4.46 for As2C3, Sb2C3, and Bi2C3, respectively. Obviously, the Bi2C3 exhibits larger Grüneisen parameter than that of As2C3 and Sb2C3, indicating that Bi2C3 possesses stronger anharmonicity. To further explore the thermal transport properties, we show the phonon relaxation time in Figure 4d. The phonon relaxation time of Bi2C3 is much shorter than that of the As2C3 and Sb2C3 at 300 K, which is a significant factor for the low κl of Bi2C3. More phonon transport details are summarized in Table 2.

Thermoelectric Figure of Merit
Based on these transport coefficients, we estimate the ZT of M2C3 at three typical temperatures. The electronic relaxation time τ is employed from a previous report [25] and  Table 2. Summary of the lattice thermal conductivity κ l (W/mK), the averages of the acoustic group velocity v (km/s) and Grüneisen parameters γ, and the lowest optical frequency ω o (THz) for the As 2 C 3 , Sb 2 C 3 , and Bi 2 C 3 , respectively.

Thermoelectric Figure of Merit
Based on these transport coefficients, we estimate the ZT of M 2 C 3 at three typical temperatures. The electronic relaxation time τ is employed from a previous report [25] and calculated by deformation potential theory [53]. The calculated ZT of M 2 C 3 is presented in Figure 5. We can see that the two peaks move gradually towards the Fermi energy level (µ = 0 eV), which originate from the decrease of band gap. Unlike the Seebeck coefficient, the ZT value has an enhancing behavior with the increasing temperature. Meanwhile, we can note that the ZT of n-type is always higher than that of the p-type for M 2 C 3 , which mainly stems from their difference in electrical conductivity. In general, if the ZT value approach 1, it can be considered as a good TE material [54]. Due to the excellent transport performance of electrons, the maximum ZT value can reach 0.93, 1.17, and 1.41 for As 2 C 3 , Sb 2 C 3 , and Bi 2 C 3 , respectively. Nanomaterials 2019, 9, x FOR PEER REVIEW 8 of 11 calculated by deformation potential theory [53]. The calculated ZT of M2C3 is presented in Figure 5. We can see that the two peaks move gradually towards the Fermi energy level (μ = 0 eV), which originate from the decrease of band gap. Unlike the Seebeck coefficient, the ZT value has an enhancing behavior with the increasing temperature. Meanwhile, we can note that the ZT of n-type is always higher than that of the p-type for M2C3, which mainly stems from their difference in electrical conductivity. In general, if the ZT value approach 1, it can be considered as a good TE material [54]. Due to the excellent transport performance of electrons, the maximum ZT value can reach 0.93, 1.17, and 1.41 for As2C3, Sb2C3, and Bi2C3, respectively.

Conclusions
In summary, we report on the TE properties of monolayer M 2 C 3 , through using the first-principles calculations and solving Boltzmann transport equation. The results indicate the monolayer M 2 C 3 exhibits low κ l at room temperature, especially for the Bi 2 C 3 (~4.31 W/mK). Compared to most TE materials, the Bi 2 C 3 exhibits analogous intrinsic properties, such as small group velocity, short relaxation time, and large Grüneisen parameters. A twofold degeneracy and stair-like PDOS are observed, which can lead to a high Seebeck coefficient. Finally, combining with available transport parameters, the ZT is found to be about 0.93, 1.17, and 1.41 for the As 2 C 3 , Sb 2 C 3 , and Bi 2 C 3 at 700 K, respectively. This work indicates that the M 2 C 3 is very promising for TE applications.