1. Introduction

Layered transition metal chalcogenides exhibit diverse physical properties^[1-3], as well as promising optoelectronic properties^[4-6], catalytic activity^[7], and thermoelectric performance^[8,9]. As a result, they have been extensively explored for applications in devices such as nanoelectronics^[10], optoelectronic devices^[4], electrochemical energy conversion^[11], and energy storage^[12,13]. Among all the transition metal chalcogenides, layered transition-metal sulfides have emerged as a particularly promising class of materials due to their tunable physical and chemical properties. By modifying the chemical composition and crystal structure, the electronic mobility and thermal conductivity of these materials can be substantially altered, providing a powerful strategy for designing high-performance materials with tailored functionalities.

The X₂PdY₆ (X = Nb, Ta; Y = S, Se) family represents a novel class of layered transition metal chalcogenides, characterized by a unique monoclinic crystal structure with the space group C2/m. The monolayer structure consists of PdY₄ tetrahedra and XY₇ polyhedra, which are stacked together through van der Waals interactions^[14]. This layered arrangement not only facilitates excellent exfoliation properties, enabling the preparation of monolayer or few-layer two-dimensional materials, but also imparts outstanding mechanical flexibility and unique electronic, optoelectronic, and thermal transport properties. For instance, Cho et al. demonstrated that a 4.64 nm thick Nb₂PdS₆ film exhibits a breakdown current density of 52 MA cm^-2 under high electric fields^[15]. Additionally, bulk Ta₂PdS₆ displays semiconducting behavior at ambient temperature and pressure, with a bandgap close to 0 eV^[4]. Remarkably, the bandgap of these materials can be tuned by varying the number of layers, transitioning from ~0 eV in the bulk to approximately 1 eV in the monolayer (1L) form^[4]. This tunability makes X₂PdY₆ materials highly promising for applications in electronic and optoelectronic devices^[16,17]. Recent studies have also uncovered intriguing superconducting properties in Ta₂PdS₆ and Ta₂PdSe₆. For example, Liu et al. reported superconducting behavior in Ta₂PdS₆ at a pressure of 82.0 GPa, with a critical temperature of up to 5.2 K^[18], while Yang et al. demonstrated superconductivity in Ta₂PdSe₆ at a critical pressure of P_c ~18.3 Gpa^[19]. Beyond these superconducting studies, Ta₂PdSe₆ has also demonstrated high-performance, polarization-sensitive terahertz photodetection^[20], while subsequent ultrahigh-pressure measurements revealed an anomalous enhancement of the superconducting critical temperature in both Ta₂PdS₆ and Ta₂PdSe₆^[21]. Furthermore, the layered structure of these materials, combined with their layer-dependent thermal and electrical transport properties, positions them as promising candidates for thermoelectric and thermal management applications^[9,22].

As a promising candidate for micro- and nano-electronic devices, efficient heat dissipation in X₂PdY₆ is critical for reliable device operation. Previous transport measurements have reported the room-temperature total TC and electrical transport properties of single-crystalline Ta₂PdS₆ and Ta₂PdSe₆ along the b-axis ([010]) direction^[8,23,24]. However, the phonon-mediated heat transport mechanisms, the full anisotropic lattice thermal conductivity (LTC), and the dimensional crossover from bulk to few-layer and monolayer X₂PdY₆ remain insufficiently understood. Therefore, a systematic investigation of the intrinsic lattice thermal transport behavior of ideal X₂PdY₆ and their layered nanostructures is still needed for evaluating their potential in thermal management and related applications.

In this study, we address this gap by training a machine learning potential (MLP) using the neuroevolution potential (NEP) model based on the dataset obtained by density functional theory (DFT) calculations. This NEP model accurately describes the interatomic interactions in bulk and layered X₂PdY₆. Using this model, we performed molecular dynamics (MD) simulations to investigate the thermal transport properties of bulk and layered X₂PdY₆. Our results demonstrate that the bulk X₂PdY₆ family exhibits anisotropic thermal transport properties, with the maximum and minimum LTCs along the in-plane [010] and the cross-plane [102] crystal directions, respectively. The LTC anisotropy ratios at 300 K for Nb₂PdS₆, Nb₂PdSe₆, Ta₂PdS₆, and Ta₂PdSe₆ reach 10.4, 8.1, 14.3, and 9.0, respectively. In layered structures, the LTC of X₂PdY₆ decreases with increasing layer thickness and approaches the corresponding bulk values when the number of layers reaches 5. These findings provide valuable insight into the tunability of thermal transport in the X₂PdY₆ family, paving the way for their potential application in thermal management.

2. Methodology

The bulk crystal of X₂PdY₆ belongs to the monoclinic crystal system and adopts a layered structure with the C2/m space group. As shown in Figure 1a, the unique layered crystal structure consists of PdY₄ tetrahedra and XY₇ polyhedra prisms. These structural units are alternately arranged, forming a stable layered lattice. The Y atoms can be categorized into four types: Y₁ and Y₂ are bonded to both Pd and X atoms, while Y₃ and Y₄ are bonded only to X atoms. In a unit cell, there are two atoms of types Y₁, Y₂, and Y₃, and one atom of Y₄. Figure 1b presents the projection of a single X₂PdY₆ layer, where the layers in the three-dimensional structure correspond to the (20$\overline{1}$) plane, and the crystallographic orientations within the layer align with the [20$\overline{1}$] and [010] directions. The periodic lengths along the [20$\overline{1}$] direction are 18.03 Å for Nb₂PdS₆, 18.84 Å for Nb₂PdSe₆, 18.27 Å for Ta₂PdS₆, and 18.94 Å for Ta₂PdSe₆ (Figure 1b).

Display Full Size

Download

Figure 1. (a) Crystal structure of bulk X₂PdY₆ (X = Nb, Ta ; Y = S, Se); (b) Projected structure of monolayer X₂PdY₆ (X = Nb, Ta; Y = S, Se).

To accurately describe the interatomic interactions in X₂PdY₆, we trained a MLP based on the NEP framework. MLPs are widely adopted for atomistic simulations^[25]. The training datasets were constructed through a combination of ab initio molecular dynamics (AIMD) simulations, random atom displacement, lattice deformation, and active learning Figure 2. The detailed composition of training datasets is listed in Table S1. We included a small number of high temperature AIMD configurations above the experimentally accessible thermal stability range of some materials (1,200 K and 1,600 K)^[15]. These structures are helpful for sampling larger displaced local environments and short-range repulsive regions of the potential energy surface, thereby reducing extrapolation during MD simulations at the target temperatures.

Display Full Size

Download

Figure 2. Training set workflow chart. DFT: density functional theory; AIMD: ab initio molecular dynamics; NEP: neuroevolution potential; NPT: number of particles, pressure, and temperature; NVT: number of particles, volume, and temperature.

All DFT calculations for the training datasets were performed using the vienna ab initio simulation package (VASP)^[26] with the projector augmented wave (PAW) method^[27]. The generalized gradient approximation (GGA) in the form of the Perdew-Burke-Ernzerhof (PBE) functional was employed as the exchange-correlation functional^[28]. A plane-wave basis set with a cutoff energy of 500 eV was used, and the total energy convergence criterion for electronic self-consistent iterations was set to 1 × 10^-7 eV. The Brillouin zone was sampled with a k-point density of 0.2 Å^-1, and a Gaussian broadening of 0.05 eV was applied. For AIMD simulations, we incorporate the third-generation D3 dispersion correction method proposed by Grimme et al.^[29,30] to account for interlayer van der Waals interactions. However, in the self-consistent calculations for the training datasets, the D3 dispersion correction was intentionally omitted, i.e., the NEP model was trained without dispersion correction for energy, forces, and virial stress. After training, the D3 energy correction was directly added to the NEP model to serve as the potential function for MD simulations^[31]. This approach not only improves the accuracy of the NEP model but also ensures a precise description of long-range van der Waals interactions. Phonon dispersion calculations were performed using the finite displacement method combined with the phonopy^[32] package, with 3 × 3 × 2 and 1 × 5 × 1 supercells employed for the bulk and two-dimensional structures X₂PdY₆, respectively.

After constructing the training datasets, we trained the MLP for X₂PdY₆ using the fourth-generation NEP framework within the graphics processing units molecular dynamics (GPUMD) package^[33]. The precision of the model was evaluated using the root mean square error (RMSE), which quantifies the deviation between the predicted values of the model and the reference data. The hyperparameters used for training the NEP model are provided in Table 1.

To simulate the LTCs of X₂PdY₆, we employed the GPUMD package in conjunction with the homogeneous nonequilibrium molecular dynamics (HNEMD) method. The LTC is calculated as^[34,35]:

(1)$\frac{{⟨J_{\mu }(t)⟩}_{ne}}{TV}=\sum _v{\kappa }_{\mu v}{F}_e^v$

where, κ_μv is the LTC tensor, T and V are the temperature and volume of the system, respectively, and $⟨J{⟩}_{ne}$ is the non-equilibrium heat current induced by the external driving force F_e^[36]. The HNEMD method allows a spectral decomposition of LTC to determine the contributions of phonons at different frequencies:

(2)$\kappa (\omega )=\frac{2}{TV{F}_e}{\int }_{-\infty }^{+\infty }dt{e}^{i\omega t}K(t)$

where $K(t)=\sum _i⟨W_i(0)\cdot v_i(t)⟩$ is the force-velocity correlation function, with W_i and v_i being the virial tensor and the velocity of atom i, respectively. Based on the spectral decomposition, we further calculated the phonon mean free path (MFP):

(3)$\lambda (\omega )=\frac{\kappa (\omega )}{G(\omega )}$

where G(ω) is the spectral thermal conductance in the ballistic limit, obtained from nonequilibrium molecular dynamics (NEMD) simulations^[37-39]. The setup for NEMD simulations is illustrated in Figure S1 and the corresponding descriptions.

Based on the relationship between the MFP and the diffusive transport LTC, the variation of the LTC with length can be obtained for each frequency:

(4)$\frac{1}{\kappa (\omega ,L)}=\frac{1}{\kappa (\omega )}(1+\frac{\lambda (\omega )}{L})$

The total LTC κ(L) as a function of system length was then obtained by integrating over all frequencies:

(5)$\kappa (L)={\int }_0^{\infty }\frac{d\omega }{2\pi }\kappa (\omega ,L)$

For MD simulations, we employed the velocity-Verlet algorithm with a time step of 1 fs to update the positions and velocities of the particles. For bulk X₂PdY₆, we constructed a 5 × 25 × 6 supercell (approximately 5.9 nm × 8.2 nm × 5.9 nm) by expanding the monoclinic unit cell along the a, b, and c lattice directions. For the two-dimensional structures, we built a 5 × 30 supercell (approximately 9.1 nm × 9.8 nm) by expanding the cell along the [20$\overline{1}$] and [010] crystal directions. These supercell sizes were chosen to ensure the convergence of the calculated LTC with respect to system size. Prior to the calculation of LTC, we equilibrated the system under the isothermal-isobaric ensemble (NPT ensemble) for 1 ns, followed by an additional 1 ns of equilibration under the canonical ensemble (NVT ensemble) to ensure the system was in an equilibrium state. Subsequently, we performed a 10 ns HNEMD simulation under the NVT ensemble to compute the LTC, with a driving force of 6 × 10^-5 Å^-1. The temperature of the system was controlled using the Nosé-Hoover Chain thermostat^[40]. The external driving force used in HNEMD enables efficient statistical sampling, and the running LTCs from independent trajectories show closely consistent convergence (Figure S2). We therefore averaged the results from three independent 10 ns simulations. Each trajectory was divided into two consecutive 5 ns blocks, yielding six block-averaged LTC estimates. The error bars represent the standard deviation across these six block averages.

The present NEP based HNEMD simulations describe only phonon mediated lattice thermal transport in ideal stoichiometric crystals. Electronic thermal conductivity (TC), electron-phonon scattering, and extrinsic scattering by defects or sample boundaries are not included. Therefore, the calculated LTC should not be equated with the experimentally measured total TC, particularly for the semimetallic compounds for which the electronic contribution may be appreciable.

3. Results and Discussion

Figure 3a presents the evolution of RMSE for energy, forces, and virial during the training of the Ta₂PdS₆ NEP model as a function of training generation. The RMSE values for all three metrics decrease significantly in the initial stages of training and stabilize after reaching 10⁵ generation, indicating that the model's fitting accuracy improves progressively and converges effectively. The precision of the NEP model is further illustrated by the parity plots for energy, forces, and virial stress, shown in Figure 3b,c,d. These plots show excellent diagonal alignment, confirming the model’s high fitting accuracy for these physical quantities. The corresponding RMSE values are 1.95 meV/atom for energy, 114.95 meV/Å for forces, and 18.73 meV/atom for virial stress. Similar results are observed for Nb₂PdS₆, Nb₂PdSe₆, and Ta₂PdSe₆, as shown in Figures S3,S4,S5. For these materials, the RMSE values for energy, forces, and virial stress are all below 2.1 meV/atom, 120.0 meV/Å, and 18.8 meV/atom, respectively. These results collectively demonstrate the high accuracy of the trained NEP models in describing the interatomic interactions in X₂PdY₆.

Display Full Size

Download

Figure 3. (a) Ta₂PdS₆ NEP model training iterative process of energy, force and virial; (b-d) Parity plot of energy, force and virial. NEP: neuroevolution potential; DFT: density functional theory; RMSE: root mean square error.

To further validate the reliability of the NEP model, we compared the lattice constants of X₂PdY₆ predicted by the NEP model with those obtained from DFT calculations and experimental values (Table 2). The lattice constants predicted by the NEP model for the four materials show excellent agreement with both the DFT and experimental results, with only minor deviations, confirming the accuracy of the NEP model training. This also underscores the reliability of the DFT method in describing the structural properties of X₂PdY₆. Figure S6 compares the energy-volume relationship predicted by the NEP model with the results from DFT calculations. The close agreement between the two further validates the accuracy of the trained NEP model.

Additionally, we predicted the phonon spectra of bulk X₂PdY₆ based on the NEP model and compared them with the results of DFT calculations. Figure 4 presents the phonon dispersion relations for the four systems in their bulk state. The phonon spectra predicted by the NEP model agree well with those obtained from DFT calculations, with no significant discrepancies. Although minor differences are observed in certain optical phonon branches, these local variations have a negligible impact on the LTC, as confirmed by our subsequent spectral decomposition of the LTC, where the high-frequency optical phonons contribute minimally to the LTC. These results demonstrate that the NEP model reliably predicts the phonon properties and LTC of X₂PdY₆ materials.

Display Full Size

Download

Figure 4. Comparison of phonon dispersion curves calculated by DFT and predicted by NEP. (a) Nb₂PdS₆; (b) Nb₂PdSe₆; (c) Ta₂PdS₆; (d) Ta₂PdSe₆. The red circle indicates phonon hybridization. DFT: density functional theory; NEP: neuroevolution potential.

The unit cell of X₂PdY₆ contains 9 atoms, resulting in 27 dispersion branches. Due to the similar crystal structures of the four materials, their dispersion curves exhibit similar shapes and characteristics. However, the phonon cutoff frequencies vary due to differences in atomic masses. Nb₂PdS₆ has the highest cutoff frequency of 12 THz. The replacement of Nb atoms with heavier Ta atoms slightly reduces the cutoff frequency, while substituting S atoms with Se atoms leads to a more pronounced decrease. Specifically, Nb₂PdSe₆ and Ta₂PdSe₆ exhibit cutoff frequencies of approximately 8.6 and 7.9 THz, respectively. The greater reduction in cutoff frequencies upon replacing S with Se, compared to replacing Nb with Ta, arises because the bands near the cutoff frequency are primarily dominated by vibrations involving S and Se atoms. Furthermore, all four materials exhibit longitudinal acoustic (LA) phonon hybridization with optical phonons, as indicated by the red circles in Figure 4. The hybridized LA mode contains significant optical phonon components (Figure S7), suggesting strong phonon-phonon scattering in these materials. This strong scattering contributes to their relatively low LTC.

After validating the precision of the NEP model, we used the HNEMD method (Equation 1) to calculate the LTCs of bulk X₂PdY₆. X₂PdY₆ belongs to the monoclinic crystal system with the C2/m space group, where the b axis is perpendicular to the a and c axes, but the a and c axes form an angle of approximately 115°. In our model, the a and b axes are aligned with the x and y axes, respectively, while the c axis is oriented at an angle to the z direction. Due to the lower symmetry of the monoclinic crystal system, the LTC tensor contains non-zero off-diagonal elements, reflecting anisotropic thermal transport characteristics and coupling between different crystallographic directions. Using the HNEMD method, we calculated the full LTC tensor by applying the driving forces sequentially along three orthogonal axes. This approach allows us to determine all components of the LTC tensor, significantly improving the precision of the calculations and enabling the determination of the LTC in any direction. As illustrated in Figure 5a, LTC in any direction (θ, ϕ) can be derived from the LTC tensor using the following transformation^[41,42]:

Display Full Size

Download

Figure 5. (a) Calculation of the polar coordinates of LTC; (b) Polar plot of LTC for Ta₂PdS₆ at 300 K; (c-e) Two-dimensional projection maps of LTC on the crystal planes {001} (c), {010} (d), and {100} (e). The units for LTC are all (Wm^-1 K^-1). LTC: lattice thermal conductivity.

(6)$\kappa (\theta ,\varphi )=\alpha \left(\begin{array}{lll}{\kappa }_{xx}& {\kappa }_{xy}& {\kappa }_{xz}\\ {\kappa }_{yx}& {\kappa }_{yy}& {\kappa }_{yz}\\ {\kappa }_{zx}& {\kappa }_{zy}& {\kappa }_{zz}\end{array}\right){\alpha }^T$

where α = (sinϕcosθ, sinϕsinθ, cosϕ). All crystallographic LTCs reported below were obtained from the above full tensor transformation without approximating a crystallographic direction by a Cartesian axis.

Figure S2 shows the variation of the average LTC of Ta₂PdS₆ at 300 K as a function of production time. Among the off-diagonal elements of the LTC tensor, κ_xy, κ_yx, κ_zy, and κ_yz are nearly zero, indicating that the heat flow in the y direction is largely decoupled from the x and z directions. In contrast, κ_xz and κ_zx exhibit significant non-zero values, implying strong coupling between heat flows along the x and z directions. This behavior is consistent with the symmetry characteristics of the monoclinic crystal system. The LTC tensor at 300 K for the four compounds is shown in Table S2. Based on the LTC tensor, we calculate the spatial variation of LTC for Ta₂PdS₆ (Figure 5b) and for Nb₂PdS₆, Nb₂PdSe₆, and Ta₂PdSe₆ (Figures S8,S9,S10). The polar plots of LTC reveal significant anisotropy in the LTC of these materials, with substantial variations along different crystallographic directions. The anisotropy of LTC is further illustrated by the projected LTC on the crystal planes {001}, {010}, and {100} (Figure 5c,d,e) and Figures S8,S9,S10). At 300 K, the ratios of maximum-to-minimum LTC on the {001}, {010}, and {100} planes for Ta₂PdS₆ are 12.8, 5.7, and 2.6, respectively. Given the layered structure of X₂PdY₆, we calculated the LTC along two in-plane directions ([010] and [20$\overline{1}$]) and the cross-plane direction ([102]), as shown in Table 3. In all four materials, the LTC is the largest along the in-plane [010] direction and the smallest along the cross-plane [102] direction due to the weak van der Waals interactions between layers. The LTC along the in-plane [20$\overline{1}$] direction lies between these two extremes.

Experimentally, the room-temperature total TC values measured along the b-axis ([010]) are approximately 14.0 W m^-1 K^-1 for Ta₂PdS₆ and 17.0 W m^-1 K^-1 for Ta₂PdSe₆^[8]. For Ta₂PdS₆, the calculated [010] LTC (27.59 ± 2.86 W m^-1 K^-1) is higher than the measured total TC. In contrast, the calculated [010] LTC of Ta₂PdSe₆ (15.58 ± 0.78 W m^-1 K^-1) is numerically close to, but lower than, the measured total TC. Since Ta₂PdSe₆ is semimetallic, its measured total TC contains a non-negligible electronic contribution^[23], and thus does not provide a direct experimental benchmark for the calculated LTC. The difference between the calculated and measured values may arise from both experimental and theoretical factors. Our simulations describe ideal stoichiometric and defect-free crystals and include only phonon-mediated heat transport. In contrast, residual point defects, nonstoichiometry, dislocations, and stacking imperfections in experimentally grown single crystals may introduce additional phonon scattering. In addition, the calculated absolute LTC is subject to the approximations in the DFT reference data, the machine-learning potential, and the classical MD treatment of anharmonic phonon transport. Nevertheless, the focus of this work is the intrinsic directional and thickness-dependent trends, including the pronounced thermal anisotropy and its phonon-level origin. The possible electronic contribution also depends strongly on the electronic structure, carrier concentration, transport direction, and sample quality. For semiconducting Ta₂PdS₆, the electronic contribution is expected to be negligible at room temperature^[24]. By contrast, semimetallic Ta₂PdSe₆ and Nb₂PdSe₆ can exhibit an electronic thermal conductivity of the order of several W m^-1 K^-1, potentially approaching 10 W m^-1 K^-1 in highly conducting crystals^[23]. Nb₂PdS₆ has been reported to exhibit metallic electrical transport in nanowire devices, suggesting a non-zero but sample-dependent electronic contribution^[15].

From the perspective of the atomic arrangement, the direction [010] exhibits the most regular periodicity, with a period ranging between 3.2 and 3.4 Å for the four materials. In contrast, the atomic arrangement along the [20$\overline{1}$] direction has a much larger periodicity (18-19 Å) and possesses asymmetric bonding between Ta/Nb and S atoms (Figure 1). These factors contribute to stronger anharmonicity in phonon transport along the [20$\overline{1}$] direction, resulting in lower LTC compared to the [010] direction. Among the four materials, the differences in LTC along the cross-plane direction are relatively small, whereas the in-plane LTC values vary significantly. Nb₂PdS₆ and Ta₂PdS₆ exhibit larger in-plane LTCs, with values exceeding 27.0 Wm^-1 K^-1 along the [010] direction and reaching 14.3 and 11.0 Wm^-1 K^-1 along the [20$\overline{1}$] direction, respectively. These differences in in-plane LTC also lead to notable variations in the LTC anisotropy ratios of the materials. Ta₂PdS₆ exhibits the highest anisotropy of LTC, with a maximum-to-minimum LTC ratio of 14.3, while Nb₂PdSe₆ exhibits the lowest LTC anisotropy, with a ratio of 8.1. These maximum-to-minimum LTC ratios are defined between the in-plane [010] and cross-plane [102] directions. This overall anisotropy primarily originates from the contrast between strong intralayer bonding and weak interlayer interactions.

The four materials share identical crystal structures, and the differences in their LTC are mainly attributed to variations in atomic masses and bond strengths. Both factors affect phonon dispersion curves, subsequently, changing phonon group velocities and scattering rates. Within the layers, X₂PdY₆ is bonded through the X-Y and Pd-Y bonds, and the bond strengths of these bonds vary depending on the choice of X and Y elements. To quantify the relative bond strengths in the four materials, we calculated the bond lengths and integrated crystal orbital Hamiltonian population (ICOHP) values for the X-Y and Pd-Y bonds. Based on bonding environments, the Y atoms were categorized into four types: Y₁, Y₂, Y₃, and Y₄. Y₁ and Y₂ are bonded to both X and Pd, while Y₃ and Y₄ are bonded only to X atoms (Figure 1). The ICOHP values were calculated using the Local Orbital Basis Suite Towards Electronic-Structure Reconstruction (LOBSTER)^[43] package, which provides a quantitative measure of bond strength; a more negative ICOHP value indicates a stronger bond. Table 4 lists the bond lengths and ICOHP values for the X-Y and Pd-Y bonds in Nb₂PdS₆, Nb₂PdSe₆, Ta₂PdS₆, and Ta₂PdSe₆. For a given material, the bond lengths generally follow the trend X-Y₁ < X-Y₂ < X-Y₃ < X-Y₄, and the ICOHP values exhibit the same trend, indicating that the bond strengths decrease in the order X-Y₁ > X-Y₂ > X-Y₃ > X-Y₄. All Y atoms form three chemical bonds: Y₁ and Y₂ atoms are bonded to two X atoms and one Pd atom, while Y₃ and Y₄ atoms are bonded to three X atoms. Since the Y-Pd bond is significantly weaker than the Y-X bond (as reflected by their less negative ICOHP values), the X-Y₁ and X-Y₂ bonds are stronger than the X-Y₃ and X-Y₄ bonds. These bonds of varying strengths are distributed symmetrically along the in-plane [010] direction, but asymmetrically along the [20$\overline{1}$] direction. This asymmetry results in greater anharmonicity in the vibrations of X atoms along the [20$\overline{1}$] direction, resulting in a relatively lower LTC compared to the [010] direction.

Among the different materials, the Pd-S bonds have shorter bond lengths and more negative ICOHP values than the Pd-Se bonds, indicating stronger Pd-S bonding. This difference arises because Se has more electron shells than S, resulting in lower electronegativity and weaker bonding with Pd. This trend is further supported by the Bader charges of Pd. In Nb₂PdS₆ and Ta₂PdS₆, Pd atoms lose 0.212 and 0.222 electrons, respectively, while in Nb2PdSe6 and Ta2PdSe6, Pd atoms lose only 0.013e and 0.034e electrons, respectively. The reduced electron transfer in selenides reflects weaker bonding between Pd and Se compared to Pd and S. The inferred X-Y bond strengths follow the order: Nb-S > Ta-S > Nb-Se > Ta-Se. This trend indicates that, for a given atom, the bond strength decreases as the atomic number of the other atom in the same group increases. This is because within the same group, larger atomic numbers correspond to weaker nuclear attraction for the outermost electrons. Electron localization analysis further confirms that the S and Se atoms, which have stronger electron-attracting abilities, exhibit stronger electron localization compared to other atoms, while Pd atoms show the strongest electron delocalization (Figure S11).

Differences in chemical bond strengths and atomic masses among the materials also lead to variations in phonon group velocities, which directly affect the LTC. Stronger chemical bonds result in larger interatomic force constants, leading to larger phonon group velocities. Conversely, larger atomic masses reduce phonon frequencies and, consequently, group velocities. Based on phonon dispersion relations, we calculate the phonon group velocities for the four materials, as shown in Figure 6. Nb₂PdS₆ exhibits the highest phonon group velocity, with the value of its low-frequency acoustic phonons reaching up to 6 km/s. This is attributed to its strongest Nb-S and Pd-S bonds, as well as its relatively small atomic mass. As chemical bonds weaken and atomic masses increase, the group velocities decrease in the order: Nb₂PdS₆ > Ta₂PdS₆ > Nb₂PdSe₆ > Ta₂PdSe₆. This trend aligns with the relative strengths of the chemical bonds and the observed LTC values. In materials with identical crystal structures, a decrease in group velocity corresponds to flatter phonon bands, which increases the phonon scattering rates. Flatter phonon bands provide more channels for phonon-phonon scattering, as it becomes easier to satisfy the energy and momentum conservation condition for scattering events. This increased scattering probability further reduces the LTC.

Display Full Size

Download

Figure 6. Phonon group velocities for (a) Nb₂PdS₆; (b) Nb₂PdSe₆; (c) Ta₂PdS₆; (d) Ta₂PdSe₆.

Figure 7 presents the normalized in-plane and cross-plane LTC as a function of temperature for all the four materials, with values normalized to the corresponding LTCs at 300 K. Interestingly, although the absolute values of LTC in different materials and crystal orientations differ significantly, their normalized LTCs are close to each other. The LTCs of the four materials follow a T^-1 relationship with temperature, indicating that three-phonon scattering dominates thermal transport in these systems.

Display Full Size

Download

Figure 7. The reduced LTC of bulk Nb₂PdS₆, Nb₂PdSe₆, Ta₂PdS₆, and Ta₂PdSe₆ along three crystal directions (corresponding to in-plane and out-of-plane directions) as a function of temperature. All LTC values have been normalized by their respective values at 300 K. LTC: lattice thermal conductivity.

To gain deeper insight into the phonon transport mechanisms, we performed spectral decomposition of the LTC and calculated the phonon MFP along the Cartesian axes x, y, and z using Equation 3 (Figure S12)^[44]. Due to limitations in the spectral decomposition method for non-diagonal elements, we restricted our analysis to the principal axes. While we can directly determine the MFP along the [010] direction (corresponding to the y axis), the [20$\overline{1}$] and [102] directions deviate from the z- and x- axes by 10°. The spectral decomposition procedure yields Cartesian-axis MFP spectra. We therefore report λ_x, λ_y, and λ_z explicitly. Because the x and z axes are close to [102] and [20$\overline{1}$], respectively, these spectra provide qualitative proxies for the corresponding crystallographic MFP trends. This qualitative correspondence is supported by the relatively small deviations between the Cartesian and tensor-transformed LTC values (Table S3). These deviations bound the error in the integrated LTC but do not constitute a frequency-resolved error bound for the MFP spectra. At 300 K, low-frequency phonons exhibit MFPs on the order of 10²-10³ nm. For each material, the MFP along the y-axis exceeds those along the z and x directions. The MFP hierarchy (y > z > x) mirrors the LTC trend ([010] > [20$\overline{1}$] > [102]), confirming that the anisotropy in thermal transport originates from fundamental differences in phonon scattering strengths along different crystallographic directions. This directional dependence reflects the underlying structural anisotropy of the monoclinic layered lattice, where strong in-plane bonding (particularly along [010]) allows longer phonon propagation compared to other directions.

Based on the phonon MFP derived from spectral decomposition, we calculated the length dependence of LTC at 300 K for bulk X₂PdY₆ through Equations 4 and 5 (Figure 8). The LTC components in the x and z directions, which are mainly contributed from the [102] and [20$\overline{1}$] directions respectively, show similar convergence behavior with increasing length. While the absolute values differ slightly from direct calculations along specific crystal directions, the numerical agreement supports the qualitative use of the Cartesian-axis MFPs as proxies for the [102] and [20$\overline{1}$] directions. Remarkably, the cross-plane (x direction) LTC requires approximately 1.2 μm to reach 90% of the converged LTC (dashed line, Figure 8). This characteristic length scale is comparable to that of the in-plane [010] direction and slightly exceeds the corresponding length for the in-plane 20$\overline{1}$] direction. The observed behavior suggests that while the interlayer van der Waals interactions are indeed much weaker than the in-plane covalent bonds, they are sufficient to provide an interlayer vibrational pathway through which low frequency phonons with relatively long MFPs contribute to cross-plane heat transport. The observed length dependence reflects finite MFP effects and convergence toward the intrinsic diffusive LTC in X₂PdY₆.

Display Full Size

Download

Figure 8. The dependence of the LTC of bulk X₂PdY₆ (X = Nb, Ta; Y = S, Se) on length at 300 K, with dashed lines indicating the cumulative LTC reaching 90% of the total LTC. (a) Along the x-axis; (b) Along the y-axis; (c) Along the z-axis. LTC: lattice thermal conductivity.

Given that all four materials are layered structures, we investigated how dimensionality affects thermal transport properties by examining the LTC from multilayer to monolayer transitions. In 2D structure simulations, the in-plane [20$\overline{1}$] and [010] crystal orientations (bulk notations) were set as the x and y directions, respectively. As shown in Figure 9a, Ta₂PdS₆ exhibits characteristic 2D thermal transport behavior where reducing layer thickness from bulk to monolayer increases LTC significantly due to suppressed interlayer phonon scattering and reduced Umklapp process^[45-47]. The LTCs of monolayer Ta₂PdS₆ along the [20$\overline{1}$] and [010] directions (bulk lattice notations) reach 30.7 ± 2.1 and 45.0 ± 3.3 Wm^-1 K^-1, respectively. This thickness dependence follows consistent trends across the four materials, with all systems approach the bulk limit within approximately 4-5 layers. The LTC along the [20$\overline{1}$] direction for monolayer Nb₂PdS₆, Nb₂PdSe₆, and Ta₂PdSe₆ reaches 50.7 ± 4.2, 25.5 ± 3.1, and 18.0 ± 4.2 Wm^-1 K^-1, respectively, while along the [010] direction, the LTC reaches 57.4 ± 4.7, 31.7 ± 1.3, and 25.6 ± 1.2 Wm^-1 K^-1, respectively. Notably, the greater LTC enhancement along the in-plane [20$\overline{1}$] direction compared to in-plane [010] reduces in-plane anisotropy ratios from bulk to monolayer. The anisotropy values in Nb₂PdS₆, Nb₂PdSe₆, Ta₂PdS₆, and Ta₂PdSe₆ are reduced from 2.0, 1.9, 2.5, and 2.7 in bulk to 1.1, 1.2, 1.5, and 1.4 in monolayer state, respectively.

Display Full Size

Download

Figure 9. (a) The variation of Ta₂PdS₆ LTC with different layer numbers, with bulk LTC included for comparison; (b,c) Spectral decomposition of LTC along the in-plane [20$\overline{1}$] and [010] crystal orientation (bulk notations) for 1-5 L. LTC: lattice thermal conductivity.

To complement our understanding of thickness-dependent thermal transport, we analyzed the spectral contributions of phonons to LTC in Ta₂PdS₆ from 1L to 5L (Figure 9b,c). Although the maximum frequency of Ta₂PdS₆ is as high as 12 THz, phonons above 4 THz (in-plane [20$\overline{1}$] direction) and 5 THz (in-plane [010] direction) contribute negligibly to the LTC in monolayer and multilayer structures, with all LTC contributions originating from low-frequency modes. Besides, reducing layer thickness preferentially enhances the contribution of these low-frequency phonons, particularly below 2 THz. This phenomenon is also observed in the other three materials (Figures S13,S14,S15). This frequency-selective enhancement directly correlates with our earlier observation of increasing LTC in thinner structures, revealing that suppressed interlayer scattering specifically benefits low-frequency phonon transport. The persistent T^-1 temperature dependence in monolayers (Figure S16) confirms that three-phonon scattering remains the dominant thermal resistance mechanism regardless of dimensionality, though interlayer van der Waals coupling in thicker samples introduces additional scattering that disproportionately affects low-frequency modes and reduces their MFPs.

The phonon MFP calculations reveal striking differences between monolayer and bulk structures (Figure S17). In monolayer X₂PdY₆, the MFP along both in-plane [20$\overline{1}$] and [010] directions show a remarkable order-of-magnitude increase compared to the bulk values. Around 0.1 THz, the in-plane [20$\overline{1}$] direction exhibits similar MFPs (~5 μm) across all four materials, while the in-plane [010] direction shows composition-dependent variations: Nb₂PdS₆ and Ta₂PdS₆ demonstrate comparable MFPs (6-7 μm) around 0.1 THz, whereas Nb₂PdSe₆ and Ta₂PdSe₆ show smaller but still substantial values (3 and 2 μm, respectively).

To quantify the origin of the reduced in-plane anisotropy, we evaluated the frequency-resolved LTC enhancement factor upon exfoliation. We defined a frequency-resolved enhancement factor along specific crystallographic α between monolayer and bulk structures, $R_{\alpha }(\omega )=\frac{{\kappa }_{\alpha }^{1L}(\omega )}{{\kappa }_{\alpha }^{bulk}(\omega )}$, for both in-plane directions, as shown in Figure 10. For the bulk [20$\overline{1}$] direction, we used the spectral LTC along the Cartesian z axis as a practical proxy. This approximation is supported by the full LTC tensor analysis: κ_zz differs from the tensor transformed ${\kappa }_{[20\bar{1}]}$ by only 1.4-2.5% across the four bulk compounds, with a maximum deviation of 2.5% (Table S3). In all four compounds, the monolayer-to-bulk enhancement is concentrated below 1 THz. Notably, this enhancement is substantially stronger along [20$\overline{1}$] than along [010], indicating that interlayer decoupling preferentially releases low-frequency phonon transport along the originally less conductive in-plane direction. Consequently, the LTC along [20$\overline{1}$] increases more rapidly than that along [010], leading to a reduced in-plane anisotropy in the monolayer limit.

Display Full Size

Download

Figure 10. Frequency-resolved spectral LTC enhancement factor, $R_{\alpha }(\omega )=\frac{{\kappa }_{\alpha }^{1L}(\omega )}{{\kappa }_{\alpha }^{\text{bulk}}(\omega )}$, along the [20 mathsrtart mathend] and [010] directions for (a) Nb₂PdS₆; (b) Nb₂PdSe₆; (c) Ta₂PdS₆; (d) Ta₂PdSe₆. For bulk [20$\overline{1}$], the Cartesian z-axis spectral LTC is used as a proxy; its integrated LTC differs from the exact tensor-transformed [20$\overline{1}$] value by no more than 2.5%.

The spectral energy density (SED) analysis^[44,48] for monolayer Ta₂PdS₆ along [20$\overline{1}$] and [010] directions demonstrates that the resolved LA and TA modes along [20$\overline{1}$] generally exhibit shorter lifetimes than their counterparts along [010] (Figure S18), indicating stronger scattering of both longitudinal and transverse acoustic phonons in this direction. The group velocities along [20$\overline{1}$] are generally lower than those along [010] (Figure S19), although the difference is considerably more pronounced for the optical branches than for the low frequency acoustic branches. Because the spectral LTC is dominated by low-frequency acoustic phonons, the shorter LA and TA lifetimes are likely the primary dynamical origin of the reduced LTC along [20$\overline{1}$], whereas the group-velocity difference makes a secondary contribution. These findings are consistent with enhanced anharmonic scattering associated with the longer structural periodicity and asymmetric X-Y bonding along [20$\overline{1}$] analyzed above.

The length dependence of LTC in monolayers further demonstrates exceptional phonon transport characteristics (Figure 11). The convergence length for LTC reaches 50-100 μm in both [20$\overline{1}$] and [010] directions across all systems, with Nb₂PdS₆ showing the longest convergence length. The characteristic length for achieving 90% of total LTC reaches or even exceeds 10 μm, several times larger than the corresponding bulk values. These extended length scales directly correlate with the dramatically increased MFPs in monolayers, confirming that low-frequency phonons with exceptionally long MFPs dominate thermal transport in these 2D systems. This behavior contrasts sharply with bulk samples where interlayer scattering severely limits phonon propagation distances at low-frequencies.

Display Full Size

Download

Figure 11. The dependence of the LTC of monolayer X₂PdY₆ (X = Nb, Ta; Y = S, Se) on length at 300 K, with dashed lines indicating the cumulative LTC reaching 90% of the total LTC. (a) Along the x-axis; (b) Along the y-axis. LTC: lattice thermal conductivity.

4. Conclusions

We developed an accurate NEP model trained on DFT data to investigate thermal transport in layered monoclinic X₂PdY₆ (X=Nb and Ta, Y=S and Se). Molecular dynamics simulations reveal strong anisotropic LTC in bulk crystals, with the maximum LTC occurring along the in-plane [010] direction (28.0, 15.8, 27.6, and 15.6 Wm^-1 K^-1 for Nb₂PdS₆, Nb₂PdSe₆, Ta₂PdS₆, and Ta₂PdSe₆ at 300 K, respectively) and the minimum occurring along the cross-plane [102] crystal direction (1.7-2.7 Wm^-1 K^-1), yielding anisotropy ratios of 10.4, 8.1, 14.3, and 9.0 for Nb₂PdS₆, Nb₂PdSe₆, Ta₂PdS₆, and Ta₂PdSe₆, respectively. Along the in-plane directions, the LTC along the [010] direction is much larger than that along the [20$\overline{1}$] direction. The smaller LTC along the [20$\overline{1}$] direction arises from asymmetric X-Y bonding and large periodicity along this crystallographic direction; both factors significantly enhance the phonon anharmonicity. The X-Y bond strength in X₂PdY₆ is inversely proportional to the atomic number of the X and Y elements (Nb-S > Ta-S > Nb-Se > Ta-Se) and directly correlated with LTC magnitudes. When exfoliated to 2D structures, all materials show enhanced LTC and reduced anisotropy due to suppressed low-frequency phonon scattering. Compared to their bulk structures, the MFP of low-frequency phonons in monolayers is increased by 5-8 times, causing the LTC to converge at system lengths as long as 50-100 μm. Our results establish intrinsic LTC trends, anisotropy, and phonon-scattering mechanisms in layered X₂PdY₆. These findings are relevant to the phonon component of heat dissipation and provide a basis for thermal management design. However, the present simulations include only phonon-mediated transport in ideal crystals. Electronic thermal transport and extrinsic scattering mechanisms are not included, which might be important in semimetallic compounds such as Ta₂PdSe₆ and Nb₂PdSe₆.

Authors contribution

Wu X: Investigation, methodology, formal analysis, writing-review & editing.

Huang J, Li W, Guo Y, Zhang W: Investigation, formal analysis.

Cheng Y: Conceptualization, supervision.

Zhang G: Project administration, supervision.

Xiong S: Project administration, conceptualization, supervision, writing review & editing.

Conflict of interest

Gang Zhang is an Editorial Board Member and Shiyun Xiong is a Youth Editorial Board Member of Thermo-X. The other authors declare no conflicts of interest.

Ethical approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Funding

This work is supported by the National Natural Science Foundation of China (Grant No. 12304059), the Basic and Applied Basic Research Foundation of Guangdong Province (Grant Nos. 2024A1515012635, 2024A1515010521 and 2022A1515110572), Guangzhou Municipal Science and Technology Program (Grant No. 2025A04J3566).

Copyright

© The Author(s) 2026.