Recently, the cold source field-effect transistor (CSFET) has emerged as a promising solution to overcome Boltzmann tyranny in its ballistic regime, offering a steep-slope subthreshold swing (SS) of less than 60 mV/decade. However, challenges arise due to scattering, particularly from inelastic scattering, which can lead to significant degradation in SS through cold carrier thermalization. In this study, we delve into the theoretical investigation of the electronic excitation/relaxation dynamic process using the state-of-the-art nonadiabatic molecular dynamics (NAMD) method. The mixed quantum-classical NAMD proves to be a powerful tool for comprehensively analyzing cold carrier thermalization and transfer processes in semiconductor Si, as well as metallic silicides (NiSi $ 2$ and CoSi $ 2$). The approach of mixed quantum-classical NAMD takes into account both carrier decoherence and detailed balance, enabling the calculation of thermalization factors, relaxation times, scattering times, and scattering rates at various energy levels. The thermalization of carriers exhibits a gradual increase from low to high energy levels. Achieving partial thermalization from the ground state to reach the thermionic current window occurs within a sub-100 fs time scale. Full thermalization across the entire energy spectrum depends sensitively on the barrier height, with the scattering rate exponentially decreasing as the energy of the out-scattering state increases. Notably, the scattering rate of NiSi $ 2$ and CoSi $ 2$ is two orders of magnitude higher than that of Si, attributed to their higher density of states compared to Si. This study not only provides insights into material design for low-power tunnel field-effect transistors but also contributes valuable information for advancing CSFET in emerging technologies.

## I. INTRODUCTION

The evolution of integrated circuits in accordance with Moore’s Law has not only propelled advancements in the traditional information industry but has also played a pivotal role in shaping emerging fields, such as artificial intelligence, the Internet of Things, and big data. As devices continue to undergo miniaturization and integration, their physical dimensions are approaching the limits set by fundamental physics. Moreover, a significant obstacle arises in the form of a power consumption wall due to the scaling limit of supply voltage, which must maintain a sufficient on/off ratio. This scaling limit becomes a critical challenge in transistor technology, as discussed by Datta *et al*.^{1} To address these challenges and surpass the existing limits, researchers are exploring new transistor switching mechanisms aimed at achieving a steep subthreshold swing (SS). Notable approaches include tunnel FETs (TFETs)^{2–5} and negative-capacitance FETs (NCFETs).^{6–8} A recent addition to this repertoire is the novel steep-slope cold source field-effect transistors (CSFETs),^{9,10} including Dirac source FET as a special case of CSFET,^{11–14} that has been proposed. These innovative CSFETs incorporate various designs of density of states engineering schemes to effectively filter cold carriers in the source region. Predictions of sub-60 mV/decade SS have been made across different material systems, including graphene/MoS $ 2$,^{15–17} graphene/Cd $ 3$C $ 2$/T-VTe $ 2$/H-VTe $ 2$/H-TaTe $ 2$-InSe,^{18} metal Heusler alloy-WS $ 2$/MoSe $ 2$,^{19} KTlCl $ 4$/Au $ 2$S,^{20} artificial type-III heterojunction of Si,^{21–23} and more, without accounting for any scattering effects. This research signifies a promising direction in the pursuit of advanced transistor technologies to overcome existing barriers and continue the progress of integrated circuits.

In the path-finding stage, various types of steep-slope CSFETs have been predicted using the non-equilibrium Green’s function method with density functional theory (NEGF/DFT) in the ballistic transport regime, excluding Dirac source FETs, which have been experimentally demonstrated.^{11–14} However, the carrier transport mechanisms and the comprehensive performance analysis of versatile CSFETs remain largely unexplored. In advanced ultra-scaled transistors, even when the transport length is below 20 nm, diverse scattering effects can degrade performance.^{24–26} Scattering processes typically fall into two categories: elastic scattering, such as Coulomb scattering, which alters the carrier momentum distribution without affecting the energy distribution, and inelastic scattering, including electron–phonon and electron–electron interactions, leading to rethermalization by transitioning cold carriers to higher energy states. Understanding the thermal, electrical, and thermoelectric transport properties of materials considering scattering effects has become increasingly crucial. While first-principles electron–phonon interaction calculations and accurate carrier relaxation time estimations have been reported,^{27,28} applying these methods to large systems poses computational challenges. An alternative is the nonadiabatic molecular dynamics (NAMD) approach,^{29} widely used for studying carrier dynamic processes. To alleviate simulation burdens, several NAMD algorithms are implemented in a mixed quantum-classical (MQC) fashion, such as mean-field Ehrenfest dynamics (MFE)^{30,31} and fewest switches surface hopping (FSSH).^{32,33} However, both MFE and FSSH have limitations in addressing the disparities between quantum and classical mechanics. MFE lacks detailed balance and the inclusion of quantum-mechanical wave functions, making it challenging for studying equilibrium properties and energy relaxation processes. FSSH, in its original form, lacks proper decoherence. Recently, a novel NAMD algorithm, incorporating a density state matrix formalism, has been proposed by Kang *et al.* This modified conventional density matrix approach enables the MQC-NAMD method to account for detailed balance and decoherence effects.^{34} With this state-of-the-art NAMD approach, it becomes feasible to accurately and efficiently calculate carrier dynamics for large systems with hundreds of atoms at the first-principles density functional theory level.

Previous investigations employing the NAMD method have predominantly focused on hot carrier relaxation,^{34} transfer dynamics,^{35} and carrier cooling processes.^{36–38} However, there has been limited exploration of the carrier dynamics in the reverse process excitation from the ground state to high energy states. This aspect is particularly critical in influencing the performance of TFET and CSFET. In this study, our emphasis is on the thermalization process to gain deeper physical insights into carrier transport in Si-CSFET.^{21} The Si-CSFET incorporates an artificial type-III heterojunction of p-Si/Metal/n-Si in the source to achieve density of states engineering. The bandgap in the source region is strategically designed to cut off high-energy carriers, while a metallic silicide layer, typically composed of NiSi $ 2$ and CoSi $ 2$, is introduced to reduce the tunneling barrier, thereby enhancing the driven current. Our investigation into the carrier thermalization process in Si-CSFET utilizes the state-of-the-art NAMD method to study inelastic scatterings, such as electron–phonon scattering,^{39} in semiconductor Si and two types of metallic silicides (NiSi $ 2$ and CoSi $ 2$). We obtain key parameters, such as the thermalization factor, relaxation time, scattering time, and scattering rate, with a focus on their energy-level dependence. The relationships established in this study not only allow us to predict the thermalization time for cold carriers in semiconductors and metallic silicides but also provide valuable insights into the determinant factors influencing the thermalization process. This understanding is crucial for optimizing the design and performance of advanced semiconductor devices, such as TFETs and CSFETs.

## II. METHODOLOGY

In the adiabatic approximation, the separation of electron and nuclear motions is justified by the significant difference in their velocities; electrons move much faster than atoms. In this approximation, the much heavier nuclei are considered fixed due to their mass being three orders of magnitude greater than that of electrons. However, numerous crucial physical phenomena involve the dynamics of excited states of carriers, introducing scenarios, such as the cooling of hot carriers, including photo-generated hot carrier cooling,^{40,41} hot carrier injection in transistors,^{42} and electron–hole recombination.^{43} In such cases, where the dynamics of excited states are essential, the adiabatic approximation falls short. This is particularly true at specific temperatures where the real-time changes in atomic position due to lattice vibrations and the dynamic interactions of energy and momentum between electrons and phonons, involving electron–phonon coupling, become significant. Lattice vibrations can be effectively described by a combination of harmonic oscillators. However, the non-harmonic components of phonons in many materials necessitate the consideration of molecular dynamics simulations.^{44} To address these complexities, the NAMD method, coupled with semi-classical and semi-quantum calculation approaches, is introduced. This integration enables the calculation of electron–phonon coupling effects on carrier distribution and the time-dependent evolution process of the carrier distribution. For the cold carrier thermalization process, the energy of the cold carrier may undergo changes through subsequent interactions, such as electron–phonon scattering or electron–electron scattering, within the thermalization region. It is crucial to recognize that the distribution of cold carriers extends from a narrow low-energy range to higher energy levels, eventually reaching equilibrium over a finite time. Consideration of these intricate processes provides a comprehensive understanding of the thermalization dynamics and its impact on carrier distribution in various material systems.

In this computational framework, the behavior of electrons is determined by solving the time-dependent Schrödinger equation, while the motion of the nucleus is governed by Newton’s second law. The nuclear motion is initially calculated using Born–Oppenheimer molecular dynamics (BOMD). This method combines first principles with molecular dynamics calculations, providing a comprehensive understanding of both nuclear motion and the associated adiabatic state’s wave function and energy. By employing this approach, researchers can capture the intricate interplay between electronic and nuclear dynamics, enabling a detailed examination of the system’s behavior.

^{45}

^{,}

^{45}given by direct calculation or as an input parameter.

In this study, we employ the NAMD module within the PWmat package^{46} to explore the cold carrier thermalization process in semiconductor Si and metallic silicides (NiSi $ 2$ and CoSi $ 2$). Our calculations unfold in two distinct steps. First, we compute the Hamiltonian of the ground state electron to acquire the wave function and eigenenergy. Concurrently, the overlap of adiabatic states at various times is output within a specified nonadiabatic energy window. Additionally, the initial vibrations of each atom are randomly generated based on the temperature using $ k BT$. Subsequently, in the second step, we delve into the non-adiabatic dynamics of carriers. This entails scrutinizing the evolution process of carriers within the non-adiabatic window over time. Crucial parameters, such as the initial time, the initial occupation distribution, and the range of the non-adiabatic energy window, are thoughtfully selected to comprehensively study non-adiabatic carrier dynamics, providing valuable insights into the cold carrier thermalization process in these materials.

## III. RESULTS AND DISCUSSIONS

### A. The carrier thermalization process in bulk Si

The primitive cell of bulk Si in the diamond structure is computed using the norm-conserving pseudopotential (NCPP) in conjunction with the Heyd–Scuseria–Ernzerhof hybrid functional (HSE), and the detailed results are available in the supplementary material. Given that the NAMD module is designed to calculate carrier dynamics solely at the $\Gamma $ point, it becomes imperative to expand the unit cell to fold the Brillouin zone and generate a quasi-continuous energy-level distribution at $\Gamma $. In PWmat, the relaxed primitive cell is utilized to intercept the unit cell along the [1 1 0], [1 $ 1 \xaf$ 0], and [0 0 1] directions, as illustrated in Fig. S2(a) of the supplementary material. The supercell structure for NAMD calculations is then generated by expanding the unit cell along the a, b, and c axes by $2\xd72\xd720$ times, respectively. This results in a supercell with a volume of 6405.5 $ \xc5 3$, and its density of states (DOS) is depicted in Fig. S2 of the supplementary material.

#### 1. The carrier energy relaxation in silicon

The Si supercell comprises a total of 320 atoms, with 1280 valence electrons. The BOMD simulation is configured with the following settings: 1500 molecular dynamics steps, each with a step size of 1 fs. The total number of energy bands is set to 1600, and the exchange correlation functional is parameterized by PBE under the Generalized Gradient Approximation (GGA). The convergence of energy is achieved with a relative error less than $1\xd7 10 \u2212 8$. The results are depicted in Fig. 1(a). To minimize the tiny vibrations of the eigenlevels, a more stringent wavefunction convergence criterion of $1\xd7 10 \u2212 6$ is employed. During the BOMD calculation, 101 energy levels at and above the conduction band edge are output for subsequent NAMD calculations. Figure 1(b) illustrates the total energy and potential energy variations at 300 K, while Figs. 1(c) and 1(d) depict the atomic temperature and kinetic energy changes. The results reveal that the total energy remains constant under the convergence condition, with kinetic energy and potential energy interchanging. The initial assumption sets the kinetic energy at twice its value. After an intensive oscillation period of approximately 100 fs, the kinetic energy and the potential energy reach equipartition, resulting in a total energy around $320\xd7 3 2 k BT=12.43$ eV. Even after 200 fs, the temperature remains stable at around 300 K, indicating that the lattice vibration tends to become fairly regular. Considering these observations, the initial 200 fs, characterized by intensive vibration, is excluded in the subsequent NAMD calculations.

To explore the impact of different initial atomic vibrations, NAMD calculations are performed with an initial time range from 200 to 400 fs, with steps of 20 fs. The initial electron distribution in energy is set at the lowest energy level, indicating an occupancy of 1 on the Conduction Band Minimum (CBM) and an occupancy of 0 on other energy levels. It is important to note that carrier relaxation in this calculation does not sensitively depend on the decoherence time, and for consistency, a value of 80 fs is set.^{34} Figure 2 illustrates the change in each eigenenergy with time within the specified window at different temperatures, calculated by NAMD. The red line and the blue line represent the average energy of electrons and the Fermi level, respectively. The non-adiabatic window is defined to include 101 states above the CBM.

At low temperatures (100 K), where lattice vibrations and potential field changes are minimal, the wave function and eigenenergy levels change little over time. The overlap between energy states is small due to the rapid decay of distribution probability with increasing energy. Electron occupation near the band edge dominates, leading to a change in the average energy of electrons from initially concentrated at the band edge to slightly higher than the band edge after scattering. At 0 K, where all atoms are frozen, the potential field remains unchanged, and the wave function and energy levels are obtained through the stationary Schrödinger equation. Electrons remain in the initial energy level without scattering to other levels. As the temperature increases to 100, 200, 300, and 400 K [Figs. 2(b)–2(d)], it can be seen that the vibration of the Fermi level ( $ E F$) and eigenenergies is increased. The greater overlap between different eigenstates intensifies electron scattering between states. The average energy of electrons is higher at higher temperatures, reflecting the broader energy distribution of electrons due to reduced probability decay with energy at elevated temperatures.

To provide a clearer understanding of the thermal process at 300 K, our focus shifts to the specific occupancy probabilities of the lowest six energy levels near the band edge over time, as depicted in Fig. 3(a). Initially, electrons exclusively occupy the first energy level of the conduction band, and the occupancy of higher energy levels increases with the scattering process. Given the constant lattice vibration, each energy level’s occupancy continually oscillates. To better elucidate the overall behavior, we extract the change in the average occupation with the average intrinsic energy level within a given period, as shown in Fig. 3(b). From the initial time, electron occupation gradually transitions from lower to higher energy levels, with the change of the occupancy probability of each energy level in 100 fs. It can be seen that the fit line is different from the tendency of the Boltzmann distribution. This meanly caused by the system is under non-equilibrium conditions, and carriers may not have enough time to respond to these rapid energy fluctuations to reach the Boltzmann distribution.^{47,48} Figures 3(c) and 3(d) illustrate the coherent relationship between the intrinsic energy-level oscillation and the occupation probability variation within the 200–300 fs range. When the wave function of each energy state overlaps, intrinsic energy levels approach each other, facilitating electron jumps from lower to higher energy levels, leading to a gradual increase in the occupation probability of higher energy levels. Ultimately, a dynamic balance of scattering from low energy levels to high energy levels and vice versa is achieved. In terms of energy, the average energy of electrons increases due to lattice vibration, and the occupation of different energy levels undergoes a gradual thermalization process, involving the scattering of cold carriers to high-energy carriers. It is important to note that the BOMD calculation determines how intrinsic energy levels change, and carrier scattering depends on the specific energy-level oscillation. As different initial processes result in different scattering processes, a series of distinct initial times is set, and the final statistical thermalization process captures the system’s characteristics.

#### 2. The carrier thermalization factor and the scattering rate in Si

Initially, for energies beyond the starting energy level where there is no occupation, the thermalization factor is 0.0. As electrons gradually undergo thermalization induced by electron–phonon coupling, low-energy electrons scatter into higher energy states, causing the thermalization factor to increase. Over time, with the progression of thermalization, the thermalization factor tends to approach 1.0 after a sufficient duration. The observed oscillation arises from fluctuations in the occupation probability of high-energy levels, which can be higher or lower than the average occupation probability due to the changing relative energy between intrinsic levels with lattice vibrations. The average occupation probability is computed at thermal equilibrium after 1000 fs. As electrons can only acquire limited energy under electron–phonon interaction, the thermalization process gradually progresses from lower to higher energy levels. Thermalization at higher energy levels commences after electrons occupy their adjacent lower energy levels. The relaxation time is defined as the total time from the initial time to a certain degree of thermalization within a specific energy window. A plot of the relaxation time for different degrees of thermalization ( $\alpha $ = 0.01, 0.1, 0.5, and 1.0) as a function of the energy window, starting from the initial energy level, can be extracted, as shown in Fig. 4(b). Generally, the thermalization time increases as the energy window expands since the relaxation of higher energy levels includes the relaxation time of lower energy windows. Additionally, the error bars in Fig. 4(b) represent the standard deviation of the relaxation time.

The thermalization process toward a higher degree of thermalization (e.g., 1.0) exhibits larger fluctuations in relaxation time compared to a lower degree of thermalization. It gradually converges to thermal equilibrium through a sequence of scattering equilibria between high and low energy levels, exemplified by some of the higher energy windows in Fig. 4(a). A noteworthy phenomenon is the relatively short time at the onset of thermalization. The timescale of this initial rising edge is on the order of tens of femtoseconds. This brief duration is attributed to most lower energy levels being fully thermalized before lattice vibrations excite carriers to higher energy levels. For instance, in Fig. 4(a), when the energy window of 7.97 to 9.0 eV begins, the thermalization factor of the lower energy window is already greater than 0.1. Based on the above discussion, the physical depiction shown in Fig. 4(c) is established. Initially, cold carriers are set at the lowest energy level, and after a relaxation time period, the low-energy range is fully thermalized, while the high-energy part is not effectively thermalized. As thermalization initiates in an energy range, carriers within this range rapidly reach higher energies through faster scattering, progressively participating in the energy relaxation process. Finally, after sufficient mutual scattering between high and low energy levels, the thermal equilibrium distribution is achieved.

^{49}To emulate an atomistic structure, the device structured as 1 nm-thick Si nanowire FETs with a channel length of 10 nm keeps the same scale as the atomistic structure, as illustrated in Fig. 6(a). Under conditions of low doping concentration, the scattering rate in the lowest conduction bands in k-space remains well below $ 10 13$ [1/s], approaching an average of $2\xd7 10 12$, as shown in Fig. 6(b). Without the influence of surface roughness, scattering is primarily attributed to electron–phonon coupling, subsequently contributing to the thermalization process. The scattering time is approximately 500 fs, as converted for the result in Fig. 6(b). This aligns well with the scattering time calculated by NAMD in Fig. 5(b). Our NAMD calculation simulates carrier dynamics using a time-dependent method, which cannot be captured by TCAD software. By comparison, it is an effective approach to study other materials, which will be discussed in Sec. III B. It is evident that a higher thermalization factor corresponds to a longer scattering time. Electrons in the high-energy range portion exhibit a shorter scattering time than their relaxation time. Furthermore, akin to the relaxation time, the fluctuation of the scattering time for a higher thermalized degree ( $\alpha \u2a7e0.5$) is also more pronounced. Furthermore, the scattering rate ( $ S E$) in the energy window can be calculated by the change in the number of electrons in the window,

### B. The carrier thermalization process in metallic silicide

#### 1. Electronic structure of NiSi_{2} and CoSi_{2}

In the proposed CSFET design, a metallic material is introduced as an inserted layer to reduce the tunneling barrier compared to a regular TFET.^{10,21} The formation of a silicide/Si interface in a highly controlled manner ensures defect control and favorable electrical properties at the metal/Si junction.^{50} These advantages can be extended to Si nanowire CSFETs.^{22} The results for NiSi $ 2$^{51,52} and CoSi $ 2$^{53,54} are presented in Fig. S3 of the supplementary material, which shows similar trends. The main difference lies in the position of the Fermi level. In CoSi $ 2$, the Fermi level is lower because the Co atom lacks one $d$-orbital electron compared to Ni. For NiSi $ 2$, the minimum of DOS is below the Fermi level, and the DOS near the Fermi level increases with energy, while for CoSi $ 2$, the minimum of DOS is above the Fermi level, and the DOS near the Fermi level decreases with energy. The DOS at the $\Gamma $ point calculated by the supercell shows a discrete distribution in Figs. S4(a) and S4(c) of the supplementary material. As shown in Figs. S4(b) and S4(d) of the supplementary material, the density of states of the metal silicide itself is relatively high, and it can be better resolved after cell expansion. The total energy and temperature changes during the BOMD process of NiSi $ 2$ and CoSi $ 2$ are calculated under the supercell structure. The temperature changes with time are presented in Fig. 7. It fluctuates significantly within the initial 500 fs and then tends to 300 K under equilibrium. Similarly, to eliminate the effect of the initial violent shock on the scattering, the initial time from 500 to 980 fs is selected, and 25 groups of NAMD processes are calculated. The initial electron is set at the first energy level above the Fermi level, the occupancy rate of the remaining energy levels is 0, and the decoherence time is 80 fs.

#### 2. Carrier energy relaxation of NiSi_{2} and CoSi_{2}

In the thermalization process of NiSi $ 2$ and CoSi $ 2$ at 300 K, the eigenenergies and average energies of the lowest 61 levels above the Fermi level are shown in Figs. 8(a) and 8(b), with the initial time set at 500 fs. For NiSi $ 2$, the DOS increases above the Fermi level, leading to a progressively denser distribution of intrinsic levels above $ E F$ [Fig. S3(b) in the supplementary material]. In contrast, for CoSi $ 2$, there is a minimum above 1 eV above the Fermi level, resulting in a relatively large gap above 0.75 eV for those intrinsic levels above $ E F$ [Fig. S3(d) in the supplementary material]. Additionally, the average energy experiences a slight increase after the thermalization time. The average energy reflects the changes in occupation of the lower energy levels, whereas the occupation of higher energy levels cannot be fully captured by the average energy. This information is analyzed by examining the occupation probability of specific energy levels. The specific occupation probability of the 1st to 19th states above the Fermi level varies with time, as shown in Figs. 8(c) and 8(d). The average energy increases as the lowest energy levels reach thermal equilibrium and begin to oscillate after being thermalized for hundreds of fs (state 1 to state 7). Higher energy levels (state 10 and above) continue thermalization after 400 fs, but their equilibrium occupancy probability decays exponentially with energy. The thermalization of the higher energy component makes a limited contribution to the average energy, resulting in a relatively small increase. Consequently, it becomes challenging to reflect the thermalization effect of some energy levels through changes in the average energy. Instead, a more direct and quantitative index of the thermalization process is obtained through the occupation probability of each energy level.

#### 3. The carrier thermalization factor and the scattering rate of NiSi_{2} and CoSi_{2}

Continuingly, the thermalization factors of NiSi $ 2$ and CoSi $ 2$ in each energy window can be extracted from the changes in their occupation probabilities over time, as shown in Figs. 9(a) and 9(b). The relaxation time from the initial time to a certain degree of thermalization is presented in Figs. 9(c) and 9(d). The thermalization process of NiSi $ 2$ and CoSi $ 2$ follows a pattern similar to that of Si, where the relaxation time depends on the position of the energy window relative to the initial energy level. Specifically, the relaxation time for partial thermalization within 0.15 eV of the initial energy is within 200 fs, while the relaxation time above 0.2 eV of the initial energy exceeds 200 fs. Furthermore, the relaxation time increases with the widening of the energy window. Examining the relaxation time, it can be observed that after thermalization, the time difference between reaching thermalization factors of 0.01 and 0.1 is very small. From 0.1 to 0.5, the relaxation time difference is approximately 50 fs, and from 0.5 to 1.0, the relaxation time difference increases. This suggests that the speed of partial thermalization (thermalization factor from 0 to 0.5) is significantly higher than that of complete thermalization (thermalization factor from 0.5 to 1.0), indicating an increased probability of high-energy window scattering back to low energy levels at higher thermalization degrees. Consequently, the net scattering rate from thermalization factor 0.5 to 1.0 decreases. Meanwhile, due to the constant oscillation of energy levels, the relaxation time exhibits large fluctuations. Due to the properties of metallic silicides, the DOS within 0.5 eV near the Fermi level for NiSi $ 2$ and CoSi $ 2$ is about 0.4–0.6 eV $ \u2212 1$ [shown in Figs. 10(a) and 10(b)], which is significantly higher than that of Si [0–0.1 eV $ \u2212 1$, Fig. 10(c)]. As a result, the number of electrons at the Fermi level in the primitive cell of NiSi $ 2$ and CoSi $ 2$ (0.0089 and 0.010 68, respectively) is much higher than that in the Si conduction band ( $2.648\xd7 10 \u2212 4$) under heavy doping (assuming $ E F$ is on the CBM). For the same degree of thermalization, the number of carriers in NiSi $ 2$ and CoSi $ 2$ will be larger than that in Si. Figures 11(a) and 11(b) present the increase in the number of electrons in the energy windows of NiSi $ 2$ and CoSi $ 2$ over time. It can be found that the number of electrons in the same energy window of NiSi $ 2$ and CoSi $ 2$ is 1 to 2 orders of magnitude higher than that in Si [Fig. 5(a)]. Moreover, the scattering time of the thermalization rising edge was extracted [Figs. 11(c) and 11(d)]. It was found that the scattering time to thermalization factor of 0.5 was still in the order of sub-100 fs. By calculating the scattering rate, the curves of the scattering rate changing with the energy window are shown in Figs. 11(e) and 11(f). It can be observed that the scattering rate of metallic silicides decreases exponentially with the increase of the energy window at high energy, and the scattering rate from partial thermalization to complete thermalization is slightly lower, similar to the result in Si. However, its scattering rate is higher than that in Si, reaching 3–4 $\xd7$ $ 10 33 s \u2212 1 cm \u2212 3$ near E = 0 eV, which is nearly 2 orders of magnitude higher than that in Si (6–8 $\xd7$ $ 10 31 s \u2212 1 cm \u2212 3$). Therefore, the scattering rate in the metal may have a more significant effect on the thermalization of the cold source than that of Si.

### C. Discussion on the cold carrier transport process with thermalization

For the scattering problem under the electron–phonon coupling, the perturbation potential introduced by the lattice vibration determines the transition rate. $W( k \u2032,k)$ reflects the strength of the electron–phonon coupling. Meanwhile, the scattering rate is also related to the DOS and the occupation probability of the initial and final states. The DOS is integrated over $ k \u2032$ space, and the more initial states satisfy the transition condition, the greater the rate of scattering to $K$, so the scattering rate and the DOS are positively correlated.^{55,56} The occupation probability reflects the occupation of the initial state $ k \u2032$ and the final state $k$. The higher the occupation probability of the initial state, the lower the occupation probability of the final state, the higher the transition probability of the initial state to the final state. For the reverse process, the lower the transition probability of scattering from the $k$ to $ k \u2032$ state. In the energy relaxation process of carrier thermalization, due to the exponential decay of the occupation probability at higher energy, the scattering rate at high energy will decrease exponentially with the occupation probability.

Therefore, at the beginning of thermalization, the occupation of the initial state is high, and the occupation of the final state approaches zero. The scattering at a low degree of thermalization primarily involves transitions from the initial state at low energy to the final state at high energy, with the rate of scattering from high energy to low energy in the inverse process being very small. As the high energy level becomes partially thermalized, the scattering process from a low energy to high energy level weakens, while the inverse process scattering is enhanced. Ultimately, the scattering between high and low energy levels becomes equal, leading to the attainment of thermal equilibrium. In light of these observations, the scattering time of the high energy window in NAMD from thermalization at 0 to partial thermalization ( $\alpha =0.5$) is lower than the time from partial thermalization to full thermalization. In summary, the primary factors influencing thermalization are the DOS and the strength of electron–phonon coupling, both of which are material-dependent. The speed of thermalization can be preliminarily assessed through the scattering time, representing the pace of restoring thermal equilibrium between adjacent energy states in the system. Specifically, the scattering time for thermalization to degree 0.01/0.1 in Si is approximately 10/40 fs, while the corresponding time for degree 0.1 in NiSi $ 2$ and CoSi $ 2$ is about 30 fs. The impact of electron–phonon coupling on the scattering process appears to be negligible at low degrees of thermalization. A detailed examination of Fig. 9 reveals that the DOS for NiSi $ 2$ and CoSi $ 2$ is substantially higher than that of Si. This qualitative observation suggests that thermalization in NiSi $ 2$ and CoSi $ 2$ is more pronounced compared to Si. Furthermore, leveraging the energy relaxation scattering times obtained via NAMD for metallic silicide and silicon, it becomes possible to calculate the influence of thermalization on the transport characteristics of CSFET. This analysis would be conducted within the framework of the Boltzmann transport equation. However, this detailed investigation exceeds the scope of the current study.

## IV. CONCLUSION

In this study, the DFT-NAMD method was employed to investigate the band structure, density of states (DOS), and carrier dynamics of bulk Si, NiSi $ 2$, and CoSi $ 2$. The thermalization process of carriers under electron–phonon coupling was analyzed, highlighting the dependencies of thermalization factor, relaxation time, scattering time, and scattering rate on energy levels. Under the influence of lattice vibrations, cold carriers undergo energy-level transitions, shifting from low-energy occupation to high-energy occupation. This process leads to energy relaxation, ultimately reaching an equilibrium distribution characterized by an exponential decrease with energy. The initial thermalization time of high levels is contingent on the thermalization of low levels, resulting in a gradual increase in relaxation time with rising energy levels. The scattering rate within the energy window exhibits an exponential decrease with increasing energy levels. Consequently, the thermal leakage of the cold source channel is predominantly determined by the thermal scattering rate near the barrier. Notably, the scattering rate of NiSi $ 2$ and CoSi $ 2$ is nearly two orders of magnitude higher than that of Si, owing to the high DOS near the Fermi level. This implies that the scattering rate of metallic silicide has a more significant impact on the degradation of thermal leakage in the cold source region. In summary, for Si-based devices, metallic materials characterized by low mismatch, low DOS, and low electron–phonon coupling play a pivotal role in achieving high performance and a steep subthreshold device.

## SUPPLEMENTARY MATERIAL

See the supplementary material for more detailed computational results in this work, including the electronic properties calculated by PWmat, band structure, and density of states with Si, NiSi 2, and CoSi 2 for a primitive cell, a unit cell, and a supercell.

## ACKNOWLEDGMENTS

This work was supported by the Ministry of Science and Technology (Grant No. 2021YFA1200502), the National Natural Science Foundation of China (Grant Nos. 12174423 and 62174040), and the 13th batch of outstanding Young Scientific and Technological Talents Project in Guizhou Province [2021]5618.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Kun Luo:** Writing – original draft (equal). **Weizhuo Gan:** Investigation (equal). **Zhaozhao Hou:** Investigation (equal). **Guohui Zhan:** Investigation (equal). **Lijun Xu:** Investigation (equal). **Jiangtao Liu:** Investigation (equal). **Zhenhua Wu:** Writing – review & editing (equal).

## DATA AVAILABILITY

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

## REFERENCES

*Proceedings of 2010 International Symposium on VLSI Technology, System and Application*(IEEE, 2010), pp. 72–73.

*2018 IEEE International Electron Devices Meeting (IEDM)*(IEEE, 2018), pp. 33.2.1–33.2.4.

*2020 IEEE International Electron Devices Meeting (IEDM)*(IEEE, 2020), pp. 12–15.

*et al.*, “

*2020 IEEE International Electron Devices Meeting (IEDM)*(IEEE, 2020), pp. 22.4.1–22.4.4.

*2015 IEEE International Electron Devices Meeting (IEDM)*(IEEE, 2015), pp. 5.1.1–5.1.4.

_{2}epitaxial structures

_{2}on Si (111): The effect of a surface p-n junction

_{2}on Si (111)

_{2}/Si epitaxial structures