Ultra-low thermal conductivity induced by significant anharmonic phonon scattering in two-dimensional graphitic metal (Sb and Bi) carbides

Abstract

Two-dimensional (2D) thermoelectric (TE) materials have attracted much attention in recent years. One major factor that limits the TE performance of 2D materials is their relatively high thermal conductivity. Searching for 2D materials with low thermal conductivity is therefore one of the central goals in theoretical and computational studies of 2D TE materials. Here, we report a comprehensive first-principles study of the TE properties of 2D graphitic metal carbides (g-MCs), recently proposed highly stable 2D crystals made of carbon and metal. After theoretically examining the entire family of 2D g-MCs, we found that 2D g-MCs (M = Sb, Bi, Al, Ga, Rh, Ir, and W) possess intrinsically low thermal conductivity and the two members among them, 2D g-SbC and g-BiC, exhibit ultra-low lattice thermal conductivities of 0.78 and 0.90 W m⁻¹ K⁻¹, respectively, at room temperature. Detailed analysis shows that the underlying ultra-low thermal conductivities of 2D g-SbC and g-BiC are the significant anharmonic effects in acoustic phonon modes originating from the extended metal–carbon bonds in their unique structure, which lead to strong anharmonic phonon scattering and highly reduced phonon group velocities. With a combination of low thermal conductivity and favorable electron conductivity, p-type Sb₂C₁₂ shows high TE figures-of-merit of around 0.98 at 300 K and 2.55 at 900 K. These findings offer new opportunities for future applications of 2D TE materials.

---

1 Introduction

Thermoelectric (TE) materials have the capability of generating electricity directly from the temperature gradient, holding significant promise for sustainable energy applications. In principle, the high performance of a TE material requires a high power factor (PF = S²σ, where S and σ are the electrical conductivity and the Seebeck coefficient, respectively) and low thermal conductivity (κ).^1–5 Due to the interdependence between σ, S, and κ, it is challenging to improve these parameters simultaneously in a TE material.^6,7 Two-dimensional (2D) materials have been shown to have great potential for TE applications because of their unique structures that often lead to high electron conductivity due to the lack of impurity scattering in the direction perpendicular to the 2D plane. Various 2D materials, including group IVA/VA–VIA heterostructures, transition metal dichalcogenides (TMDCs), 2D carbon allotropes and their derivatives (graphene and graphyne families), have been intensively studied for TE applications in experiments and theoretical predictions.^8–20 One of the main factors that limits the TE efficiency of 2D materials has been reported to be their high thermal conductivity.^21–23 As a reference, MoS₂, which has been widely studied as a 2D TE material, possesses a thermal conductivity around of 34.5 W m⁻¹ K⁻¹.²¹ Finding 2D materials with ultra-low thermal conductivity is therefore essential for future applications of 2D TE materials.

2D graphitic metal carbides (g-MCs) are a recently proposed new family of atomically thin crystals with metal-C₃ (MC₃) moieties periodically distributed in a graphenic lattice.^24,25 Although 2D g-MCs have not yet been synthesized experimentally, several recent advances naturally suggest a plausible realization route in which metal precursors are loaded and trapped, rapidly reduced under confinement, and then driven to undergo surface lattice reconstruction on γ-graphyne templates.^26–30 2D g-MCs have a chemical formula of M₂C₁₂. Because of the relatively strong long-distance metal–metal interactions originating from their unique structures, 2D g-MCs are extremely stable with highly tunable electronic and magnetic properties.²⁵ Consequently, 2D g-MCs have been considered in various applications including catalysis and electronic/spintronic devices.^{24,25,31–35} The materials’ TE properties, however, have not been discussed so far. Using first-principles methods based on density functional theory (DFT) and the Boltzmann transport approach, we theoretically investigated the TE performance of this intriguing family of 2D materials. We found that the highly tunable electronic properties of 2D g-MCs offer unique opportunities for greatly lowering the thermal conductivity and significantly improving the figure-of-merit, the ZT value (ZT = S²σT/κ) of the materials. In particular, all semiconducting 2D g-MCs are found to have low intrinsic thermal conductivities, with two members of the family, 2D g-SbC and g-BiC, exhibiting ultra-low lattice thermal conductivities of 0.78 and 0.90 W m⁻¹ K⁻¹, respectively, at room temperature. Detailed analysis reveals that anharmonic three-phonon and four-phonon scattering play a key role in the calculated low thermal conductivity of 2D g-SbC and g-BiC. Together with its favorable electron conductivity, p-type 2D g-SbC is identified as the best TE material among all g-MCs with ZT values of around 0.98 at 300 K and 2.55 at 900 K.

2 Computational methods

We performed first-principles calculations based on DFT using the Quantum Espresso (QE) package.³⁶ Standard solid-state pseudopotentials³⁷ were employed to model the ion core potentials and valence electrons. The Perdew–Burke–Ernzerhof functional³⁸ within the generalized gradient approximation was used to describe the exchange–correlation energy. For structural relaxation and electronic structure calculations, the Kohn–Sham orbitals were expanded in a plane-wave basis with kinetic energy cutoffs of 80 Ry for wavefunctions and 320 Ry for charge density and the potential. A 12 × 12 × 1 Monkhorst–Pack grid was used to sample the k-points across the entire Brillouin zone. Convergence criteria for the total energy and force were set to 10⁻¹² Ry and 10⁻⁶ Ry per Bohr, respectively. Phonon dispersion, dynamical matrices, and harmonic interatomic force constants (IFCs) were calculated using density functional perturbation theory,³⁹ as implemented in the QE package.

We adopted the Boltzmann transport theory with the relaxation time approximation (RTA) to predict the lattice thermal conductivity of 2D g-MCs. Both three-phonon and four-phonon interactions were considered as the primary scattering processes in the linearized Boltzmann transport equation (BTE). The corresponding third-order and fourth-order IFCs were calculated using the finite displacement approach. A 3 × 3 × 1 supercell was used for the finite difference method with the THIRDORDER⁴⁰ and FOURTHORDER packages,⁴¹ incorporating interactions up to the sixth-nearest neighbors for the third-order IFCs and the third-nearest neighbors for the fourth-order IFCs. In the RTA framework, the phonon relaxation time τ_λ for each phonon mode λ was calculated as follows:

andwhere Γ^+/−_{λλ′λ′′} denotes the three-phonon scattering rate for the combination/splitting process and Γ^{++/+−/−−}_{λλ′λ′′λ′′′} represents the four-phonon scattering rate for the combination/redistribution/splitting process. Explicit expressions for Γ are provided in the SI. N is the total number of q points used to sample the Brillouin zone, where a 15 × 15 × 1 grid was used for both three-phonon and four-phonon processes. To solve the linearized BTE, we employed the ShengBTE package,⁴² which offers both iterative solutions to the BTE and results in the RTA. In our calculations, three-phonon and four-phonon processes were all addressed using the iterative method. The lattice thermal conductivity tensor κ^αβ_L contributed by phonons was calculated using the following expression:where N is the total number of modes in the Brillouin zone and V₀ represents the volume of the unit cell. F_λ denotes the linear deviation of the non-equilibrium distribution function from the equilibrium Bose–Einstein distribution function f⁰_λ = (e^ħωλ/k_BT − 1)⁻¹. Within the RTA, F_λ was calculated simply as the phonon mean free path F_λ = τ_λv_λ; however, for the iterative method, a correction term Δ_λ representing the error introduced by the RTA was added to the phonon group velocity, such that F_λ = τ_λ(v_λ + Δ_λ).

To predict electronic transport driven by a temperature gradient and evaluate the associated thermoelectric properties, we calculated the electron–phonon interaction (EPI) to accurately determine the electron relaxation time (τ_e). Unlike traditional deformation potential theory,^43,44 EPI calculations yield τ^n,k_e as a function of band index n and wave vector k, rather than as a global constant. The EPW code⁴⁵ was used to compute the electron–phonon matrix elements g_mn,ν(k, q) and the electron self-energy Σ_nk(ω, T) with interpolation on a 72 × 72 × 1 dense k-point and a q-point mesh. The electron relaxation time τ^n,k_e is given by the imaginary part of the electron self-energy:

With τ^n,k_e, we calculated the electrical conductivity, Seebeck coefficient, and electronic thermal conductivity by solving the BTE using the BoltzTraP2 package:⁴⁶andwhere

3 Results and discussion

3.1 Crystal structure and stability

After the first proposal of monolayer 2D g-MCs, Yam et al.²⁵ conducted a comprehensive study on the ground-state structures of g-MCs and demonstrated that 33 types of 2D g-MCs with different metal elements are thermally stable at room temperature. In this work, we focused on the 7 types of stable monolayer 2D g-MCs that exhibit semiconducting behavior (M = Sb, Bi, Al, Ga, Rh, Ir, and W). As illustrated in Fig. 1(a), the ground-state structure of selected g-MCs is a hexagonal lattice with P3̄ group symmetry in the unit cell. In the main text, we focus on the properties of monolayer Sb₂C₁₂ due to its superior TE performance among monolayer 2D g-MCs. Detailed analysis of the other 6 types of 2D g-MCs are shown in Fig. S1–S9.

Fig. 1 Top and side views of (a) monolayer Sb₂C₁₂ and (b) the structure during the AIMD process at 900 K and 12 ps, with the structural unit indicated by the black dashed line. (c) Phonon dispersion of monolayer Sb₂C₁₂. (d) Projected phonon density of states. (e) Variation in the Sb–C bond length (red line) and free energy (blue line) of Sb₂C₁₂ during the AIMD simulation at 900 K.

The optimized lattice constant of monolayer Sb₂C₁₂ is 6.74 Å, and the average Sb–C bond length is 2.18 Å. To verify the dynamic stability of the crystal, we calculated the phonon dispersion of monolayer Sb₂C₁₂. As shown in Fig. 1(c), no imaginary phonon modes are observed, indicating the dynamic stability of the ground-state structure. Notably, a clear phonon band gap exists between 85 cm⁻¹ and 140 cm⁻¹ (A–O band gap), separating the acoustic and optical phonon modes. The element-decomposed phonon density of states (PDOS) for monolayer Sb₂C₁₂ is shown in Fig. 1(d), where collective vibrations of Sb atoms dominate the acoustic and low-frequency optical regions, indicating the importance of metal atoms in phonon transport and scattering processes. Furthermore, ab initio molecular dynamics (AIMD) calculations were conducted to assess the thermal stability of monolayer Sb₂C₁₂ at elevated temperatures. As shown in Fig. 1(b) and (e), the average Sb–C bond length oscillates around its static value of 2.18 Å throughout a 6000-step MD simulation with a time step of 2 fs. Besides, the regular variation in free energy further confirms the thermal stability of monolayer Sb₂C₁₂ at 900 K.

3.2 Phonon transport

An accurate lattice thermal conductivity is crucial to evaluate the TE performance of 2D g-MCs. For materials with strongly anharmonic phonon interactions, high-order anharmonicities should be considered to obtain a reliable κ_L and its corresponding temperature dependence. The harmonic approximation (HA)^47,48 and self-consistent phonon (SCP) theory^49–51 are two approaches for quantifying anharmonic effects on phonon transport. Based on the quasiparticle (QP) approximation, the state-of-the-art SCP theory effectively captures the impact of anharmonic lattice dynamics on phonon dispersion at high temperatures.^52,53 However, its predictions for κ_L still require further validation against the corresponding experimental measurements. In contrast, the HA, which simply employs 0 K harmonic phonon dispersion, neglects temperature-induced renormalization of phonon energies. Nevertheless, HA has been shown to be a reasonable approximation for materials where high-order anharmonicity cannot be overlooked.^54–56 For instance, Ruan et al.⁵⁵ demonstrated that HA, when considering 3-phonon (3-ph) and 4-phonon (4-ph) scattering, slightly overestimates the κ_L value of III–V semiconductors up to 900 K in contrast to the experimental results. In our study, monolayer Sb₂C₁₂ is a semiconductor with a direct band gap of 1.17 eV and contains the group VA element (Sb). Thus, we used the HA as an appropriate approximation to calculate the κ_L value of monolayer Sb₂C₁₂.

As shown in Fig. 2(a), the lattice thermal conductivity κ_3ph (only with 3-ph scattering) and κ_3ph,4ph (including both 3-ph and 4-ph scattering) along the lattice a-axis was calculated over a temperature range from 200 K to 900 K. At 300 K, the κ_3ph and κ_3ph,4ph values reach 1.23 W m⁻¹ K⁻¹ and 0.78 W m⁻¹ K⁻¹, indicating a 36% reduction in thermal conductivity due to 4-ph scattering. Additionally, as the temperature increases, the impact of quartic anharmonicity becomes more pronounced: at 600 K and 900 K, the κ_3ph,4ph values are 0.38 and 0.23 W m⁻¹ K⁻¹, which are 48% and 56% lower than those of κ_3ph, respectively. In contrast to other 2D materials, the lattice thermal conductivity of monolayer Sb₂C₁₂ is exceptionally low—around 0.15% of the value for silicon-supported monolayer graphene (κ_L ≈ 500 W m⁻¹ K⁻¹).⁵⁷ To further investigate the importance of quartic anharmonicity in 2D g-MCs, we fitted the temperature dependence of lattice thermal conductivity with T^−α. In general, the intrinsic κ_L of semiconductors follows Debye's T³ law in the low temperature limit, but as the temperature is raised, it typically follows a T⁻¹ relationship due to strong Umklapp scattering.⁵⁸ As shown in Fig. 2(a), when only cubic anharmonicity is considered, κ_3ph exhibits an anomalous temperature dependence with an α value of 0.72, deviating from the universal T⁻¹ law. However, the inclusion of 4-ph scattering improves the dependence and increases the α value to 1.03, suggesting inevitable quartic anharmonicity in 2D g-MCs.

Fig. 2 (a) Calculated temperature-dependent κ_L values of Sb₂C₁₂ from 200 K to 900 K, with consideration of only 3-phonon scattering and 3-phonon plus 4-phonon scattering. The dashed lines represent the curve fitting of κ_L with T^−α. Frequency-dependent anharmonic scattering rates (ASRs) of Sb₂C₁₂ (b) at 300 K and (d) at 900 K, contributed by 3-phonon and 4-phonon processes, respectively. The purple solid line indicates the intrinsic vibrational frequency of phonons. (c) Frequency cumulative κ_L and its 1^st derivative (inset) at 300 K. The shaded region between the gray dashed lines represents the A–O band gap of Sb₂C₁₂.

A key reason for the ultra-low κ_L value of Sb₂C₁₂ is its high phonon scattering rates (SRs) in the 3-ph and 4-ph scattering processes. As shown in Fig. 2(b) and (d), the 4-ph SRs are large enough that they cannot be neglected in the total SRs, with their proportion increasing as the temperature increases. At 300 K, the 4-ph SRs are obviously lower than the 3-ph SRs across the entire frequency range, but in the acoustic region, the averaged 4-ph SR is almost in the same order of magnitude as the averaged 3-ph SR, with a ratio of around 10.5. In contrast, when the temperature increases, the 4-ph SRs are significantly enhanced and the same ratio decreases to around 3.5 at 900 K, demonstrating the effect of 4-ph scattering on the phonon relaxation process. Additionally, to verify the validity of the quasiparticle (QP) approximation, we compared the intrinsic phonon vibrational frequencies with phonon SRs for both 3-ph and 4-ph processes. As shown in Fig. 2(b), the dominant SRs fall below the purple solid line, indicating that phonon vibrational periods are much longer than the phonon relaxation times. This ensures that phonons undergo several full oscillations before scattering events occur, which is required in the QP approximation. To illustrate the frequency-dependent behavior of 3-ph and 4-ph SRs, it is essential to consider the influence of selection rules on these scattering processes.^59–61 Due to the selection rules, the total number of allowed 3-ph patterns in the phase space is considerably lower than that of the 4-ph processes, with a ratio of around 0.7% in the Sb₂C₁₂ case. Among the allowed 3-ph processes, as the frequency range of the acoustic phonon (Δω_A ∼85cm⁻¹) is slightly larger than the A–O phonon band gap (Δω_A−O ∼55 cm⁻¹), acoustic–acoustic–optical (AAO) 3-ph patterns are greatly reduced in the low-frequency region, resulting in a marked decrease in the 3-ph SRs below 55 cm⁻¹. However, the restrictions imposed by the selection rules on the 4-ph processes are much fewer, so the drop of 4-ph SRs is not as violated as that of 3-ph SRs below 55 cm⁻¹. Consequently, within the acoustic region, the 4-ph SRs are comparable to the 3-ph SRs. Furthermore, as shown in Fig. 2(c), the frequency cumulative κ_L and its derivative indicate that the primary contribution to the total κ_L arises from acoustic phonons, where 4-ph scattering largely suppresses the increase of κ_L. From the analysis of the phonon PDOS, we determined that the low-frequency phonons mainly originate from the vibrations of metal atoms, so the ‘rattling’ of metal atoms in the cage of surrounding carbons becomes a direct incentive for the strong anharmonicity.

To provide a deeper understanding of the mechanisms underlying the 3-ph and 4-ph scattering, we analyzed the decomposed phonon SRs for each type of process. Fig. 3 comprises the decomposed 3-ph and 4-ph SRs at 300 K. As illustrated in Fig. 3(a), the 3-ph SRs are mainly dominated by the splitting process at high frequencies above 800 cm⁻¹ and the combination process at low frequencies below 85 cm⁻¹, with a relatively balanced contribution from both processes in the intermediate frequency range. Due to the presence of the A–O band gap, the SRs from the combination process for low-energy acoustic phonons remain relatively low. This limitation arises because the excitation of an acoustic phonon into the optical region requires at least one participant phonon with a frequency exceeding 55 cm⁻¹. In contrast, as shown in Fig. 3(b), the A–O band gap has minimal impact on the 4-ph SRs in the predominant redistribution process. The AOAO pairing mechanism is not restricted by the band gap, with the SR values derived from the redistribution process being almost an order of magnitude higher than those of the other processes and comparable to those of the 3-ph SRs in the acoustic region. Fig. 3(c) and (d) further illustrate the SRs associated with normal and Umklapp processes in the 3-ph and 4-ph interactions. As shown in the insets, the normal process is the primary scattering mechanism in 3-ph scattering, while the Umklapp process predominates in 4-ph scattering. This distinction aligns with the higher α factor in the temperature dependence of κ_L ≈ T^−α when incorporating both 3-ph and 4-ph processes. Moreover, we also examined whether electron–phonon scattering may become another potential source of thermal resistance. Unfortunately, as shown in Fig. S10(b), the EPI contribution to the phonon SRs is orders of magnitude lower than that derived from phonon–phonon anharmonic scattering, so we neglected this contribution in our calculations.

Fig. 3 (a) Decomposed 3-phonon scattering rate into the combination (λ + λ′ → λ′′) and splitting (λ → λ′ + λ′′) processes at 300 K. (b) Decomposed 4-phonon scattering rate into the combination (λ + λ′ → λ′′), splitting (λ → λ′ + λ′′ + λ′′′) and redistribution (λ + λ′ → λ′′ + λ′′′) processes at 300 K. (c) Decomposed 3-ph and (d) 4-ph scattering rate into the normal and Umklapp processes at 300 K. The inset shows the ASR ratio between the Umklapp and normal processes.

We also considered other key factors to illustrate the ultra-low κ_L value of monolayer Sb₂C₁₂. As shown in Fig. 4(a) and (b), we calculated the Grüneisen parameter (γ) and phonon group velocity (|V_g|). Compared to optical phonons, the significantly higher |γ| values of acoustic phonons indicate a stronger anharmonic effect within the acoustic branches, which emphasizes the critical role of metal atoms in inducing lattice anharmonicity. Besides, in Table S1, we list the γ values of 2D g-MCs and other representative 2D materials. Compared with graphene and MoS₂, 2D g-MCs exhibit higher |[small gamma, Greek, macron]_a| values and much lower acoustic phonon group velocities and relaxation times, leading to their characteristically low thermal conductivity. Meanwhile, the relatively low group velocities of phonons directly contribute to the low κ_L value, as shown in eqn (4). We observed that the low-frequency optical branches are almost dispersion-less and evaluated their impact on κ_L by using the corresponding group velocity |V_g|. The average group velocity of low-frequency optical modes |V̄_low-o| = 0.66 km s⁻¹ is less than one-third of the average acoustic phonon group velocity |V̄_a| = 2.05 km s⁻¹, so its direct contribution to κ_L is much smaller than that of the acoustic branches. As shown in Fig. 2(c), the frequency-cumulative κ_L further confirms its suppressive role with a contribution of less than 0.1 W m⁻¹ K⁻¹ to κ_L. Moreover, the vibrational eigenmodes of the lowest-frequency optical branch are illustrated and discussed in Fig. S11. Besides, the decomposed weighted phase space (WPS) of 3-ph and 4-ph scattering processes is plotted in Fig. 4(c) and (d), which shows the multi-phonon scattering processes allowed by energy conservation. The WPS of allowed 4-ph scattering patterns is comparable to that of 3-ph scattering, indicating that 4-ph scattering significantly increases the total scattering rate in the frequency domain. For both 3-ph and 4-ph scattering, the substantial distribution of WPS in the acoustic region aligns well with the Grüneisen parameter—acoustic phonons predominate in the allowed anharmonic scattering patterns. Thus, with these favorable factors, monolayer Sb₂C₁₂ naturally exhibits a low lattice thermal conductivity. In real experimental settings, some other factors, such as boundary scattering and the substrate effect, may further reduce lattice thermal conductivity, potentially yielding a lower κ_L value than our theoretical predictions. As illustrated in Fig. S12, we evaluated the lattice thermal conductivity of three representative defect configurations in Sb₂C₁₂, and all of them exhibit similarly low κ_L values.

Fig. 4 (a) Grüneisen parameter γ and (b) phonon group velocity |V_g| of Sb₂C₁₂ at 300 K. (c) Decomposed weighted phase space for 3-ph scattering processes and (d) 4-ph scattering processes.

3.3 Electronic structure

Thermoelectric properties, such as the Seebeck coefficient and electrical conductivity, are closely related to the electronic structure of materials. The band structure of monolayer Sb₂C₁₂ is depicted in Fig. 5(a), where the valence band maximum (VBM) and the conduction band minimum (CBM) are both located at the M point in k-space, forming a direct band gap of approximately 1.17 eV. Besides, at a point between Γ and K, another pair of direct VB and CB edges lies near the VBM and CBM, with a difference of around 0.2 eV on both sides. Extensive research has shown that multi-valley effects can significantly enhance the Seebeck effect.^62–64 Therefore, the presence of multiple valleys in the band structure provides a promising direction for further improving the thermoelectric performance of monolayer Sb₂C₁₂. Techniques such as doping and strain engineering, aimed at reducing the energy difference between the VB and CB edges, may further enhance the Seebeck coefficient of this 2D g-MC. Furthermore, as plotted in Fig. 5(b), the electron DOS possesses a steep slope at the VBM and CBM, which is advantageous for achieving a high Seebeck coefficient. The element-projected DOS reveals that both the VBM and CBM are mainly contributed by light C atoms, and this is consistent with the absence of notable spin–orbital coupling splitting near the Fermi level. To visualize the electronic structure around the band edges, we present the band-decomposed charge densities of the VBM and CBM in Fig. 5(c) and (f). With a greater contribution of metal atoms in the PDOS, the VBM charge density is localized more around the metal atoms, whereas the CBM charge density is more extensive on the carbon skeleton. Consequently, metal–carbon bonds easily form bridges in hole transport for p-type Sb₂C₁₂, while electronic transport is more efficient through the carbon skeleton in the case of n-type doping. The electron localization function (ELF) is also shown in Fig. 5(e) and its value at the midpoint between Sb and C is around 0.77, indicating that the covalent bond between the metal atom and carbon is not so strong as to hinder the rattling of metal atoms and the corresponding lattice anharmonicity. Additionally, the effective masses of holes and electrons at the VBM and CBM are 1.07m₀ and 1.01m₀, where m₀ denotes the free electron mass. The light carrier effective masses are also beneficial for achieving high electrical conductivity and Seebeck coefficients.

Fig. 5 (a) Element-decomposed band structure and (b) projected density of states of Sb₂C₁₂, with the color bar indicating the weight of Sb in the band structure. The Fermi level was set to zero. Decomposed charge density of the (c) valence band maximum and (f) conduction band minimum at the high-symmetry point, M. (d) Momentum-dependent electron relaxation times τ^n,k_e of Sb₂C₁₂ at 300 K, 600 K and 900 K. (e) 2D ELF at a plane perpendicular to the z-axis between the carbon skeleton and the protruding Sb atom in the positive z direction.

Another vital factor to determine the TE performance is the electron relaxation time τ_e. Based on the electronic and phonon structure, we calculated the electron–phonon matrix and the electron self-energy in the EPI. To interpolate the electron–phonon matrix on a denser k-grid, maximum localized Wannier functions (MLWFs) are employed as an intermediary for transformation. As shown in Fig. S10(a), MLWFs accurately reproduce the band structure in close agreement with the plane-wave basis, validating the reliability of Wannierization in our calculations. The momentum-dependent electron relaxation time τ^n,k_e is plotted in Fig. 5(d), where the τ_e value at the VBM and CBM is on the order of 10⁻¹ ps at 300 K. As the temperature increases, the τ_e value decreases significantly due to enhanced electron–phonon scattering. In the following subsection, as shown in eqn (6)–(9), the momentum-dependent τ_e is utilized to solve the electron BTE and provide accurate predictions for the electrical conductivity and Seebeck coefficient of Sb₂C₁₂.

3.4 Thermoelectric performance

First, we examined the Seebeck coefficient (S) and electrical conductivity (σ) for both p-type and n-type monolayer Sb₂C₁₂. The doping level or the carrier concentration is artificially introduced by using the rigid band approximation, where the additional carriers shift the Fermi level while keeping the intrinsic band structure unchanged. In our analysis, we focused on the scalar values of S and σ along two identical lattice axes. The Seebeck coefficient is positive for p-type doping and negative for n-type doping, and it is antisymmetric across the band gap, with similar |S| values at the same electron or hole concentration, as shown in Fig. 6(a). At 300 K, the |S| values are approximately 200 μV K⁻¹ when the doping level approaches the VBM and CBM, but a slightly lower doping level of 10¹² cm⁻² increases the S values to around 375 μV K⁻¹ for both doping types. Besides, within the doping level range of 5 × 10¹¹ to 10¹⁴ cm⁻², |S| exhibits an excellent linear relationship with the logarithm of carrier concentration and a positive correlation with the temperature, which is a general feature of the Seebeck coefficient in semiconductors. As another significant TE parameter, the electrical conductivity σ is plotted in Fig. 6(b). Due to the similar values of τ_e and DOS at the VBM and CBM, the σ values of p-type and n-type Sb₂C₁₂ are almost equivalent when the doping level approaches the band edges. However, its dependence on the carrier concentration and temperature is more sensitive and opposite to that of |S|. At higher carrier concentrations, the increased number of charge carriers enhances electrical transport, leading to a higher σ value. Conversely, as the temperature increases, stronger electron–phonon scattering significantly reduces the τ_e and σ values. At 300 K, the σ values for both p-type and n-type Sb₂C₁₂ are around 10⁵ S m⁻¹ at a carrier concentration of about 10¹³ cm⁻². Consequently, the high levels of S and σ contribute to the favorable TE performance of monolayer Sb₂C₁₂.

Fig. 6 (a) Absolute Seebeck coefficient and (b) electrical conductivity of p-type and n-type Sb₂C₁₂ at 300 K, 600 K and 900 K, respectively. (c) and (d) Power factors of p-type and n-type Sb₂C₁₂.

In TE devices, the power factor (PF = S²σ) is a key parameter for assessing the energy conversion efficiency. Due to the opposing dependence of S and σ on the carrier concentration, a high PF results from a good balance between these factors. As shown in Fig. 6(c) and (d), the PF peak for p-type doping near the VBM reaches approximately 3.1 × 10^–3 W m⁻¹ K⁻² at 300 K, while the n-type PF peak near the CBM is slightly lower, with a value of 2.8 × 10^–3 W m⁻¹ K⁻². At higher temperatures, the PF for both doping types decreases due to the stronger temperature dependence of σ. The PF of monolayer Sb₂C₁₂ is comparable to those of of other predicted 2D materials with good TE performance, such as monolayer MoS₂ (around 2 × 10^–3 W m⁻¹ K⁻²) and monolayer SnSe₂ (around 1 × 10^–2 W m⁻¹ K⁻²).^65,66

The electronic thermal conductivity (κ_e) is also an outcome of solving the electron BTE, as shown in Fig. 7(a). In contrast to the lattice thermal conductivity, κ_e makes a minor contribution to the total thermal conductivity at low doping levels, but its effect becomes significant when the carrier concentration exceeds 1 × 10¹³ cm⁻², as plotted in Fig. 7(b). Using the PF and total thermal conductivity, as shown in Fig. 7(c) and (d), we evaluated the figure-of-merit, the ZT value, for monolayer Sb₂C₁₂ at various doping levels and temperatures. Given that thermal conductivity is more sensitive to temperature than the PF, the ZT value increases with increasing temperature. At 300 K, the p-type ZT peak reaches 0.98 at a hole concentration of 5.35 × 10¹² cm⁻², while the n-type ZT peak is approximately 0.86 at an electron concentration of 5.05 × 10¹² cm⁻². At higher temperatures, the maximum ZT values for p-type and n-type peaks reach 1.77 and 1.58 at 600 K, and 2.55 and 2.26 at 900 K, respectively. The ZT value of monolayer Sb₂C₁₂ at room temperature surpasses those of many other predicted and measured 2D materials, such as graphyne and few-layer TMDs.^67,68 This exceptional TE performance is primarily attributed to the ultralow lattice thermal conductivity and high PF, stemming from the strong anharmonic effects of metal atoms and the steep slope of the DOS near the band edges.

Fig. 7 (a) Electronic thermal conductivity and (b) total thermal conductivity of Sb₂C₁₂. (c) ZT values for p-type and (d) n-type Sb₂C₁₂ as a function of carrier concentration.

To evaluate the general TE properties of g-MCs, we compared the TE performance of monolayer Sb₂C₁₂ with those of six other types of 2D monolayer g-MCs. As shown in Fig. 8(a), the Seebeck coefficients of monolayer 2D g-MCs are mainly in the range of 180–220 μV K⁻¹ at their optimal ZT values and their electrical conductivities are within the 10⁴–10⁵ S m⁻¹ range. Monolayer W₂C₁₂ is an exception with a very low σ value due to its narrow band gap of less than 0.25 eV, which is comparable to the energy of high-frequency phonons. This leads to strong EPI and a reduced electron relaxation time near the band edges, whereas all other candidates have a band gap larger than 0.40 eV. Consequently, as shown in Fig. 8(b), the PF of monolayer g-MCs is generally on the order of 10⁻³ W m⁻¹ K⁻² and monolayer Rh₂C₁₂ exhibits the highest PF value of around 3.7 × 10⁻³ W m⁻¹ K⁻² among the 2D g-MCs. Furthermore, all monolayer 2D g-MCs exhibit low lattice thermal conductivity at 300 K, underscoring the universality of the strong anharmonicity effects induced by the metal atoms.

Fig. 8 (a) Seebeck coefficient and electrical conductivity of 2D monolayer g-MCs at their optimal ZT values at 300 K. (b) Power factor and lattice thermal conductivity of g-MCs. (c) and (d) Maximum ZT values for p-type and n-type g-MCs and the corresponding carrier concentrations (doping levels) at 300 K.

To elucidate the ranking of κ_L among monolayer 2D g-MCs, we listed the averaged metal–carbon bond length and ICOHP in Table 1. The g-MCs were divided into three groups according to the element group of the metal atoms: (1) Sb and Bi, (2) Al and Ga, and (3) Rh and Ir. The κ_L values of monolayer Sb₂C₁₂ and Bi₂C₁₂ are lower than those of the others due to their remarkably longer metal–carbon bond length. The weaker bonds allow Sb and Bi atoms to vibrate more freely, resulting in stronger lattice anharmonicity in these two g-MCs. For 2D g-MCs with metal atoms from the same element group, their κ_L values are positively correlated with the ICOHP. A small absolute ICOHP signifies that the metal atom is weakly bound to the nearby carbons, which induces anharmonicity and leads to a low lattice thermal conductivity. Additional details on the κ_L ordering of the studied 2D g-MCs are provided in Table S1. We found that the ordering also closely tracks the corresponding acoustic phonon group velocities and relaxation times. Besides, as shown in Fig. 8(c) and (d), even though the ZT value of monolayer Sb₂C₁₂ is the highest among all 2D g-MCs, some other 2D g-MCs also show competitive TE performance. For example, both p-type and n-type Bi₂C₁₂ achieve ZT values of above 0.85 at room temperature. At 600 K, the ZT values of Bi₂C₁₂ increase to 1.24 for p-type and 1.25 for n-type, which are comparable to that of the well-known 3D TE material PbTe (around 1.30 at 600 K).⁶⁹ Overall, 2D g-MCs represent promising candidates for TE applications and their ultralow lattice thermal conductivity and high power factor are key factors in achieving high ZT values.

Table 1 The relationship between κ_L at 300 K, metal–carbon bond length and integrated crystal orbital Hamilton population (ICOHP) for the metal–carbon bond of 2D g-MCs

	Sb	Bi	Al	Ga	Rh	Ir
d M–C (Å)	2.18	2.28	1.91	1.94	1.94	1.93
–ICOHP (eV)	4.17	3.67	4.04	3.15	3.51	4.21
κ L (W m^−1 K^−1)	0.78	0.62	2.86	0.83	6.48	14.32

---

4 Conclusions

In summary, using first-principles calculations based on DFT and the semi-classical BTE, we theoretically investigated the electrical and thermal transport properties of 2D g-MCs and assessed their TE performance. We showed that the strong anharmonic effects induced by metal atoms in their unique structures lead to an ultra-low lattice thermal conductivity of 2D g-SbC and g-BiC. Further analysis of phonon scattering rates, weighted phase space, and the Grüneisen parameter reveals that 4-ph scattering in g-MCs cannot be ignored for phonon relaxation processes. The multi-valley electronic band structure and the steep slope of the DOS near the band edges contribute to a high Seebeck coefficient and power factor. Due to its ultra-low thermal conductivity and favorable electrical conductivity, the ZT values of p-type monolayer Sb₂C₁₂ reach 0.98 at 300 K and 2.55 at 900 K. The present work provides a systematic understanding of the TE performance of 2D g-MCs, highlighting the effectiveness of tuning the thermal conductivity of 2D materials with anharmonic multi-phonon scattering and offering new opportunities for future design of high-performance 2D TE materials.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the Supplementary Information. Data for this article, including the expressions for 3-phonon and 4-phonon scattering rates, detailed thermoelectric properties of other six types of 2D g-MCs, and the effect of defect on the lattice thermal conductivity, are available at DOI: https://doi.org/10.1039/d5nr02299a.

Acknowledgements

We acknowledge the support of the NUS academic research fund (A-8002944-00-00). The calculations were performed on the computational facilities of NUS HPC and the National Supercomputing Centre (NSCC), Singapore.

References

1. J. R. Sootsman, D. Y. Chung and M. G. Kanatzidis, Angew. Chem., Int. Ed., 2009, 48, 8616–8639.
2. M. S. Dresselhaus, G. Chen, M. Y. Tang, R. G. Yang, H. Lee, D. Z. Wang, Z. F. Ren, J. P. Fleurial and P. Gogna, Adv. Mater., 2007, 19, 1043–1053.
3. H. Goldsmid, Introduction to Thermoelectricity, Springer, 2009, vol. 1.
4. G. J. Snyder and E. S. Toberer, Nat. Mater., 2008, 7, 105–114.
5. J. He and T. M. Tritt, Science, 2017, 357, eaak9997.
6. G. Tan, L. D. Zhao and M. G. Kanatzidis, Chem. Rev., 2016, 116, 12123–12149.
7. X. L. Shi, J. Zou and Z. G. Chen, Chem. Rev., 2020, 120, 7399–7515.
8. R. Venkatasubramanian, E. Siivola, T. Colpitts and B. O'Quinn, Nature, 2001, 413, 597–602.
9. B. Zhu, X. Liu, Q. Wang, Y. Qiu, Z. Shu, Z. Guo, Y. Tong, J. Cui, M. Gu and J. He, Energy Environ. Sci., 2020, 13, 2106–2114.
10. Z. H. Zheng, X. L. Shi, D. W. Ao, W. D. Liu, M. Li, L. Z. Kou, Y. X. Chen, F. Li, M. Wei, G. X. Liang, P. Fan, G. Q. Lu and Z. G. Chen, Nat. Sustainability, 2022, 6, 180–191.
11. M. Buscema, M. Barkelid, V. Zwiller, H. S. J. van der Zant, G. A. Steele and A. Castellanos-Gomez, Nano Lett., 2013, 13, 358–363.
12. K. Hippalgaonkar, Y. Wang, Y. Ye, D. Y. Qiu, H. Zhu, Y. Wang, J. Moore, S. G. Louie and X. Zhang, Phys. Rev. B, 2017, 95, 115407.
13. M. Kayyalha, J. Maassen, M. Lundstrom, L. Shi and Y. P. Chen, J. Appl. Phys., 2016, 120, 134305.
14. M. Yoshida, T. Iizuka, Y. Saito, M. Onga, R. Suzuki, Y. Zhang, Y. Iwasa and S. Shimizu, Nano Lett., 2016, 16, 2061–2065.
15. Y. M. Zuev, W. Chang and P. Kim, Phys. Rev. Lett., 2009, 102, 096807.
16. S. G. Nam, D. K. Ki and H. J. Lee, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 245416.
17. P. Wei, W. Bao, Y. Pu, C. N. Lau and J. Shi, Phys. Rev. Lett., 2009, 102, 166808.
18. T. Li, A. D. Pickel, Y. Yao, Y. Chen, Y. Zeng, S. D. Lacey, Y. Li, Y. Wang, J. Dai, Y. Wang, B. Yang, M. S. Fuhrer, A. Marconnet, C. Dames, D. H. Drew and L. Hu, Nat. Energy, 2018, 3, 148–156.
19. J. Kang, Z. Wei and J. Li, ACS Appl. Mater. Interfaces, 2019, 11, 2692–2706.
20. Y.-M. Zhao, Z. Wang, J. Zhou, C. Zhang, S. Shin and L. Shen, J. Mater. Chem. C, 2024, 12, 14496.
21. R. Yan, J. R. Simpson, S. Bertolazzi, J. Brivio, M. Watson, X. Wu, A. Kis, T. Luo, A. R. Hight Walker and H. G. Xing, ACS Nano, 2014, 8, 986–993.
22. J. U. Lee, D. Yoon, H. Kim, S. W. Lee and H. Cheong, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 081419.
23. S. Chen, A. L. Moore, W. Cai, J. W. Suk, J. An, C. Mishra, C. Amos, C. W. Magnuson, J. Kang, L. Shi and R. S. Ruoff, ACS Nano, 2011, 5, 321–328.
24. S. Li, K. M. Yam, N. Guo, Y. Zhao and C. Zhang, npj 2D Mater. Appl., 2021, 5, 52.
25. K. M. Yam, Y. Zhang, N. Guo, Z. Jiang, H. Deng and C. Zhang, Nanotechnology, 2023, 34, 465706.
26. Z. Zheng, L. Qi, X. Luan, S. Zhao, Y. Xue and Y. Li, Nat. Commun., 2024, 15, 7331.
27. S. Zhang, Y. Kong, Y. Gu, R. Bai, M. Li, S. Zhao, M. Ma, Z. Li, L. Zeng, D. Qiu, Q. Zhang, M. Luo, L. Gu, Y. Yu, S. Guo and J. Zhang, J. Am. Chem. Soc., 2024, 146, 4433–4443.
28. S. He, B. Wu, Z. Xia, P. Guo, Y. Li and S. Song, Nanoscale Adv., 2023, 5, 2487–2492.
29. C. Bolding, W. Martin, T. Haraniya, J. I. Deagueros, X. Li, L. Heck, E. Kone, V. Desyatkin, Z. Huang, W. Linthicum, X. Gao, X. Ma and V. Rodionov, 2025.  DOI:10.26434/chemrxiv-2025-3qbbv.
30. V. G. Desyatkin, W. B. Martin, A. E. Aliev, N. E. Chapman, A. F. Fonseca, D. S. Galvão, E. R. Miller, K. H. Stone, Z. Wang, D. Zakhidov, F. T. Limpoco, S. R. Almahdali, S. M. Parker, R. H. Baughman and V. O. Rodionov, J. Am. Chem. Soc., 2022, 144, 17999–18008.
31. H. Wang, Y. Zhang, K. M. Yam, X. Tang, X.-S. Wang and C. Zhang, Mater. Today Electron., 2023, 6, 100073.
32. Y. Xiang, W. Chen, M. Wang, Z.-Z. Zhu, S. Wu and X. Cao, ACS Appl. Mater. Interfaces, 2024, 16, 23199–23208.
33. S.-L. Li, Y. Song, G. Tian, Q. Liu, L. Qiao, Y. Zhao and L.-Y. Gan, Fuel Process. Technol., 2024, 261, 108106.
34. S.-L. Li, G. Tian, Y. Chen, Y. Zhao, F. Cai and L. Qiao, Int. J. Hydrogen Energy, 2025, 144, 1043.
35. Y. Xiang, L. Li, Y. Li, Z.-Z. Zhu, S. Wu and X. Cao, Appl. Surf. Sci., 2022, 602, 154380.
36. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. De Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
37. G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet and N. Marzari, npj Comput. Mater., 2018, 4, 72.
38. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
39. S. Baroni, S. De Gironcoli, A. Dal Corso and P. Giannozzi, Rev. Mod. Phys., 2001, 73, 515–562.
40. W. Li, L. Lindsay, D. A. Broido, D. A. Stewart and N. Mingo, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 174307.
41. Z. Han, X. Yang, W. Li, T. Feng and X. Ruan, Comput. Phys. Commun., 2022, 270, 108179.
42. W. Li, J. Carrete, N. A. Katcho and N. Mingo, Comput. Phys. Commun., 2014, 185, 1747–1758.
43. J. Bardeen and W. Shockley, Phys. Rev., 1950, 80, 72–80.
44. C. Herring and E. Vogt, Phys. Rev., 1956, 101, 944–961.
45. S. Poncé, E. R. Margine, C. Verdi and F. Giustino, Comput. Phys. Commun., 2016, 209, 116–133.
46. G. K. H. Madsen, J. Carrete and M. J. Verstraete, Comput. Phys. Commun., 2018, 231, 140–145.
47. J. M. Ziman, Electrons and phonons: the theory of transport phenomena in solids, Oxford University Press, 2001.
48. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
49. T. Tadano and W. A. Saidi, Phys. Rev. Lett., 2022, 129, 185901.
50. O. Hellman, I. A. Abrikosov and S. I. Simak, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 180301.
51. I. Errea, M. Calandra and F. Mauri, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 064302.
52. J. Yue, J. Zheng, J. Li, X. Shen, W. Ren, Y. Liu and T. Cui, Mater. Today Phys., 2024, 46, 101517.
53. X. Ding, X. Jin, Z. Chang, D. Li, X. Zhou, X. Yang and R. Wang, Phys. Rev. B, 2024, 110, 054304.
54. T. Feng, L. Lindsay and X. Ruan, Phys. Rev. B, 2017, 96, 161201.
55. X. Yang, T. Feng, J. Li and X. Ruan, Phys. Rev. B, 2019, 100, 245203.
56. H. Yu, L. C. Chen, H. J. Pang, P. F. Qiu, Q. Peng and X. J. Chen, Phys. Rev. B, 2022, 105, 245204.
57. J. H. Seol, I. Jo, A. L. Moore, L. Lindsay, Z. H. Aitken, M. T. Pettes, X. Li, Z. Yao, R. Huang, D. Broido, N. Mingo, R. S. Ruoff and L. Shi, Science, 2010, 328, 213–216.
58. B. Abeles, Phys. Rev., 1963, 131, 1906–1911.
59. N. K. Ravichandran and D. Broido, Phys. Rev. X, 2020, 10, 021063.
60. L. Lindsay, D. A. Broido and T. L. Reinecke, Phys. Rev. Lett., 2013, 111, 025901.
61. S. Mukhopadhyay, L. Lindsay and D. S. Parker, Phys. Rev. B, 2016, 93, 224301.
62. T. Yue, Y. Zhao, J. Ni, S. Meng and Z. Dai, npj Comput. Mater., 2023, 9, 17.
63. Y. Pei, X. Shi, A. Lalonde, H. Wang, L. Chen and G. J. Snyder, Nature, 2011, 473, 66–69.
64. N. Wang, M. Li, H. Xiao, Z. Gao, Z. Liu, X. Zu, S. Li and L. Qiao, npj Comput. Mater., 2021, 7, 18.
65. J. Hong, C. Lee, J. S. Park and J. H. Shim, Phys. Rev. B, 2016, 93, 035445.
66. G. Li, G. Ding and G. Gao, J. Phys.: Condens. Matter, 2017, 29, 015001.
67. P. H. Jiang, H. J. Liu, L. Cheng, D. D. Fan, J. Zhang, J. Wei, J. H. Liang and J. Shi, Carbon, 2017, 113, 108–113.
68. D. Wickramaratne, F. Zahid and R. K. Lake, J. Chem. Phys., 2014, 140, 124710.
69. Y. Xiao and L. D. Zhao, npj Quantum Mater., 2018, 3, 55.

---