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.

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 3C 2/T-VTe 2/H-VTe 2/H-TaTe 2-InSe,18 metal Heusler alloy-WS 2/MoSe 2,19 KTlCl 4/Au 2S,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.

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.

Second, we take the classic path approximation, which assumes that the excited carrier has little influence on the nuclear motion within the carrier dynamic processes. The time evolution of the wave function and the overlap between the wave functions of the adiabatic states can be obtained from the BOMD, which reflects the influence of the lattice vibration on the electronic structure. On this basis, the time evolution of the carrier excited states can be obtained from the NAMD simulation of large space and time scale. At each time, the adiabatic state Φ α ( t ) and eigenenergy ε α ( t ) are determined by the Hamiltonian of single electron H ( t ), which only depends on the nuclear position R ( t ). The non-adiabatic state Ψ ( t ) can be expanded by each adiabatic states, which satisfied the time-dependent Schrödinger equation,
H ( t ) Φ α ( t ) = ε α ( t ) Φ α ( t ) ,
(1)
Ψ l ( t ) = α C α l ( t ) Φ α ( t ) .
(2)
The charge density of the collection of electrons in the system can be expressed by the density matrix as
ρ ( r ) = α , β D α , β Φ α ( r ) Φ β ( r ) ,
(3)
D α , β ( t ) = l ω l ( t ) C α l ( t ) C β l ( t ) ,
(4)
where ω l ( t ) is the statistical weight of Ψ l ( t ). A decoherence effect and a detailed balance condition can be introduced by a P-matrix. The time evolution of the density matrix satisfies the following equations:45,
t D α , β ( t ) = i k [ V α , k ( t ) D k , β ( t ) D α , k ( t ) V k , β ( t ) ] = i [ V , D ] α , β ,
(5)
with
V α , β ( t ) = δ α , β ε α ( t ) i Φ α ( t ) | Φ β ( t ) t ,
(6)
where δ α , β is the delta function. Both the decoherence effect and the detailed balance condition are contained in the NAMD45 given by direct calculation or as an input parameter.

In this study, we employ the NAMD module within the PWmat package46 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 B T. 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.

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 Γ point, it becomes imperative to expand the unit cell to fold the Brillouin zone and generate a quasi-continuous energy-level distribution at Γ. In PWmat, the relaxed primitive cell is utilized to intercept the unit cell along the [1 1 0], [1 1 ¯ 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 × 2 × 20 times, respectively. This results in a supercell with a volume of 6405.5  Å 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 × 10 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 × 10 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 × 3 2 k B T = 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.

FIG. 1.

BOMD calculation of an Si supercell with 1500 fs (a) total energy and charge density relative error, (b) the change of total energy and potential energy with time, and (c) the temperature and (d) kinetic energy change with time.

FIG. 1.

BOMD calculation of an Si supercell with 1500 fs (a) total energy and charge density relative error, (b) the change of total energy and potential energy with time, and (c) the temperature and (d) kinetic energy change with time.

Close modal

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.

FIG. 2.

Thermalization process of an Si supercell at (a) 100, (b) 200, (c) 300, and (d) 400 K in 1500 fs calculated by NAMD. The gray lines are the eigenvalues of different energy levels. The red line is the average energy E a v g of each level according to the occupation weight, and the blue line is the Fermi level E F.

FIG. 2.

Thermalization process of an Si supercell at (a) 100, (b) 200, (c) 300, and (d) 400 K in 1500 fs calculated by NAMD. The gray lines are the eigenvalues of different energy levels. The red line is the average energy E a v g of each level according to the occupation weight, and the blue line is the Fermi level E F.

Close modal

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.

FIG. 3.

The NAMD process of an Si supercell from 200 to 1500 fs at 300 K. (a) The change of occupancy over time at six energy levels near the band edge. (b) The change of an average occupancy rate at different energy levels within 100 fs. (c) The change of eigenenergies and (d) the change of occupancy of the six energy levels.

FIG. 3.

The NAMD process of an Si supercell from 200 to 1500 fs at 300 K. (a) The change of occupancy over time at six energy levels near the band edge. (b) The change of an average occupancy rate at different energy levels within 100 fs. (c) The change of eigenenergies and (d) the change of occupancy of the six energy levels.

Close modal

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

In order to quantitatively characterize the thermalization degree of the cold carrier in the thermalization process, the thermalization factor can be expressed as the following equation:
α ( t ) = ε i E w i n d o w O c c . [ ε i ( t ) ] ε i E w i n d o w O c c . [ ε i e q ] ,
(7)
where ε i is the eigenvalue of each energy level, O c c . [ ε i ( t ) ] is the occupation probability of a certain energy level at a certain time, and O c c . [ ε i e q ] is the occupation probability of the intrinsic energy level at thermal equilibrium. Figure 4 extracts the time line of the thermalization factor for each energy window from 200 to 1000 fs. In the entire energy space ranging from 7.77 to 9.0 eV, the total occupation probability is normalized, ensuring the conservation of particle number, and as a result, the thermalization factor consistently remains at 1.0.
FIG. 4.

Si supercell at 300 K. (a) Thermalization factor vs time for each energy window from 200 to 1000 fs. (b) Relaxation time extracted from the initial time to different thermalization levels, NAMD processes with an initial time of 200–400 fs, and a step of 20 fs are counted. (c) Schematic diagram of the occupation probability of cold carriers at each energy level as a function of the thermalization process.

FIG. 4.

Si supercell at 300 K. (a) Thermalization factor vs time for each energy window from 200 to 1000 fs. (b) Relaxation time extracted from the initial time to different thermalization levels, NAMD processes with an initial time of 200–400 fs, and a step of 20 fs are counted. (c) Schematic diagram of the occupation probability of cold carriers at each energy level as a function of the thermalization process.

Close modal

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 ( α = 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.

The total occupancy can be obtained by summing the occupancy of each state in the specific energy window E w i n d o w. The occupancy in the whole NAMD window is normalized to i O c c . ( ε i ( t ) ) = 1.0. The number of electrons ( N E ( t )) in E w i n d o w can be expressed as the following equation:
N E ( t ) = ε i E w i n d o w O c c . [ ε i ( t ) ] 1.0 ε 0 + g ( ε ) f ( ε ) d ε ,
(8)
where ε 0 is the eigenenergy of the lowest energy of the NAMD window. In the case of heavy doping for bulk Si, it is assumed that the Fermi level is at the band edge, and the number of electrons in each energy window above the CBM can be calculated based on the DOS of the primitive cell, as depicted in Fig. 5(a). It is evident that the number of electrons is conserved across the entire window, ranging from 7.77 to 9.0 eV. The relaxation time for higher energy ranges is longer, and the number of electrons decreases exponentially upon reaching equilibrium. This behavior aligns with the thermalization process illustrated in Figs. 4(a) and 4(b). Moreover, the extracted scattering time is presented in Fig. 5(b). The scattering time is defined as the time from the initiation of thermalization to a certain degree of thermalization within a specific energy window, representing the rising edge time of the number of electrons in that window. It is important to note that the scattering time and the relaxation time are distinct. The scattering time does not exhibit a clear dependence on the energy position. Instead, the scattering process depends on the interaction between several neighboring energy levels, rather than on the initial energy.
FIG. 5.

Bulk Si. (a) The number of electrons in different energy windows varies with time, the initial time is 200 fs. (b) The scattering time from the beginning of thermalization to different thermalized degrees. Electron scattering (c) from initial thermalization to different degrees of thermalization and (d) partial thermalization.

FIG. 5.

Bulk Si. (a) The number of electrons in different energy windows varies with time, the initial time is 200 fs. (b) The scattering time from the beginning of thermalization to different thermalized degrees. Electron scattering (c) from initial thermalization to different degrees of thermalization and (d) partial thermalization.

Close modal
Carrier scattering is computed using the Poisson–Schrödinger equation and the Fermi Golden Rule within the commercial Global TCAD Solutions (GTS) software.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 × 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 ( α 0.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,
S E = Δ N E V C Δ t ,
(9)
where V C is the volume of the cell, δ N E is the change in the number of electrons under different degrees of thermalization, and Δ t is the scattering time under different degrees of thermalization. Figure 5(c) depicts the scattering rates from the beginning to different thermalization factors ( α= 0.01/0.1/0.5/1.0), and Fig. 5(d) illustrates the scattering rates at different thermalization degrees. Some common characteristics can be discerned: the scattering rate decreases exponentially with the increase of the energy window. The high-energy part can be occupied by DOS, so the scattering rate will not be truncated. The occupation probability of the high-energy part decreases exponentially with energy because the probability of electron scattering to the high-energy window also decreases exponentially. According to the data, the S E of the entire thermalized process is slightly smaller than the scattering rate of partial thermalization. This discrepancy arises because, throughout the thermalization process from the beginning to complete thermalization, the occupation of low-energy states decreases while the occupation of high-energy states increases. The scattering rate from high-energy to low-energy states also increases, resulting in a net decrease in the scattering rate within the energy window. However, the scattering rate of lower thermalization degrees (such as 0.5 and 0.1) does not change much. This is because the occupation of low-energy states remains high, and the occupation of high-energy states remains low, leading the net scattering rate to be mainly dominated by scattering from low energy to high energy windows. Similarly, the extraction of partial thermalization scattering rates also shows the same phenomenon, as shown in Fig. 5(d). There is no distinct difference in scattering rates corresponding to lower thermalization degrees. However, for rather high thermalization degrees, the thermalization scattering rate decreases significantly as the degree of thermalization increases, e.g., from 0.5 to 1.0. Based on the extracted scattering rate, a semi-quantitative understanding of the transport process of the CSFET with thermalization and the influence of thermalization on the leakage current can be derived, which will be compared with the thermalization influence of the metal contact layer.
FIG. 6.

(a) Schematic of a device structure of 1 nm-thick Si nanowire FET and (b) the scattering rate in the conduction band in k-space.

FIG. 6.

(a) Schematic of a device structure of 1 nm-thick Si nanowire FET and (b) the scattering rate in the conduction band in k-space.

Close modal

1. Electronic structure of NiSi2 and CoSi2

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 251,52 and CoSi 253,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 Γ 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.

FIG. 7.

Supercells of (a) NiSi 2 and (b) CoSi 2 as a function of temperature (lattice vibrational kinetic energy) over 3000 fs at a 300 K.

FIG. 7.

Supercells of (a) NiSi 2 and (b) CoSi 2 as a function of temperature (lattice vibrational kinetic energy) over 3000 fs at a 300 K.

Close modal

2. Carrier energy relaxation of NiSi2 and CoSi2

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.

FIG. 8.

Thermalization process of supercells of (a) NiSi 2 and (b) CoSi 2 at 300 K in 3000 fs. Time dependence of the occupation probability of energy levels above the Fermi level for the NAMD process from 500 fs initiation to 3000 fs for the supercells of (c) NiSi 2 and (d) CoSi 2, with state 1 as the initial lowest energy level.

FIG. 8.

Thermalization process of supercells of (a) NiSi 2 and (b) CoSi 2 at 300 K in 3000 fs. Time dependence of the occupation probability of energy levels above the Fermi level for the NAMD process from 500 fs initiation to 3000 fs for the supercells of (c) NiSi 2 and (d) CoSi 2, with state 1 as the initial lowest energy level.

Close modal

3. The carrier thermalization factor and the scattering rate of NiSi2 and CoSi2

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 1 [shown in Figs. 10(a) and 10(b)], which is significantly higher than that of Si [0–0.1 eV 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 × 10 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 × 10 33 s 1 cm 3 near E = 0 eV, which is nearly 2 orders of magnitude higher than that in Si (6–8 × 10 31 s 1 cm 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.

FIG. 9.

Thermalization factor vs time at 300 K for supercells of (a) NiSi 2 and (b) CoSi 2 for each energy window. Relaxation times for supercells of (c) NiSi 2 and (d) CoSi 2 to reach different degrees of thermalization from the initial time. NAMD processes from 200 to 800 fs with a step of 20 fs are counted.

FIG. 9.

Thermalization factor vs time at 300 K for supercells of (a) NiSi 2 and (b) CoSi 2 for each energy window. Relaxation times for supercells of (c) NiSi 2 and (d) CoSi 2 to reach different degrees of thermalization from the initial time. NAMD processes from 200 to 800 fs with a step of 20 fs are counted.

Close modal
FIG. 10.

DOS near the Fermi level for (a) NiSi 2, (b) CoSi 2, (c) Si with a primitive cell, and (d) the primitive cell and the supercell structure of Si, NiSi 2, and CoSi 2.

FIG. 10.

DOS near the Fermi level for (a) NiSi 2, (b) CoSi 2, (c) Si with a primitive cell, and (d) the primitive cell and the supercell structure of Si, NiSi 2, and CoSi 2.

Close modal
FIG. 11.

Time dependence of the number of electrons in different energy windows for (a) NiSi 2 and (b) CoSi 2, scattering time from the onset of thermalization to different degrees of thermalization for (c) NiSi 2 and (d) CoSi 2, and the scattering rate of (e) NiSi 2 and (f) CoSi 2 at different thermalization degrees.

FIG. 11.

Time dependence of the number of electrons in different energy windows for (a) NiSi 2 and (b) CoSi 2, scattering time from the onset of thermalization to different degrees of thermalization for (c) NiSi 2 and (d) CoSi 2, and the scattering rate of (e) NiSi 2 and (f) CoSi 2 at different thermalization degrees.

Close modal
In this NAMD simulation, it is a dynamic evolution process with time and is set to the initial occupation distribution at the lowest energy. However, the real device will not transition from an initial carrier distribution to the carrier distribution after thermalization at the DC operating point, and the current generated by thermalization and the leakage current of the channel will be in dynamic equilibrium, as shown in Fig. 12. Taking the n-type MOSFET as an example, the cold source injects low-energy electrons, and the channel leakage of the cold source device is generated by the scattering of cold carriers to high energy in the thermalized region. When the current generated by thermalized scattering is greater than the current under the original channel thermal balance, the cold source filtering effect fails. When the current generated by thermalized scattering is less than the current under the channel thermal equilibrium, the leakage current is limited by the injection of the source and the current generated by thermalization. The cold source injection is not the limiting factor of the leakage current because the resistance of the cold source injection junction is very small and meets the requirement of high on-state current. Because the scattering rate from low energy to high energy decays exponentially with the increase of the energy window, the thermal scattering rate under the channel barrier is much higher than that near the channel barrier. The leakage current ultimately depends on the rate of scattering from the energy level near the barrier to above the barrier. The energy range above the barrier here corresponds to the energy window range calculated by NAMD, and different energy ranges represent different barrier heights. The scattering rate S E can be qualitatively explained from theory, and the scattering rate for the transition to the K-state is given by W ( k , k ) f ( k ) [ 1 f ( k ) ] W ( k , k ) f ( k [ 1 f ( k ) ] ) integrating in the k space. While W ( k , k ) is the transition rate from the k state to the k state, calculated by the Fermi Golden Rule,
W ( k , k ) = 2 π | M k , k | 2 δ ( ω k ω k ± ω q ) ,
(10)
M k , k = Ψ k | H | Ψ k ,
(11)
where ω k is the angular frequency, M k , k is the transition matrix element, H is the time-dependent perturbation potential, and ω q is the angular frequency of the perturbation potential.
FIG. 12.

(a) Schematic of an atomic structure of silicon CSFET. (b) Current transport diagram of CSFET in an off state. Taking an n-type transistor as an example, a p-type semiconductor bandgap truncation is considered an electron cold source, and the current on the channel barrier is injected through the thermal zone.

FIG. 12.

(a) Schematic of an atomic structure of silicon CSFET. (b) Current transport diagram of CSFET in an off state. Taking an n-type transistor as an example, a p-type semiconductor bandgap truncation is considered an electron cold source, and the current on the channel barrier is injected through the thermal zone.

Close modal

For the scattering problem under the electron–phonon coupling, the perturbation potential introduced by the lattice vibration determines the transition rate. W ( k , 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 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 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 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 ( α = 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.

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.

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.

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.

The authors have no conflicts to disclose.

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).

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

1.
S.
Datta
,
W.
Chakraborty
, and
M.
Radosavljevic
, “
Toward attojoule switching energy in logic transistors
,”
Science
378
,
733
740
(
2022
).
2.
A. C.
Seabaugh
and
Q.
Zhang
, “
Low-voltage tunnel transistors for beyond CMOS logic
,”
Proc. IEEE
98
,
2095
2110
(
2010
).
3.
A.
Ionescu
and
H.
Riel
, “
Tunnel field-effect transistors as energy-efficient electronic switches
,”
Nature
479
,
329
337
(
2011
).
4.
H.
Lu
and
A.
Seabaugh
, “
Tunnel field-effect transistors: State-of-the-art
,”
IEEE J. Electron Devices Soc.
2
,
44
49
(
2014
).
5.
R. J.
Prentki
,
M.
Harb
,
L.
Liu
, and
H.
Guo
, “
Nanowire transistors with bound-charge engineering
,”
Phys. Rev. Lett.
125
,
247704
(
2020
).
6.
S.
Salahuddin
and
S.
Datta
, “
Use of negative capacitance to provide voltage amplification for low power nanoscale devices
,”
Nano Lett.
8
,
405
410
(
2008
).
7.
A. K.
Yadav
,
K. X.
Nguyen
,
Z.
Hong
,
P.
García-Fernández
,
P.
Aguado-Puente
,
C. T.
Nelson
,
S.
Das
,
B.
Prasad
,
D.
Kwon
,
S.
Cheema
,
A. I.
Khan
,
C.
Hu
,
J.
Íñiguez
,
J.
Junquera
,
L.-Q.
Chen
,
D. A.
Muller
,
R.
Ramesh
, and
S.
Salahuddin
, “
Spatially resolved steady-state negative capacitance
,”
Nature
565
,
468
471
(
2019
).
8.
S.
Guo
,
R. J.
Prentki
,
K.
Jin
,
C.-l.
Chen
, and
H.
Guo
, “
Negative-capacitance FET with a cold source
,”
IEEE Trans. Electron Devices
68
,
911
918
(
2021
).
9.
K. P.
Cheung
, “On the 60 mV/dec @300 K limit for MOSFET subthreshold swing,” in Proceedings of 2010 International Symposium on VLSI Technology, System and Application (IEEE, 2010), pp. 72–73.
10.
F.
Liu
,
C.
Qiu
,
Z.
Zhang
,
L.-M.
Peng
,
J.
Wang
,
Z.
Wu
, and
H.
Guo
, “First principles simulation of energy efficient switching by source density of states engineering,” in 2018 IEEE International Electron Devices Meeting (IEDM) (IEEE, 2018), pp. 33.2.1–33.2.4.
11.
Q.
Chenguang
,
L.
Fei
,
X.
Lin
,
D.
Bing
,
X.
Mengmeng
,
S.
Jia
,
L.
Li
,
Z.
Zhiyong
,
W.
Jian
,
G.
Hong
,
P.
Hailin
, and
P.
Lian-Mao
, “
Dirac-source field-effect transistors as energy-efficient, high-performance electronic switches
,”
Science
361
,
387
392
(
2018
).
12.
M.
Liu
,
H. N.
Jaiswal
,
S.
Shahi
,
S.
Wei
,
Y.
Fu
,
C.
Chang
,
A.
Chakravarty
,
F.
Yao
, and
H.
Li
, “Monolayer MoS 2 steep-slope transistors with record-high sub-60-mV/decade current density using Dirac-source electron injection,” in 2020 IEEE International Electron Devices Meeting (IEDM) (IEEE, 2020), pp. 12–15.
13.
Z.
Tang
,
C.
Liu
,
X.
Huang
,
S.
Zeng
,
L.
Liu
,
J.
Li
,
Y.-G.
Jiang
,
D. W.
Zhang
, and
P.
Zhou
, “
A steep-slope MoS 2/graphene Dirac-source field-effect transistor with a large drive current
,”
Nano Lett.
21
,
1758
1764
(
2021
).
14.
M.
Liu
,
H. N.
Jaiswal
,
S.
Shahi
,
S.
Wei
,
Y.
Fu
,
C.
Chang
,
A.
Chakravarty
,
X.
Liu
,
C.
Yang
,
Y.
Liu
et al., “
Two-dimensional cold electron transport for steep-slope transistors
,”
ACS Nano
15
,
5762
5772
(
2021
).
15.
Q.
Wang
,
P.
Sang
,
X.
Ma
,
F.
Wang
,
W.
Wei
,
W.
Zhang
,
Y.
Li
, and
J.
Chen
, “Cold source engineering towards sub-60 mV/dec p-type field-effect-transistors (pFETs): Materials, structures, and doping optimizations,” in 2020 IEEE International Electron Devices Meeting (IEDM) (IEEE, 2020), pp. 22.4.1–22.4.4.
16.
Q.
Wang
,
P.
Sang
,
F.
Wang
,
W.
Wei
, and
J.
Chen
, “
Tunneling junction as cold source: Toward steep-slope field-effect transistors based on monolayer MoS 2
,”
IEEE Trans. Electron Devices
68
,
4758
4761
(
2021
).
17.
L. L.
Li
,
R.
Gillen
,
M.
Palummo
,
M. V.
Milošević
, and
F. M.
Peeters
, “
Strain tunable interlayer and intralayer excitons in vertically stacked MoSe 2/WSe 2 heterobilayers
,”
Appl. Phys. Lett.
123
,
033102
(
2023
).
18.
J.
Lyu
,
J.
Pei
,
Y.
Guo
,
J.
Gong
, and
H.
Li
, “
A new opportunity for 2D van der Waals heterostructures: Making steep-slope transistors
,”
Adv. Mater.
32
,
1906000
(
2020
).
19.
F.
Liu
, “
Switching at less than 60 mV/decade with a “cold” metal as the injection source
,”
Phys. Rev. Appl.
13
,
064037
(
2020
).
20.
D.
Logoteta
,
J.
Cao
,
M.
Pala
,
P.
Dollfus
,
Y.
Lee
, and
G.
Iannaccone
, “
Cold-source paradigm for steep-slope transistors based on van der Waals heterojunctions
,”
Phys. Rev. Res.
2
,
043286
(
2020
).
21.
W.
Gan
,
R. J.
Prentki
,
F.
Liu
,
J.
Bu
,
K.
Luo
,
Q.
Zhang
,
H.
Zhu
,
W.
Wang
,
T.
Ye
,
H.
Yin
,
Z.
Wu
, and
H.
Guo
, “
Design and simulation of steep-slope silicon cold source fets with effective carrier distribution model
,”
IEEE Trans. Electron Devices
67
,
2243
2248
(
2020
).
22.
W.
Gan
,
K.
Luo
,
G.
Qi
,
R. J.
Prentki
,
F.
Liu
,
J.
Huo
,
W.
Huang
,
J.
Bu
,
Q.
Zhang
,
H.
Yin
,
H.
Guo
,
Y.
Lu
, and
Z.
Wu
, “
A multiscale simulation framework for steep-slope Si nanowire cold source FET
,”
IEEE Trans. Electron Devices
68
,
5455
5461
(
2021
).
23.
G.
Qi
,
W.
Gan
,
L.
Xu
,
J.
Liu
,
Q.
Yang
,
X.
Zhu
,
J.
Zhou
,
X.
Ma
,
G.
Hu
,
T.
Chen
,
S.
Yu
,
Z.
Wu
,
H.
Yin
, and
Y.
Lu
, “
The device and circuit level benchmark of Si-based cold source FETs for future logic technology
,”
IEEE Trans. Electron Devices
69
,
3483
3489
(
2022
).
24.
A.
Svizhenko
and
M.
Anantram
, “
Role of scattering in nanotransistors
,”
IEEE Trans. Electron Devices
50
,
1459
1466
(
2003
).
25.
Z.
Stanojević
,
O.
Baumgartner
,
L.
Filipović
,
H.
Kosina
,
M.
Karner
,
C.
Kernstock
, and
P.
Prause
, “
Consistent low-field mobility modeling for advanced MOS devices
,”
Solid-State Electron.
112
,
37
45
(
2015
).
26.
J.
Yao
,
J.
Li
,
K.
Luo
,
J.
Yu
,
Q.
Zhang
,
Z.
Hou
,
J.
Gu
,
W.
Yang
,
Z.
Wu
,
H.
Yin
, and
W. M.
Wang
, “
Physical insights on quantum confinement and carrier mobility in Si, Si 0.45Ge 0.55, Ge gate-all-around NSFET for 5 nm technology node
,”
IEEE J. Electron Devices Soc.
6
,
841
848
(
2018
).
27.
Z.
Jiawei
,
L.
Bolin
, and
C.
Gang
, “
First-principles calculations of thermal, electrical, and thermoelectric transport properties of semiconductors
,”
Semicond. Sci. Technol.
31
,
1
24
(
2016
).
28.
G.
Krishnendu
and
S.
Uttam
, “
Impact ionization in β-Ga 2O 3
,”
J. Appl. Phys.
124
,
085707
(
2018
).
29.
R.
Crespo-Otero
and
M.
Barbatti
, “
Recent advances and perspectives on nonadiabatic mixed quantum–classical dynamics
,”
Chem. Rev.
118
,
7026
7068
(
2018
).
30.
J.
Ren
,
N.
Vukmirović
, and
L.-W.
Wang
, “
Nonadiabatic molecular dynamics simulation for carrier transport in a pentathiophene butyric acid monolayer
,”
Phys. Rev. B
87
,
205117
(
2013
).
31.
P. V.
Parandekar
and
J. C.
Tully
, “
Mixed quantum-classical equilibrium
,”
J. Chem. Phys.
122
,
094102
(
2005
).
32.
J. C.
Tully
, “
Molecular dynamics with electronic transitions
,”
J. Chem. Phys.
93
,
1061
1071
(
1990
).
33.
K.
Drukker
, “
Basics of surface hopping in mixed quantum/classical simulations
,”
J. Comput. Phys.
153
,
225
272
(
1999
).
34.
J.
Kang
and
L.-W.
Wang
, “
Nonadiabatic molecular dynamics with decoherence and detailed balance under a density matrix ensemble formalism
,”
Phys. Rev. B
99
,
224303
(
2019
).
35.
F.
Zheng
and
L.-W.
Wang
, “
Ultrafast hot carrier injection in Au/GaN: The role of band bending and the interface band structure
,”
J. Phys. Chem. Lett.
10
,
6174
6183
(
2019
).
36.
Z.
Vardeny
and
J.
Tauc
, “
Hot-carrier thermalization in amorphous silicon
,”
Phys. Rev. Lett.
46
,
1223
1226
(
1981
).
37.
J. M.
Richter
,
F.
Branchi
,
F.
Valduga de Almedia Camargo
,
B.
Zhao
,
R. H.
Friend
,
G.
Cerullo
, and
F.
Deschler
, “
Ultrafast carrier thermalization in lead iodide perovskite probed with two-dimensional electronic spectroscopy
,”
Nat. Commun.
8
,
1
7
(
2017
).
38.
S.
Banerjee
,
J.
Kang
,
X.
Zhang
, and
L.-W.
Wang
, “
The effects of interstitial iodine in hybrid perovskite hot carrier cooling: A non-adiabatic molecular dynamics study
,”
J. Chem. Phys.
152
,
091102
(
2020
).
39.
S.
Giannini
,
A.
Carof
, and
J.
Blumberger
, “
Crossover from hopping to band-like charge transport in an organic semiconductor model: Atomistic nonadiabatic molecular dynamics simulation
,”
J. Phys. Chem. Lett.
9
,
3116
3123
(
2018
).
40.
S.
Banerjee
,
J.
Kang
,
X.
Zhang
, and
L.-W.
Wang
, “
The effects of interstitial iodine in hybrid perovskite hot carrier cooling: A non-adiabatic molecular dynamics study
,”
J. Chem. Phys.
152
,
091102
(
2020
).
41.
M.
Bernardi
,
D.
Vigil-Fowler
,
C. S.
Ong
,
J. B.
Neaton
, and
S. G.
Louie
, “
Ab initio study of hot electrons in GaAs
,”
Proc. Natl. Acad. Sci. U.S.A.
112
,
5291
5296
(
2015
).
42.
E.
Takeda
and
N.
Suzuki
, “
An empirical model for device degradation due to hot-carrier injection
,”
IEEE Electron Device Lett.
4
,
111
113
(
1983
).
43.
W.
Chu
,
Q.
Zheng
,
O. V.
Prezhdo
,
J.
Zhao
, and
W. A.
Saidi
, “
Low-frequency lattice phonons in halide perovskites explain high defect tolerance toward electron-hole recombination
,”
Sci. Adv.
6
,
eaaw7453
(
2020
).
44.
T.
Markussen
,
M.
Palsgaard
,
D.
Stradi
,
T.
Gunst
,
M.
Brandbyge
, and
K.
Stokbro
, “
Electron-phonon scattering from Green’s function transport combined with molecular dynamics: Applications to mobility predictions
,”
Phys. Rev. B
95
,
245210
(
2017
).
45.
K. F.
Wong
and
P. J.
Rossky
, “
Solvent-induced electronic decoherence: Configuration dependent dissipative evolution for solvated electron systems
,”
J. Chem. Phys.
116
,
8429
8438
(
2002
).
46.
W.
Jia
,
Z.
Cao
,
L.
Wang
,
J.
Fu
,
X.
Chi
,
W.
Gao
, and
L.-W.
Wang
, “
The analysis of a plane wave pseudopotential density functional theory code on a GPU machine
,”
Comput. Phys. Commun.
184
,
9
18
(
2013
).
47.
J.
Ren
,
N.
Vukmirović
, and
L.-W.
Wang
, “
Nonadiabatic molecular dynamics simulation for carrier transport in a pentathiophene butyric acid monolayer
,”
Phys. Rev. B
87
,
205117
(
2013
).
48.
Q.
Gao
and
J.
Kang
, “
Hot carrier relaxation in CsPbBr 3 nanocrystals: Electron–hole asymmetry and shape effects
,”
Phys. Chem. Chem. Phys.
24
,
9891
9896
(
2022
).
49.
Z.
Stanojevic
,
O.
Baumgartner
,
F.
Mitterbauer
,
H.
Demel
,
C.
Kernstock
,
M.
Karner
,
V.
Eyert
,
A.
France-Lanord
,
P.
Saxe
,
C.
Freeman
, and
E.
Wimmer
, “Physical modeling—A new paradigm in device simulation,” in 2015 IEEE International Electron Devices Meeting (IEDM) (IEEE, 2015), pp. 5.1.1–5.1.4.
50.
K.-N.
Tu
,
R.
Thompson
, and
B.
Tsaur
, “
Low Schottky barrier of rare-earth silicide on n-Si
,”
Appl. Phys. Lett.
38
,
626
628
(
1981
).
51.
R.
Tung
,
J.
Gibson
, and
J.
Poate
, “
Formation of ultrathin single-crystal silicide films on Si: Surface and interfacial stabilization of Si-NiSi2 epitaxial structures
,”
Phys. Rev. Lett.
50
,
429
(
1983
).
52.
R.
Tung
,
K.
Ng
,
J.
Gibson
, and
A.
Levi
, “
Schottky-barrier heights of single-crystal NiSi2 on Si (111): The effect of a surface p-n junction
,”
Phys. Rev. B
33
,
7077
(
1986
).
53.
R.
Tung
,
J.
Bean
,
J.
Gibson
,
J.
Poate
, and
D.
Jacobson
, “
Growth of single-crystal CoSi2 on Si (111)
,”
Appl. Phys. Lett.
40
,
684
686
(
1982
).
54.
E.
Rosencher
,
S.
Delage
,
Y.
Campidelli
, and
F. A.
d’Avitaya
, “
Transistor effect in monolithic Si/CoSi2/Si epitaxial structures
,”
Electron. Lett.
20
,
762
764
(
1984
).
55.
M. V.
Fischetti
and
S. E.
Laux
, “
Monte Carlo analysis of electron transport in small semiconductor devices including band-structure and space-charge effects
,”
Phys. Rev. B
38
,
9721
(
1988
).
56.
S.
Poncé
,
E. R.
Margine
,
C.
Verdi
, and
F.
Giustino
, “
EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions
,”
Comput. Phys. Commun.
209
,
116
133
(
2016
).