Machine learning-guided discovery of multifunctional Nb₂CX₂ MXenes (X = S, Se, Te): insights into mechanical properties and thermal conductivity

Abstract

The lattice thermal conductivity of Nb₂CX₂ (X = S, Se, Te) MXenes has been investigated with machine learning force fields (MLFF). To validate the precision of the MLFF method, structural parameters and elastic properties were analyzed and compared between density functional theory (DFT) and MLFF results. The results showed minimal deviations, demonstrating that the MLFF method provides reliable accuracy at the DFT level. The calculated Young's modulus is high, indicating the strong mechanical properties of these materials and their potential for energy storage applications. Interestingly, functionalization with heavier terminal groups reduces the lattice thermal conductivity (κ_L), with Nb₂CS₂, Nb₂CSe₂, and Nb₂CTe₂ exhibiting values of 16.93, 12.05, and 6.01 W m⁻¹ K⁻¹ at 300 K, respectively. This trend is further explained through phonon group velocities, phonon lifetimes, and Grüneisen parameters, which provide deeper insight into the underlying thermal transport mechanisms. Our findings highlight the key role of surface functionalization in tailoring thermal transport, paving the way for optimized MXenes in thermal management. This study also establishes MLFF as a highly efficient, accurate, and cost-effective tool for accelerating materials discovery in energy and thermal applications.

---

1. Introduction

The swift emergence of two-dimensional (2D) graphene-like materials has attracted considerable attention due to their unique electronic, magnetic, optical, catalytic, and thermal properties.^1–3 A notable advancement in this field is the discovery of MXenes, a new class of 2D hexagonal transition metal carbides and nitrides. MXenes are produced by selectively etching the “A” element from MAX phases, where M is a transition metal, A is an element from group IIIA or IVA, and X represents carbon (C) or nitrogen (N).^4,5 MXenes can be easily functionalized with surface groups like O, F, and OH. This functionalization significantly alters their properties and increases their versatility.^5,6 For example, the adsorption of O or F groups on the surface of MXenes such as Sc₂C, Ti₂C, Hf₂C, Zr₂C, and Mo₂C can transform their metallic state into a semiconducting state.⁷ The O-functionalized M₂CO₂ (M = Zr, Hf) MXenes have demonstrated great potential as photocatalysts for water splitting. Their ultrahigh carrier mobilities and suitable band structures make them particularly promising for this application.⁸ It is evident that MXenes have demonstrated a great deal of possibilities for a range of uses. Li/Na-ion batteries, hybrid electrochemical capacitors, catalysts, spintronic devices, and flexible electronics are a few of their applications.^9–16 In Li-ion batteries, surface functionalization typically increases the open-circuit voltage while decreasing the transport mobility of Li-ions in the electrode materials used.^7,17,18

Two-dimensional Nb₂CS₂ MXene is one of the most widely used and can be experimentally produced as monolayers from the bulk phase.¹⁹ Nb₂CS₂ bulk phase is a layered structure composed of 2D slabs held together by weak van der Waals (vdW) forces. Experimental studies have also demonstrated that Nb₂CS₂ can intercalate with various metal atoms and complex organic molecules.²⁰ This highlights the high feasibility of producing Nb₂CS₂ monolayers using exfoliation techniques. Interestingly, it was discovered that Nb₂CS₂ is a superconductor at low temperatures.²¹ The pristine and TM-decorated Nb₂CS₂ monolayer has also been the subject of numerous systematic investigations using ab initio calculations. These investigations have predicted that the material will make an excellent anode for metal-ion batteries,²² single-atom catalysts for a low-temperature CO oxidation reaction,²³ and chemical sensors.²⁴ Additionally, delaminated Nb₂CT_x and carbon nanotube composites had a greater capacity than multilayer Nb₂CT_x when utilized as Li-ion battery electrodes, with T_x representing OH, O, and F-based functional groups.²⁵ According to another theoretical investigation, Nb₂CS₂ exhibits a superconducting critical temperature that is extremely similar to the experimental data and the MXenes series Nb₂CT₂ (T = O, S, Se, or Te) are intrinsic phonon-mediated superconductors.²⁶ While extensive research has focused on the impact of surface functionalization on the electronic, magnetic, and electrochemical properties of MXenes, experimental measurements of their thermal conductivity remain unexplored. This is primarily due to the challenges associated with synthesizing freestanding MXenes. Thus, the influence of surface functionalization on the thermal conductivity and elasticity of MXenes is still unclear.

Thermal transport plays a crucial role in condensed matter physics, where heat is primarily carried by phonons and electrons.²⁷ At temperatures near or above room temperature, phonon–phonon interactions significantly influence lattice thermal conductivity (κ_L), making the accurate prediction of phonon anharmonic behavior a long-standing challenge.^28–30 Traditionally, the phonon Boltzmann transport equation (PBTE)^31,32 combined with density functional theory (DFT) calculations has been the standard approach for predicting κ_L. However, this method is computationally expensive, as it requires extensive DFT calculations to determine second- and third-order interatomic force constants (IFCs), typically limiting anharmonic interactions to third-order terms. While these interactions capture essential phonon anharmonicity, they may not fully describe phonon quasiparticles at elevated temperatures, where higher-order effects become increasingly important.^33,34 Moreover, this approach faces scalability challenges when applied to large or complex systems. Recent advancements in machine learning potentials, such as Gaussian approximation potentials (GAP) and deep potentials (DP), have improved the balance between computational efficiency and accuracy.^35,36 Among these, on-the-fly machine-learned force fields (MLFFs) integrated into VASP provide a highly efficient alternative, allowing real-time interpolation of DFT data for IFC computations without relying on empirical fitting parameters.³⁷ Motivated by these developments, this study leverages MLFFs to efficiently compute the phononic contribution to κ_L in Nb₂CX₂, providing a cost-effective yet reliable approach to thermal transport modeling.

In light of the aforementioned claims, we use the newly established on-the-fly MLFF in VASP in order to expedite the prediction of elastic characteristics and lattice thermal properties in MXenes Nb₂CX₂ with its diverse surface terminal groups (X = S, Se, Te). To explore the factors affecting changes in lattice thermal conductivity due to various surface terminal groups, we examine the heat capacity, phonon group velocity, and phonon lifetime of Nb₂CX₂. The outcomes validate the approach as a swift and precise way to assess the thermal conductivity of the material, which will significantly aid in the future development of thermoelectric and thermal management systems.

2. Computational method

The on-the-fly machine learning force fields (MLFFs), integrated into the Vienna Ab Initio Simulation Package (VASP), were generated during ab initio molecular dynamics (AIMD) simulations.³⁷ Using a time step of 1 fs, AIMD simulations were conducted in the canonical (NVT) ensemble in this study at a temperature of 200 K for 5 ps. Utilizing its stochastic characteristics for efficient phase space sampling, the Langevin thermostat³⁸ was used to control the temperature of the system. The choice of 200 K for MLFF training was justified by a detailed analysis of root mean square error (RMSE) and mean absolute error (MAE) values across different temperatures (300 K and 400 K) for Nb₂CS₂ (see Fig. S1 in the ESI†). The results show that RMSE and MAE values at 200 K are not only comparable to those at higher temperatures but also slightly lower, indicating improved accuracy and stability in MLFF training. Given this trend, we reasonably extend the selection of 200 K for MLFF training of the other two systems, ensuring consistency and computational efficiency in our approach. A 5 × 5 × 1 supercell of Nb₂CX₂ MXenes made up the training cell. The Bayesian error threshold on the forces was set at 2 × 10⁻² eV Å⁻¹ and it was programmed to automatically update as the accuracy of the simulation increased with the number of steps.^39–41 For the angular and radial descriptors, we established cutoff radii of 5.0 and 8.0 Å, respectively. A Gaussian function of 0.5 Å width was used to expand atomic distributions.⁴² In MLFFs, the polynomial power ξ of the kernel and the weight β for the descriptors were set to 4 and 0.1, respectively. Local reference configurations were sparsified using the CUR algorithm with a threshold of 10⁻⁹.⁴³

To ensure the accuracy of the MLFFs at the DFT level, DFT calculations in VASP served as the foundation for our computations.^44,45 The Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) was used to characterize the electron exchange–correlation functional.⁴⁶ The plane-wave expansion uses a kinetic energy cutoff of 520 eV. The vacuum space must be larger than 10 perpendicular to the sheet in order to prevent self-image interactions. The Brillouin zone is integrated using a 10 × 10 × 1 Monkhorst–Pack k-point mesh.⁴⁷ A total energy difference and a Hellmann–Feynman force of less than 10⁻⁶ eV and 0.01 eV Å⁻¹, respectively, are the convergence criterion.

The finite displacement method inside Phonopy packages was used to determine the phonon dispersion and harmonic interatomic force constants (IFCs).⁴⁸ For all Nb₂CX₂, the atomic displacement distance for harmonic IFCs computations was fixed at 0.01 eV Å⁻¹ using a supercell 5 × 5 × 1. The lattice thermal conductivity (κ_L) was determined through the direct solution of the linearized phonon Boltzmann transport equation (PBTE). The anharmonic interatomic force constants (IFCs) were computed using the same methodology as for the harmonic IFCs, as implemented in the Phono3py package.^32,49,50 For the anharmonic IFC calculations, the atomic displacement amplitude was set to 0.03 eV Å⁻¹, while maintaining the same supercell size as used for the harmonic IFC calculations. Additionally, a dense 81 × 81 × 1 q-mesh was employed to ensure convergence.

The reliability of MLFF-generated forces was further quantified by computing the root mean square error (RMSE) of forces and stresses, which remained relatively small (detailed in the Results and discussion section). The consistency of MLFF predictions with DFT results for optimized structural parameters and phonon properties further validates the accuracy of our computational setup.

The elastic properties, including Young's modulus (Y_n) and Poisson's ratio (ν_n), in an arbitrary direction (n⃑), are obtained using the computed C_ij values and the following equations:⁵¹

Here, Δ = C₁₁C₂₂ − C₁₂², c = cos θ, s = sin θ, and n⃑ = cos θx̂ + sin θŷ

The Debye temperature (θ_D), average sound velocity (V_m), and transverse and longitudinal sound velocities (V_T and V_L) have been calculated by the following equations:

Here, ρ represents the density.

where ħ is the reduced Planck constant, k_B is the Boltzmann constant, N is the number of atoms in the unit cell, and A its area.

The expression for lattice thermal conductivity (κ_L) is represented by:

Here, c_qj, v_qj, and τ_qj represent the constant-volume heat capacity, group velocity, and relaxation time associated with the phonon mode (q, j), respectively, while V₀ denotes the unit cell volume.

The gradient of the phonon dispersion curve yields these velocities, v_g, which can be written as follows:

The wave vector is represented by k, the angular frequency is indicated by ω, and the frequency associated with k by ω(k). Furthermore, for a given phonon mode, the Grüneisen parameter, which measures the shift in phonon frequencies (ω_qj) in relation to changes in material volume (V), can be expressed mathematically as follows:

3. Results and discussion

3.1 Structural properties

As seen in Fig. 1, we began by structurally optimizing the crystal structure of Nb₂CX₂ (where X = S, Se, and Te). Every optimized geometry of Nb₂CX₂ has the symmetry of the hexagonal space group p3̄m1 (no. 164). Table 1 displays the determined bond lengths and lattice parameter a from the DFT computation. Here, we have adopted the energetically most favourable configurations for Nb₂CX₂ based on the previous literature.²⁶ Our findings regarding the structural parameters for Nb₂CX₂ are also exactly in line with the earlier report.^22,24,26 The phonon dispersions of Nb₂CX₂ have been calculated in order to confirm their dynamic stability. Here, we take the phonon band structure into account and apply DFT to validate the dynamic stability of the MXenes. In the phonon spectrum across the Brillouin zone, positive frequencies represent dynamic stability. Dynamic instability is indicated with imaginary frequencies (negative values) in the spectrum, which have no restoring effect against atom displacement. Fig. 2 displays the phonon band structure that was estimated for Nb₂CX₂. All MXenes Nb₂CX₂ with X = S, Se, and Te surface functionalization are dynamically stable, as evidenced by the absence of imaginary modes across the Brillouin zone.

Fig. 1 Side view of the optimized structure: (a) unit cell and (b) of 5 × 5 × 1 supercell Nb₂CX₂, (c) the top view structure of the Nb₂CX₂ monolayer in a 5 × 5 × 1 supercell (X = S, Se, and Te).

Table 1 Comparison of structural parameters obtained using DFT and MLFF for Nb₂CX₂ (X = S, Se, Te)

System	DFT			MLFF			Ref.
	a = b (Å)	Nb–C (Å)	Nb–X (Å)	a = b (Å)	Nb–C (Å)	Nb–X (Å)
Nb2CS2	3.278	2.237	2.504	3.274	2.227	2.516	Present work
	3.28	2.24	2.51	—	—	—	24
	3.28	2.34	2.51	—	—	—	22
	3.29	2.23	2.51	—	—	—	26
Nb2CSe2	3.327	2.245	2.635	3.321	2.237	2.647	Present work
	3.34	2.26	2.64	—	—	—	26
Nb2CTe2	3.435	2.275	2.831	3.393	2.263	2.849	Present work
	3.43	2.28	2.83	—	—	—	26

---

Fig. 2 The phonon band structures of Nb₂CX₂ explored through DFT and MLF study.

The structures of Nb₂CX₂ were reoptimized using X = S, Se, and Te surface functionalization by machine learning force field (MLFF) in order to verify the computational accuracy of our MLFF. In comparison to the DFT results, Table 1 displays the structural parameters, including the lattice constants, with a negligible divergence. As determined by the MLFF and DFT approaches, the per-atom force and stress tensor of Nb₂CX₂ are shown in Fig. 3. Interestingly, there is a good linear connection (slope of about 0.5) between the force and stress tensor values from FMLP and those from DFT. Additionally shown in Fig. 3 are the associated RMSE and MAE for the forces and stress tensors, which show remarkably low values at 200 K. Because of these low RMSE and MAE values, MLFF is equivalent to traditional DFT computations in terms of its great accuracy in replicating the potential energy surface. This validates the ability of MLFF-driven accelerated computations to accurately and consistently predict a variety of physical attributes. Furthermore, the total energy plotted against volume for both DFT and MLFF is shown in Fig. S2 (ESI†), which shows that MLFF and DFT have very similar total energies. Additionally, we computed the phonon band structure of Nb₂CX₂ using MLFF, as shown in Fig. 2 and compared the findings with DFT. A fairly comparable pattern was found between the two. These results validate the accuracy of the developed MLFF in reproducing the optimized structures of Nb₂CX₂.

Fig. 3 The comparison of per-atom force (a)–(c) and stress tensor (d)–(f) of Nb₂CX₂ calculated by MLFF and DFT.

3.2 Mechanical properties

The elastic constants of a crystalline solid are the primary parameters that indicate its stiffness, stability, and anisotropy. They provide crucial information about the forces operating in the material and act as a link between its mechanical and thermal characteristics and its dynamic behavior. MXenes have three separate linear elastic constants (C₁₁, C₁₂, and C₆₆) by virtue of their 2D hexagonal crystal symmetry. C₁₂ demonstrates the stress that results along the x-axis when a uniaxial strain is applied along the crystallographic y-axis, whereas C₁₁ illustrates the resistance of the material to linear compression along the crystallographic x-axis. Isotropic conditions are enforced in trigonal and hexagonal symmetries, such as the Cauchy relation, C₆₆ = (C₁₁ − C₁₂)/2, and the relation C₁₁ = C₂₂.⁵¹ In hexagonal crystals, this symmetry criterion for C₆₆ is important. Regarding elastic characteristics, the [100] plane in the [010] direction is subject to the shear constant C₆₆, whereas the shear constant (C₁₁ − C₁₂)/2 is pertinent for the [110] plane in the [110] direction. A particular aspect of transverse isotropy in the [001] zone is that the elastic shear constant is constant in all planes and orientations. This suggests that the hexagonal crystals are elastically isotropic within the xy plane since the elastic constants reflect invariance with regard to arbitrary rotations around the z-axis. For Nb₂CX₂, the estimated elastic constants are shown in Table 2. The materials were evidently mechanically stable, as they fulfilled the stability conditions C₁₁ > 0 and C₁₁ > |C₁₂|.⁵² Across all MXenes, our results demonstrate an upward trend with few outliers for both DFT and MLFF techniques. The C₁₁ and C₁₂ values observed in Nb₂CS₂ and Nb₂CSe₂ compared to Nb₂CTe₂ can be attributed to the stronger bonding between Nb and S or Se. This stronger bonding arises from the higher electronegativity of S and Se compared to Te, resulting in greater lattice stiffness and higher elastic constants.^53,54 Te also contributes more to lattice vibrations (phonons) due to its higher atomic mass, which enhances the susceptibility of the lattice to mechanical deformation. This phonon-induced softening leads to a reduction in the elastic constants C₁₁ and C₁₂ for Nb₂CTe₂.^55,56

Table 2 Comparison of mechanical properties from DFT and MLFF calculations for Nb₂CS₂, Nb₂CSe₂, and Nb₂CTe₂

Properties	DFT			MLFF
	Nb2CS2	Nb2CSe2	Nb2CTe2	Nb2CS2	Nb2CSe2	Nb2CTe2
C 11 (N m^−1)	217.051	210.237	191.190	245.744	215.602	185.873
C 22 (N m^−1)	217.051	208.946	187.298	245.744	214.496	183.262
C 12 (N m^−1)	44.479	47.395	34.421	81.905	75.349	23.140
C 66 (N m^−1)	86.286	82.128	73.090	81.919	100.257	81.345
Y (N m^−1)	207.936	199.552	184.993	218.445	216.426	182.992
ν	0.204	0.225	0.180	0.333	0.341	0.124
ρ (g m^−2)	4.67 × 10^−3	6.16 × 10^−3	7.36 × 10^−3	4.69 × 10^−3	6.18 × 10^−3	7.54 × 10^−3
G (N m^−1)	86.286	82.128	73.0900	81.919	100.257	81.346
B (N m^−1)	130.766	128.817	112.806	163.825	147.976	104.507
V L (m s^−1)	6461.888	5559.789	4772.309	6930.423	6009.280	4685.274
V T (m s^−1)	4298.449	3651.365	3151.30	4179.322	4027.756	3284.599
V m (m s^−1)	4700.376	3999.236	3449.181	4621.053	4400.053	3565.908
Θ D (K)	933.054	782.181	653.390	918.429	862.129	653.863
M. stability	Stable	Stable	Stable	Stable	Stable	Stable

---

The Poisson ratio ν = C₁₂/C₁₁ and Young's modulus Y = (C₁₁² − C₁₂²)/C₁₁ are determined in this instance and are both independent of the angle θ, verifying planar isotropy. These numbers describe the mechanical performance of MXenes, which is compiled in Table 2. The deformation of the material in directions perpendicular to the applied load is determined by ν, whereas Y establishes the stiffness. Since wearable and flexible devices must bend and deform without breaking, MXenes with a lower Y and a higher ν are widely desired. On the other hand, larger Y and lower ν are preferable for high-strength applications, including energy storage devices, to make sure that the material withstands stress without experiencing undue lateral deformation.⁵⁷ According to our calculations, the MXenes Nb₂CS₂, Nb₂CSe₂ and Nb₂CTe₂ all exhibit potential for energy storage devices because of their higher Y and lower ν. The (n⃑) dependence of Y and ν for the anisotropic structures in polar coordinates with DFT data could be evaluated using eqn (1) and (2) in conjunction with the elastic constants from Table 2, as seen in Fig. 4. In this case, a perfectly round shape in the Y_n and ν_n graphs for Nb₂CS₂ indicates a fully isotropic elastic behavior. However, Fig. 4 shows that Nb₂CTe₂ is significantly more anisotropic than Nb₂CSe₂. The isotropic and anisotropic behavior is further confirmed by the relationship between elastic constants, where C₁₁ = C₂₂ for Nb₂CS₂, indicating isotropy, while C₁₁ ≠ C₂₂ for Nb₂CSe₂ and Nb₂CTe₂, highlighting their anisotropic nature (Table 2). Additionally, as expected, Fig. 4(b) shows that Nb₂CTe₂ has a Poisson ratio that is significantly lower than Nb₂CS₂ and Nb₂CSe₂.

Fig. 4 Polar diagram for the (a) Young's modulus and (b) Poisson ratio of Nb₂CX₂.

Moreover, the Debye temperature (θ_D), average sound velocity (V_m), and transverse and longitudinal sound velocities (V_T and V_L) are crucial factors in comprehending the thermal and elastic characteristics of materials. The mechanical and thermal behavior of the material can be better understood by analyzing them in terms of bulk modulus (B) and shear modulus (G) (eqn (3)).^58,59 The average sound velocity V_m offers a macroscopic understanding of the behavior of wave propagation in the material and can be computed by eqn (4). Again, eqn (5) can be used to determine the Debye temperature (θ_D) from the average sound velocity. Strong atomic bonds and high heat conductivity are indicated by high θ_D values.⁶⁰ These parameters were computed using DFT and MLFF data (Table 2), and both methods demonstrate a steady trend with little variations for all MXenes.

3.3 Lattice thermal conductivity

The phonon Boltzmann transport equation (PBTE) can be solved with the harmonic second-order and anharmonic third-order force constants to find the lattice thermal conductivity κ_L of the monolayers. The expression for lattice thermal conductivity (κ_L) under the single-mode relaxation time approximation (RTA) is given by eqn (6)

The phonon lifetime and phonon relaxation time are the same in RTA. Here, we used MLFF to determine the lattice thermal conductivity of the MXenes, and Fig. 5 displays the κ_L against the temperature plot. As expected, κ_L falls as the temperature rises. Temperature-induced increases in phonon–phonon scattering provide an explanation for this occurrence. According to our findings, at room temperature (300 K), Nb₂CTe₂ has a lower κ_L of 6.01 W m⁻¹ K⁻¹ than Nb₂CS₂ and Nb₂CSe₂, which have respective values of 16.93 and 12.05 W m⁻¹ K⁻¹. At 1000 K, these values have been dropped to 4.78, 3.53, and 1.75 W m⁻¹ K⁻¹ for Nb₂CS₂, Nb₂CSe₂ and Nb₂CTe₂, respectively. The greater size of Te compared to S and Se may also contribute more structural distortions and flaws into the lattice, which justifies the decreased thermal conductivity of Nb₂CTe₂. These flaws significantly lower the thermal conductivity of Nb₂CTe₂ by serving as phonon scattering centers.⁶⁰

Fig. 5 Calculated lattice thermal conductivities (κ_L) of Nb₂CX₂.

We investigated the phonon lifetime, group velocity, and Grüneisen parameters as functions of frequency to gain a deeper understanding of the mechanisms influencing κ_L. All the MXenes being assessed exhibited similar constant volume heat capacities c_qj, even while Nb₂CX₂ has a slightly lower c_qj at temperatures underneath 400 K [Fig. 6]. We will therefore be primarily concerned with the phonon group velocities, phonon lifetime, and Grüneisen parameters. Properties of heat transport are largely determined by the group velocities of phonons (eqn (7)). Lattice thermal conductivity is affected by low-lying optical phonon branches in the phonon dispersion curve, as well as the group velocity, which is represented by the slope of acoustic phonon branches.⁶¹ Lattice thermal conductivity is generally correlated with a lower group velocity and frequency. Lattice thermal conductivity further reduces with increasing temperature in materials with low-lying optical modes because of phonon–phonon interactions, especially Umklapp scattering. In the majority of materials, the acoustic branches dominate heat transport because they have larger group velocities than the optical branches, which have lower group velocities and higher frequencies. The group velocity of the optical branches is smaller than that of the acoustic branches, as seen in Fig. 7. Our analysis revealed that the highest group velocities for acoustic branches for the Nb₂CS₂, Nb₂CSe₂, and Nb₂CTe₂ are around 42.31 km s⁻¹, 38.112 km s⁻¹, and 30.438 km s⁻¹, respectively. This is in line with the higher κ_L values found in the Nb₂CS₂ as opposed to Nb₂CSe₂ and Nb₂CTe₂.

Fig. 6 Constant volume heat capacities (c_qj) of Nb₂CX₂.

Fig. 7 The phonon group velocity (v_g) of Nb₂CX₂.

Low-lying optical modes with large group velocities become more important at high temperatures, since these optical phonons scatter acoustic modes more efficiently. Lattice thermal conductivity decreases as a result of the phonon lifetimes being shortened by the resulting increased scattering rate. In essence, a phonon lifetime is the amount of time it can travel before scattering; a longer lifetime denotes fewer scattering occurrences and, hence, more effective heat transport. Fig. 8 shows that the phonon lifetime in the instance of Nb₂CX₂ changes with frequency. The long-wavelength vibrations of lower-frequency phonons contact fewer contaminants and imperfections, resulting in longer lifetimes. Our results show that the acoustic branches in the Nb₂CS₂ have maximum phonon lifetimes of about 10³ ps, which is consistent with a larger κ_L value than that of Nb₂CSe₂ and Nb₂CTe₂.

Fig. 8 Phonon lifetimes of Nb₂CX₂ calculated at 300 K.

Our study also concentrated on the Grüneisen parameter (γ_qj), which is essential to comprehending lattice thermal conductivity, especially because it provides information on anharmonic phonon interactions, which is the main cause of thermal resistance in materials (eqn (8)).⁶² The phonon lifetime is inversely related to the square of the averaged Grüneisen parameter γ² in the context of continuum theory.⁶³ This anharmonicity is measured by the Grüneisen parameter since κ_L is mostly controlled by phonon–phonon scattering processes, which are known to be largely anharmonic. Greater values of γ² indicate more robust anharmonic phonon interactions, which raise phonon–phonon scattering and lower the thermal conductivity of the material.^64,65

In order to conduct a thorough evaluation of phonon anharmonicity and its consequences for thermal transport properties, we used the methods described in ref. 66 to directly compute γ_qj from third-order force constants. Fig. 9 shows the mode-Grüneisen parameters (γ_qj) of Nb₂CX₂. A positive (γ_qj) value means that the frequency of a certain phonon mode (q,j) falls with increasing volume. This is typical behavior when force constants diminish and interatomic distances grow. However, some materials, like Cu₂O, have notable negative (γ_qj) values, which contribute to negative thermal expansion across particular temperature ranges.⁶⁷ In our study, Nb₂CS₂ predominantly shows positive γ_qj values. It is noteworthy that Nb₂CTe₂ and Nb₂CSe₂ have some fairly negative γ_qj values for phonon modes in both the low-frequency and high-frequency domains, specifically below 1 THz and over 40 THz. Furthermore, we found that, in comparison to Nb₂CS₂, Nb₂CSe₂ and Nb₂CTe₂ typically show larger γ_qj values, especially in the acoustic phonon modes below 10 THz. This pattern is consistent with the lower lattice thermal conductivity κ_L values found for Nb₂CSe₂ and Nb₂CTe₂, demonstrating that greater atomic masses lead to lower group velocities and phonon frequencies, which in turn have an impact on lattice thermal conductivity.⁶⁸

Fig. 9 The mode-Grüneisen parameters (γ_qj) of Nb₂CX₂.

4. Conclusion

In this study, we investigated the lattice thermal conductivity (κ_L) of Nb₂CX₂ MXenes (X = S, Se, Te) using DFT-based machine learning force field (MLFF) potentials. The structural parameters and elastic properties obtained from DFT and MLFF calculations show excellent agreement, confirming the reliability of MLFFs. The mechanical stability of Nb₂CX₂ (X = S, Se, Te) MXenes is validated, as all variants satisfy the stability criteria C₁₁ > 0 and C₁₁ > |C₁₂|. Our findings also reveal that these MXenes, specifically Nb₂CS₂, Nb₂CSe₂, and Nb₂CTe₂, are highly suitable for energy storage applications due to their high Young's modulus (Y). Notably, our MLFF simulations demonstrate that the lattice thermal conductivity (κ_L) decreases with the addition of heavier terminal groups (S → Te). Among the variants, Nb₂CSe₂ and Nb₂CTe₂ exhibit particularly low κ_L values at room temperature, measured at approximately 12.047 W m⁻¹ K⁻¹ and 6.014 W m⁻¹ K⁻¹, respectively. This reduction of κ_L is further substantiated by analyzing phonon group velocities, lifetimes, and Grüneisen parameters, which confirm the role of heavier functional groups in suppressing phonon-mediated heat transport.

Beyond its physical insights, this study underscores the advantages of MLFF in computational materials science. By integrating MLFF with DFT, we achieve a substantial reduction in computational cost while maintaining high accuracy in phonon transport predictions. Compared to traditional DFT-based approaches, MLFF enables efficient large-scale simulations without compromising precision, positioning it as a powerful tool for high-throughput material screening. These results highlight the transformative potential of MLFF in accelerating the discovery and optimization of materials for next-generation thermal management technologies.

Data availability

The data that support the findings of this study are available within the article and its ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This project is supported by Kasetsart University Research and Development Institute, KURDI: FF(KU-SRIU)7.67. S. K. J. greatly acknowledges the Postdoctoral Fellowship from the Kasetsart University. A. B. is funded by the 31st Science & Technology Research Grant (STRG) from the Thailand Toray Science Foundation. T. T. is supported by the National Research Council of Thailand (NRCT), No. N41A661103. We wish to thank the voucher-type2 program from NSTDA Supercomputer Center (ThaiSC), project no. pv812003, for providing computing resources for this work.

References

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D.-E. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Electric field effect in atomically thin carbon films, Science, 2004, 306, 666–669.
2. F. Bonaccorso, L. Colombo, G. Yu, M. Stoller, V. Tozzini, A. C. Ferrari, R. S. Ruoff and V. Pellegrini, Graphene, related two-dimensional crystals, and hybrid systems for energy conversion and storage, Science, 2015, 347, 1246501.
3. Y. Guo, F. Pan, M. Ye, Y. Wang, Y. Pan, X. Zhang, J. Li, H. Zhang and J. Lu, Interfacial properties of stanene–metal contacts, 2D Mater., 2016, 3, 035020.
4. B. Anasori, M. R. Lukatskaya and Y. Gogotsi, Mxenes, Jenny Stanford Publishing, 2023, pp. 677–722.
5. M. Naguib, V. N. Mochalin, M. W. Barsoum and Y. Gogotsi, 25th anniversary article:MXenes: a new family of two-dimensional materials, Adv. Mater., 2014, 26, 992–1005.
6. M. Naguib, M. Kurtoglu, V. Presser, J. Lu, J. Niu, M. Heon, L. Hultman, Y. Gogotsi and M. W. Barsoum, Mxenes, Jenny Stanford Publishing, 2023, pp. 15–29.
7. M. Khazaei, M. Arai, T. Sasaki, C.-Y. Chung, N. S. Venkataramanan, M. Estili, Y. Sakka and Y. Kawazoe, Novel electronic and magnetic properties of two-dimensional transition metal carbides and nitrides, Adv. Funct. Mater., 2013, 23, 2185–2192.
8. Z. Guo, J. Zhou, L. Zhu and Z. Sun, MXene: a promising photocatalyst for water splitting, J. Mater. Chem. A, 2016, 4, 11446–11452.
9. X. Wang, X. Shen, Y. Gao, Z. Wang, R. Yu and L. Chen, Atomic-scale recognition of surface structure and intercalation mechanism of Ti3C2X, J. Am. Chem. Soc., 2015, 137, 2715–2721.
10. S. Kajiyama, L. Szabova, K. Sodeyama, H. Iinuma, R. Morita, K. Gotoh, Y. Tateyama, M. Okubo and A. Yamada, Sodium-ion intercalation mechanism in MXene nanosheets, ACS Nano, 2016, 10, 3334–3341.
11. M. Naguib, O. Mashtalir, J. Carle, V. Presser, J. Lu, L. Hultman, Y. Gogotsi and M. W. Barsoum, Two-dimensional transition metal carbides, ACS Nano, 2012, 6, 1322–1331.
12. T.-D. T. Carbide, Cation Intercalation and High Volumetric Capacitance of, Science, 2013, 1241488, 341.
13. H. Tang, Q. Hu, M. Zheng, Y. Chi, X. Qin, H. Pang and Q. Xu, MXene–2D layered electrode materials for energy storage, Prog. Nat. Sci.:Mater. Int., 2018, 28, 133–147.
14. J. Ran, G. Gao, F.-T. Li, T.-Y. Ma, A. Du and S.-Z. Qiao, Ti3C2 MXene co-catalyst on metal sulfide photo-absorbers for enhanced visible-light photocatalytic hydrogen production, Nat. Commun., 2017, 8, 13907.
15. Z. Guo, N. Miao, J. Zhou, B. Sa and Z. Sun, Strain-mediated type-I/type-II transition in MXene/Blue phosphorene van der Waals heterostructures for flexible optical/electronic devices, J. Mater. Chem. C, 2017, 5, 978–984.
16. H. Zhang, Z. Fu, R. Zhang, Q. Zhang, H. Tian, D. Legut, T. C. Germann, Y. Guo, S. Du and J. S. Francisco, Designing flexible 2D transition metal carbides with straincontrollable lithium storage, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E11082–E11091.
17. Y. Xie, M. Naguib, V. N. Mochalin, M. W. Barsoum, Y. Gogotsi, X. Yu, K.-W. Nam, X.-Q. Yang, A. I. Kolesnikov and P. R. Kent, Role of surface structure on Li-ion energy storage capacity of two-dimensional transition-metal carbides, J. Am. Chem. Soc., 2014, 136, 6385–6394.
18. Q. Tang, Z. Zhou and P. Shen, Are MXenes promising anode materials for Li ion batteries? Computational studies on electronic properties and Li storage capability of Ti3C2 and Ti3C2X2 (X= F, OH) monolayer, J. Am. Chem. Soc., 2012, 134, 16909–16916.
19. H. Boller and K. Hiebl, Quaternary pseudo-intercalation phases Tx [Nb2S2C](T V, Cr, Mn, Fe, Co, Ni, Cu) and metastable Nb2S2C formed by topochemical synthesis, J. Alloys Compd., 1992, 183, 438–443.
20. K. Sakamaki, H. Wada, H. Nozaki, Y. Ōnuki and M. Kawai, Topochemical formation of van der Waals type niobium carbosulfide 1T-Nb2S2C, J. Alloys Compd., 2002, 339, 283–292.
21. K. Sakamaki, H. Wada, H. Nozaki, Y. Ōnuki and M. Kawai, van der Waals type carbosulfide superconductor, Solid State Commun., 2001, 118, 113–118.
22. Y. Jing, J. Liu, Z. Zhou, J. Zhang and Y. Li, Metallic Nb2S2C monolayer: a promising two-dimensional anode material for metal-ion batteries, J. Phys. Chem. C, 2019, 123, 26803–26811.
23. M. Li, T. Li and Y. Jing, Nb2S2C Monolayers with Transition Metal Atoms Embedded at the S Vacancy Are Promising Single-Atom Catalysts for CO Oxidation, ACS Omega, 2023, 8, 31051–31059.
24. S. Lakshmy, N. Kalarikkal and B. Chakraborty, Transition Metal (Cu, Pd, and Ag)-Modified Nb2S2C Monolayer for Highly Efficient Catechol Sensing: A First-Principles Investigation, Langmuir, 2024, 40, 13819–13833.
25. A. Byeon, A. M. Glushenkov, B. Anasori, P. Urbankowski, J. Li, B. W. Byles, B. Blake, K. L. Van Aken, S. Kota and E. Pomerantseva, et al., Lithium-ion capacitors with 2D Nb2CTx (MXene)–carbon nanotube electrodes, J. Power Sources, 2016, 326, 686–694.
26. S.-Y. Wang, C. Pan, H. Tang, H.-Y. Wu, G.-Y. Shi, K. Cao, H. Jiang, Y.-H. Su, C. Zhang and K.-M. Ho, et al., Straintronic Effect on Phonon-Mediated Superconductivity of Nb2CT2 (T = O, S, Se, or Te) MXenes, J. Phys. Chem. C, 2022, 126, 3727–3735.
27. T. M. Tritt, Thermal conductivity: theory, properties, and applications, Springer Science & Business Media, 2005.
28. C. Herring, Role of low-energy phonons in thermal conduction, Phys. Rev., 1954, 95, 954.
29. P. G. Klemens, Solid state physics, Elsevier, 1958, vol. 7, pp. 1–98.
30. P. Klemens, The scattering of low-frequency lattice waves by static imperfections, Proc. Phys. Soc., Sect. A, 1955, 68, 1113.
31. A. Ward, D. Broido, D. A. Stewart and G. Deinzer, Ab initio theory of the lattice thermal conductivity in diamond, Phys. Rev. B:Condens. Matter Mater. Phys., 2009, 80, 125203.
32. A. Togo, L. Chaput and I. Tanaka, Distributions of phonon lifetimes in Brillouin zones, Phys. Rev. B:Condens. Matter Mater. Phys., 2015, 91, 094306.
33. T. Feng, L. Lindsay and X. Ruan, Four-phonon scattering significantly reduces intrinsic thermal conductivity of solids, Phys. Rev. B, 2017, 96, 161201.
34. L. Xie, J. Feng, R. Li and J. He, First-principles study of anharmonic lattice dynamics in low thermal conductivity AgCrSe2: evidence for a large resonant four-phonon scattering, Phys. Rev. Lett., 2020, 125, 245901.
35. A. P. Bartók, J. Kermode, N. Bernstein and G. Csányi, Machine learning a generalpurpose interatomic potential for silicon, Phys. Rev. X, 2018, 8, 041048.
36. T. Wen, L. Zhang, H. Wang, E. Weinan and D. J. Srolovitz, Deep potentials for materials science, Mater. Futures, 2022, 1, 022601.
37. P. Liu, C. Verdi, F. Karsai and G. Kresse, Phase transitions of zirconia: Machine-learned force fields beyond density functional theory, Phys. Rev. B, 2022, 105, L060102.
38. W. G. Hoover, A. J. Ladd and B. Moran, High-strain-rate plastic flow studied via nonequilibrium molecular dynamics, Phys. Rev. Lett., 1982, 48, 1818.
39. J. Skilling, Maximum Entropy and Bayesian Methods: Cambridge, England, 1988, Springer Science & Business Media, 2013, vol. 36.
40. D. J. MacKay, Bayesian interpolation, Neural Comput., 1992, 4, 415–447.
41. R. Jinnouchi and R. Asahi, Predicting catalytic activity of nanoparticles by a DFT-aided machine-learning algorithm, J. Phys. Chem. Lett., 2017, 8, 4279–4283.
42. A. P. Bartók, M. C. Payne, R. Kondor and G. Csányi, Gaussian approximation potentials:The accuracy of quantum mechanics, without the electrons, Phys. Rev. Lett., 2010, 104, 136403.
43. R. Jinnouchi, F. Karsai, C. Verdi, R. Asahi and G. Kresse, Descriptors representing two-and three-body atomic distributions and their effects on the accuracy of machinelearned inter-atomic potentials, J. Chem. Phys., 2020, 152, 234102.
44. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169.
45. G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 1996, 6, 15–50.
46. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 1996, 77, 3865.
47. H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B, 1976, 13, 5188.
48. A. Togo and A. Seko, On-the-fly training of polynomial machine learning potentials in computing lattice thermal conductivity, J. Chem. Phys., 2024, 160, 211001.
49. A. Togo, L. Chaput, T. Tadano and I. Tanaka, Implementation strategies in phonopy and phono3py, J. Phys.: Condens. Matter, 2023, 35, 353001.
50. L. Chaput, Direct solution to the linearized phonon Boltzmann equation, Phys. Rev. Lett., 2013, 110, 265506.
51. E. Cadelano, P. L. Palla, S. Giordano and L. Colombo, Elastic properties of hydrogenated graphene, Phys. Rev. B:Condens. Matter Mater. Phys., 2010, 82, 235414.
52. V. Wang, G. Tang, Y.-C. Liu, R.-T. Wang, H. Mizuseki, Y. Kawazoe, J. Nara and W. T. Geng, High-throughput computational screening of two-dimensional semiconductors, J. Phys. Chem. Lett., 2022, 13, 11581–11594.
53. J. Häglund, A. F. Guillermet, G. Grimvall and M. Körling, Theory of bonding in transition-metal carbides and nitrides, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 48, 11685.
54. Z. Wu, X.-J. Chen, V. V. Struzhkin and R. E. Cohen, Trends in elasticity and electronic structure of transition-metal nitrides and carbides from first principles, Phys. Rev. B:Condens. Matter Mater. Phys., 2005, 71, 214103.
55. F. Zeng, W.-B. Zhang and B.-Y. Tang, Electronic structures and elastic properties of monolayer and bilayer transition metal dichalcogenides MX2 (M = Mo, W; X = O, S, Se, Te): a comparative first-principles study, Chin. Phys. B, 2015, 24, 097103.
56. N. K. Nepal, L. Yu, Q. Yan and A. Ruzsinszky, First-principles study of mechanical and electronic properties of bent monolayer transition metal dichalcogenides. Physical Review, Materials, 2019, 3, 073601.
57. C. Zhang, B. Anasori, A. Seral-Ascaso, S.-H. Park, N. McEvoy, A. Shmeliov, G. S. Duesberg, J. N. Coleman, Y. Gogotsi and V. Nicolosi, Transparent, flexible, and conductive 2D titanium carbide (MXene) films with high volumetric capacitance, Adv. Mater., 2017, 29, 1702678.
58. O. L. Anderson, A simplified method for calculating the Debye temperature from elastic constants, J. Phys. Chem. Solids, 1963, 24, 909–917.
59. R. Hill, The elastic behaviour of a crystalline aggregate, Proc. Phys. Soc., Sect. A, 1952, 65, 349.
60. Z. Tong, S. Li, X. Ruan and H. Bao, Comprehensive first-principles analysis of phonon thermal conductivity and electron-phonon coupling in different metals, Phys. Rev. B, 2019, 100, 144306.
61. E. S. Toberer, L. L. Baranowski and C. Dames, Advances in thermal conductivity, Annu. Rev. Mater. Res., 2012, 42, 179–209.
62. L. Lindsay and D. Broido, Three-phonon phase space and lattice thermal conductivity in semiconductors, J. Phys.: Condens. Matter, 2008, 20, 165209.
63. J. M. Ziman, Electrons and phonons: the theory of transport phenomena in solids, Oxford University Press, 2001.
64. E. S. Toberer, A. Zevalkink and G. J. Snyder, Phonon engineering through crystal chemistry, J. Mater. Chem., 2011, 21, 15843–15852.
65. G. A. Slack, The thermal conductivity of nonmetallic crystals, Solid State Phys., 1979, 34, 1–71.
66. N. Glebko and A. J. Karttunen, Lattice thermal conductivity of TiS 2, ZrS 2, and HfS 2: Periodic trends studied by dispersion-corrected hybrid density functional methods, Phys. Rev. B, 2019, 100, 024301.
67. J. Linnera and A. Karttunen, Ab initio study of the lattice thermal conductivity of Cu₂O using the generalized gradient approximation and hybrid density functional methods, Phys. Rev. B, 2017, 96, 014304.
68. G. A. Slack, Thermal conductivity of pure and impure silicon, silicon carbide, and diamond, J. Appl. Phys., 1964, 35, 3460–3466.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4tc05430j

---