Atomistic insights into hydrogen migration in IGZO from machine-learning interatomic potential: linking atomic diffusion to device performance

Abstract

Understanding hydrogen diffusion is critical for improving the reliability and performance of oxide thin-film transistors (TFTs), where hydrogen plays a key role in carrier modulation and bias instability. In this work, we investigate hydrogen diffusion in amorphous IGZO (a-IGZO) and c-axis aligned crystalline IGZO (CAAC-IGZO) using machine-learning interatomic potential molecular dynamics (MLIP-MD) simulations. We construct accurate phase-specific MLIPs by fine-tuning SevenNet-0, a universal pretrained MLIP, and validate the models against a comprehensive dataset covering hydrogen-related configurations and diffusion environments. Hydrogen diffusivity is evaluated over 650–1700 K, revealing enhanced mobility above 750 K in a-IGZO due to the glassy matrix, while diffusion at lower temperatures is constrained by the rigid network. Arrhenius extrapolation of the diffusivity indicates that hydrogen in a-IGZO can reach the channel/insulator interface within 10⁴ seconds at 300–400 K, likely contributing to negative bias stress-induced device degradation. Trajectory analysis reveals that long-range diffusion in a-IGZO is enabled by a combination of hydrogen hopping and flipping mechanisms. In CAAC-IGZO, hydrogen exhibits high in-plane diffusivity but severely restricted out-of-plane transport due to a high energy barrier along the c-axis. This limited vertical diffusion in CAAC-IGZO suggests minimal impact on bias instability. This work bridges the atomic-level hydrogen transport mechanism and device-level performance in oxide TFTs by leveraging large-scale MLIP-MD simulations.

---

1. Introduction

Over the past decade, amorphous indium gallium zinc oxide (a-IGZO) and related oxides have remained at the forefront of metal oxide electronics, serving as a key channel material for thin-film transistors (TFTs) in display technologies.^1–4 This prominence is attributed to its numerous advantages, including compatibility with low-temperature processing (<200 °C), relatively high electron mobility (>10 cm² V⁻¹ s⁻¹) compared to amorphous silicon (<1 cm² V⁻¹ s⁻¹), a large band gap over 3 eV that enables low power consumption and high optical transparency, a low threshold voltage (<1 V), and a small subthreshold swing (<0.1 V dec⁻¹).^1,3–5 The success of oxide TFTs in displays has recently expanded the application of a-IGZO to emerging fields such as next-generation computing and artificial intelligence (AI), including neuromorphic transistors, offering a pathway to overcome the limitations of conventional Si-based complementary metal–oxide–semiconductor (CMOS) technology.⁶ Additionally, due to its back-end-of-line (BEOL) compatibility,⁷ the a-IGZO channel has attracted growing interest for use in monolithic three-dimensional (M3D) semiconductor devices.⁸ Indium-based oxides, including ITZO and IGO, have also exhibited high electron mobility, expanding the set of viable oxide channel materials.^9,10

Despite these advances, addressing degradation and instability in IGZO TFTs remains a critical challenge. Previous studies have identified the effects of hot carriers, oxygen vacancies, and hydrogen incorporation as key sources of instability.^11–14 Among these, hydrogen is an ubiquitous element and is known to be readily incorporated into the oxide channel during the fabrication of oxide TFTs. Two primary pathways have been suggested for hydrogen incorporation into the a-IGZO channel. First, hydrogen can be introduced during channel formation via sputtering or pulsed laser deposition (PLD), mainly due to residual hydrogen-containing species in the reaction chamber, such as H₂ O and H₂.^4,15,16 The hydrogen concentration in the channel is subsequently determined by the extent of hydrogen out-diffusion during post-deposition annealing at elevated temperatures. Second, hydrogen can also be introduced during the deposition of gate dielectric or passivation over-layers, such as SiO₂ or Al₂O₃, using plasma-enhanced chemical vapor deposition (PECVD) or atomic layer deposition (ALD).^17–19 In this case, hydrogen atoms diffuse from the overlying films into the underlying a-IGZO layer during deposition and subsequent annealing processes.

Once incorporated, hydrogen acts as a critical defect that significantly influences the electrical properties of a-IGZO and, consequently, the performance of oxide TFTs based on it. According to previous studies, the hydrogen concentration must be carefully controlled to achieve an appropriate threshold voltage (V_th) and subthreshold swing (SS).¹⁴ Particularly, excessive hydrogen incorporation can result in a highly negative V_th and a large SS.²⁰ This degradation occurs because hydrogen acts as shallow n-type donors in a-IGZO,²¹ leading to a significant increase in free electron concentration—often exceeding 10¹⁹ cm⁻³—when present in large amounts.²²

Furthermore, hydrogen impurities have been identified as a key source of electrical instabilities in oxide TFTs under gate bias and illumination stress conditions, such as negative bias stress (NBS), negative bias illumination stress (NBIS), and positive bias stress (PBS).^{14,17,23–27} In particular, under NB(I)S conditions, positively charged hydrogen species are supposed to diffuse from the channel toward the gate insulator interface. This diffusion leads to the accumulation of positive charge at the interface, resulting in a negative shift in V_th and subsequent degradation in device performance.

The foregoing discussion underscores that hydrogen diffusion in a-IGZO plays a pivotal role in determining the electrical performance and bias-related instability of oxide TFTs. Therefore, a thorough understanding of its behavior is essential for the development of high-performance and reliable oxide TFTs. Despite its significance, little is known about the detailed behavior of hydrogen diffusion in a-IGZO. In particular, hydrogen is difficult to detect using conventional experimental techniques due to its exceptionally low atomic mass.^28,29 Consequently, analyzing the hydrogen diffusion process—which requires precise, time- and space-resolved quantification of hydrogen concentration—remains highly challenging and is often accompanied by substantial uncertainty. Accordingly, only a limited number of experiments have reported quantitative analysis on hydrogen diffusion in a-IGZO or similar oxide semiconductors.^15,30,31

To address this issue, several atomistic modeling studies have investigated hydrogen diffusion in a-IGZO by means of density functional theory (DFT) calculations. These studies computed activation energies (E_a) along predefined hydrogen migration pathways²³ and performed on-lattice kinetic Monte Carlo (kMC) simulations to assess hydrogen dynamics.³² However, the reliance on predefined migration paths may not capture the full diversity of diffusion mechanisms present in a highly disordered amorphous matrix. Furthermore, on-lattice kMC methods—where transition rates between artificial lattice sites are assigned based on a set of pre-calculated activation energies—neglect spatial and energetic correlations between successive diffusion events, potentially distorting the true kinetics of hydrogen transport. To gain physically meaningful insights and enable quantitative analysis of hydrogen diffusion, more direct approaches, such as molecular dynamics (MD) simulations, are required for modeling hydrogen migration processes in amorphous systems.

First-principles molecular dynamics (FPMD) simulations, based on DFT, have been widely adopted to investigate dynamical properties such as atomic diffusion.^33,34 Despite their high accuracy, the substantial computational cost of FPMD severely limits both the system size and simulation timescale, rendering it impractical for modeling long-range hydrogen migration in disordered materials like a-IGZO. As an alternative, machine-learning interatomic potentials (MLIPs) have recently emerged as a powerful approach to overcome these limitations.³⁵ By learning energies and forces from DFT datasets, MLIPs achieve near-DFT accuracy while offering linear computational scaling ([scr O, script letter O](N)), where N is the number of atoms, making them suited for simulating large-scale amorphous systems over extended timescales.

In this work, we present a theoretical study on hydrogen diffusion in a-IGZO and c-axis aligned crystalline IGZO (CAAC-IGZO) using MLIP-MD simulations. We first develop MLIPs suitable for modeling hydrogen migration in respective phases by fine-tuning SevenNet-0 (hereafter referred to as 7net-0), a universal pretrained MLIP model.³⁶ The fine-tuned MLIP models are thoroughly validated against a diverse dataset relevant to hydrogen-related configurations and diffusion environments in IGZO. We evaluate hydrogen diffusivity at elevated temperatures ranging from 650 to 1700 K by performing MD simulations using the fine-tuned MLIP models. The results reveal that hydrogen migration in a-IGZO is facilitated at temperatures above 750 K due to the glassy nature of the amorphous matrix, resulting in a lower activation energy. In contrast, at lower temperatures, hydrogen diffusion occurs through a more rigid atomic network, leading to a higher activation energy. Detailed analysis of the MD trajectories reveals that a combination of hydrogen hopping and flipping mechanisms enables long-range hydrogen diffusion in a-IGZO. In comparison, in CAAC-IGZO, hydrogen exhibits high diffusivity along the in-plane direction, while out-of-plane diffusion along the c-axis is significantly suppressed. Extrapolation of the diffusion coefficients to 300–400 K via the Arrhenius relation suggests that hydrogen in the a-IGZO channel can reach the channel/insulator interface within 10⁴ seconds under NBS conditions, thereby accelerating device degradation. In contrast, negligible hydrogen diffusion in CAAC-IGZO along the c-axis suggests limited impact on bias-induced instability. Overall, by elucidating hydrogen migration in IGZO, our study provides key insights into its impact on the performance and electrical instability of oxide TFTs.

2. Methods

2.1. Calculation setups

As a MLIP, we employ 7net-0 (version July 11, 2024),³⁶ a pretrained universal MLIP model trained on the MPtraj dataset.³⁷ 7net-0 is built upon the NequIP architecture,³⁸ which utilizes an equivariant graph neural network framework. The model has demonstrated high accuracy across a wide range of chemical systems, including oxides, and offers excellent scalability for large-scale molecular dynamics (MD) simulations on multi-GPU platforms—a key advantage over many other graph-based MLIP approaches.

7net-0 exhibits reasonable accuracy in predicting atomic structures and dynamics of IGZO in both amorphous and crystalline phases (Fig. S1). However, in hydrogen-doped IGZO systems, we observe notable deficiencies in accuracy related to hydrogen-containing configurations (Fig. S2). To address this, we fine-tune the pretrained 7net-0 model using a curated training set that includes hydrogen-related configurations with large prediction errors. Two fine-tuned MLIP models are developed: the a-FT model for a-IGZO and the c-FT model for CAAC-IGZO. The a-FT model is obtained by fine-tuning 7net-0 on a dataset containing hydrogen configurations in amorphous structures, while the c-FT model is trained using a dataset focused on hydrogen configurations in CAAC-IGZO (Table S1). To preserve the model’s accuracy for defect-free systems, both fine-tuning datasets also include FPMD snapshots for corresponding IGZO systems without hydrogen, thereby mitigating the risk of catastrophic forgetting. Details of the fine-tuning procedure and the training set are provided in the SI. The training and validation root-mean-square errors (RMSEs) are reasonably small, comparable to previous MLIP studies,³⁹ as summarized in Table S2, illustrating the successful training of the fine-tuned models.

MLIP-MD simulations are performed using the LAMMPS molecular dynamics package⁴⁰ under the NVT ensemble using a Nosé–Hoover thermostat. For faster and memory-efficient simulations, the FlashTP CUDA kernel is employed in the convolution layers of SevenNet, resulting in approximately a threefold increase in throughput.⁴¹ A time step of 2 fs is used for systems without hydrogen, while it is reduced to 1 fs for hydrogen-containing systems to ensure MD stability. For MLIP-based structural relaxations and nudged elastic band (NEB) simulations, the atomic simulation environment (ASE) interface is employed.⁴² The convergence criterion for structural relaxation is set to 0.02 eV Å⁻¹.

DFT calculations to generate the fine-tuning and test dataset are performed using the Vienna ab initio Simulation Package (VASP),⁴³ with the Perdew–Burke–Ernzerhof (PBE) functional⁴⁴ (version 52) as the exchange–correlation functional. A plane-wave cutoff energy of 520 eV is used for a-IGZO systems, while 500 eV for CAAC-IGZO systems. For structural relaxation, the convergence criterion is set to 0.02 eV Å⁻¹. FPMD simulations are conducted under the NVT ensemble using a Nosé–Hoover thermostat. The Baldereschi k-point⁴⁵ at (0.25, 0.25, 0.25) is used for Brillouin zone sampling in the FPMD simulations. To generate the fine-tuning dataset, we perform static DFT calculations on selected snapshots from FPMD simulations, using denser k-point grids of 3 × 3 × 3 for a-IGZO and 2 × 2 × 2 for CAAC-IGZO to obtain accurate atomic forces and total energies.

2.2. Generation of a-IGZO structures

Amorphous IGZO models, with a widely used stoichiometric ratio of In : Ga : Zn : O = 1 : 1 : 1 : 4 for device fabrication,⁴⁶ are generated using a melt-quench procedure.⁴⁷ An initial atomic configuration is constructed in a cubic unit cell with a lattice parameter chosen to match the experimental density of 6.10 g cm⁻³.⁴⁸ Initially, atoms are randomly distributed while avoiding unphysically short interatomic distances. The system undergoes a short 2-ps premelting step at 5000 K for initialization, followed by a 10-ps melting step at 2500 K. To generate diverse amorphous structures, additional MD simulations are performed at 2500 K, and snapshots are extracted every 2 ps. These configurations are subsequently quenched from 2500 K to 300 K at a cooling rate of 100 K ps⁻¹. Final amorphous structures are obtained through structural relaxation. Constructing the fine-tuning and test datasets from FPMD simulations is based on a-IGZO models containing 18 formula units (In₁₈Ga₁₈Zn₁₈O₇₂), with atoms added or removed as necessary to generate relevant defect configurations.

2.3. Diffusivity calculations

The self-diffusivity of hydrogen is evaluated from its mean square displacement (MSD) during MLIP-MD simulations. In this work, we compute the time-averaged MSD over an elapsed time Δt from a given reference time t, defined as:where r^(e)_i is the position of i-th hydrogen atom in the e-th supercell. t_tot denotes the total simulation time, and N_H is the number of hydrogen atoms in the supercell. To account for statistical variations, we perform an ensemble average over N_Cell independent supercell calculations. Subsequently, we evaluate H diffusivity (D_H) from the slope of MSD-Δt curve:

When calculating eqn (2), we use the portion of the MSD curve that exhibits clear linearity with respect to Δt. Specifically, as shown in Fig. 1, we exclude the initial non-linear region (up to ∼10 Ų), which includes the ballistic regime, where vibrational motion is dominant, and hopping events confined to a small spatial region. On the other hand, as Δt increases, the time-averaged MSD becomes less reliable, causing greater variance and deviation from the expected linear MSD behavior, due to insufficient statistical sampling.³³ Therefore, we extract the diffusivity from the intermediate region of the MSD data, which reflects statistically meaningful long-range hydrogen diffusion. The total simulation time t_tot is adjusted for each system to ensure an sufficiently long linear regime beyond MSD = 10 Ų for reliable diffusivity estimation, with longer simulations conducted at lower temperatures to adequately capture hydrogen hopping events (Fig. S13 and S14). The calculated temperature-dependent diffusivities are fitted to the Arrhenius relation:

where D₀ is the pre-factor obtained from the Arrhenius fit, E_a corresponds to the apparent activation energy, and k_B is the Boltzmann constant (Tables S4–S7).

Fig. 1 Time evolution of the MSD of H atoms in a-IGZO at 650 K. The MSD behavior depends on the displacement range, and diffusivity is extracted from the linear regime corresponding to long-range diffusion.

For the amorphous phase, we use four ensemble simulation cells with the composition In₁₀₈Ga₁₀₈Zn₁₀₈O₄₃₂H₁₂, corresponding to N_H = 12 and N_Cell = 4. These structures are generated using MLIP-MD simulations via a melt-quench protocol and yield hydrogen concentrations of approximately 1.2 × 10²¹ cm⁻³, which falls within the experimentally relevant range.¹⁵ We observe that decreasing the hydrogen concentration has minimal impact on the calculated diffusivity (Table S3), suggesting that hydrogen atoms are sufficiently separated and rarely interact with each other at the concentrations considered. This suggests that the resulting hydrogen diffusivity is applicable across a range of hydrogen concentrations observed in experiments. We also confirm that variations in oxygen stoichiometry—namely, oxygen-deficient and oxygen-rich environments—within the experimentally accessible range do not significantly affect hydrogen diffusivity. In the case of CAAC-IGZO, four ensemble simulation cells with the composition In₁₁₂Ga₁₁₂Zn₁₁₂O₄₄₈H₁₂ are constructed, yielding hydrogen concentrations comparable to those in the amorphous systems. In the crystalline structures, Ga and Zn atoms are randomly distributed over their respective cation sites, consistent with experimental observations.⁴⁹

3. Results and discussion

3.1. Validation of fine-tuned models

We first validate the accuracy of the two fine-tuned MLIP models—a-FT and c-FT. Both models demonstrate improved predictability and robustness relative to DFT results, as evidenced by reduced mean absolute errors (MAEs) and higher R² values across all test sets, including pristine structures and diverse hydrogen configurations (Fig. S3–S6). Among the test cases, we closely examine systems containing an interstitial hydrogen impurity forming an O–H bond (H_i), which is the most thermodynamically favorable hydrogen defect in IGZO (Fig. 2a and b). For comparison, we also present the MAEs of the original 7net-0 model. The fine-tuned a-FT MLIP achieves significantly low MAEs—less than 2 meV per atom for energy and less than 23 meV Å⁻¹ for force. These errors are more than two times lower than those obtained with the original 7net-0 model. A comparable level of improvement is also observed for the c-FT MLIP. In addition to predictive accuracy, MD simulations based on the a-FT model accurately reproduces the distribution of metal–oxygen (M–O) bond lengths in a-IGZO as obtained from FPMD simulations as shown in the radial distribution functions (RDFs) in Fig. 2c. The In–O, Ga–O, and Zn–O bond lengths estimated from the corresponding RDF peak positions are in good agreement with previous calculations and experiments.^47,50,51 The a-FT model also produces angle distributions of metal-centered local motifs in excellent agreement with DFT results (Fig. 2d). Beyond local atomic structures, the energies and gravimetric densities of a-IGZO configurations obtained by the a-FT MLIP also show favorable agreement with DFT calculations, as illustrated in Fig. 2e and f, respectively. In Fig. 2e, the relatively long tail in the DFT distribution toward the low-energy region originates from a statistical outlier, which does not indicate a limitation of the MLIP model. The lowest-energy structure in DFT, which gives rise to this tail, is also well described by the a-FT model, with an error of only 1.5 meV per atom.

Fig. 2 Evaluation of prediction performance of fine-tuned MLIP models. (a) Energy MAEs and (b) force MAEs for systems with H_i forming a O–H bond. (c) Radial distribution function (g(r)), (d) angular distribution function (P(θ)), (e) total energy per atom, and (f) gravimetric density of pristine a-IGZO structures. In (e) and (f), white dots indicate the median and black bars denote the interquartile range.

We also confirm the improved predictive accuracy of the fine-tuned models for various types of hydrogen defects, including substitutional hydrogen occupying an oxygen site (H_O), another well-established hydrogen configuration in IGZO systems,^23,27,52 as shown in Fig. S7. Notably, the force prediction accuracy for H_O configurations—where interactions between hydrogen and neighboring cations are inadequately captured by 7net-0—is significantly enhanced through fine-tuning, as demonstrated in Fig. S8.

Next, we compare the relative energies of various H_i configurations in a-IGZO and CAAC-IGZO systems between DFT and fine-tuned MLIPs. Because hydrogen migrates between different oxygen sites, accurately capturing the energy ordering among H_i configurations is essential for reliably modeling hydrogen dynamics. As shown in Fig. 3a and b, the fine-tuned MLIPs successfully reproduce the energy hierarchy for both systems, demonstrating strong linear correlation with DFT results. We further assess the energy ordering of various hydrogen defect types in diverse atomic environments, including oxygen-deficient a-IGZO (Fig. S9), which confirms the excellent predictive accuracy of the a-FT model.

Fig. 3 Validation of fine-tuned MLIP models for relative energies of IGZO with various H_i configurations and migration barriers (E_bs). Relative energies in (a) a-IGZO (15 configurations) and (b) CAAC-IGZO (5 configurations). Migration barriers of (c) a-IGZO (24 paths) and (d) CAAC-IGZO (12 paths).

We further evaluate whether the fine-tuned MLIP models can accurately reproduce hydrogen migration barriers (E_b) in both a- and CAAC-IGZO systems. Long-range hydrogen diffusion in IGZO predominantly proceeds via a combination of hopping and flipping mechanisms, which will be discussed in detail below. To assess model accuracy for these diffusion processes, we select representative hopping and flipping events in both a- and CAAC-IGZO and calculate the corresponding migration barriers using both DFT and MLIP-based nudged elastic band (NEB) methods. As shown in Fig. 3c and d, the a-FT and c-FT models successfully reproduce the DFT-calculated migration barriers for their respective systems, demonstrating the high accuracy of the fine-tuned MLIP models in capturing hydrogen diffusion kinetics.

Of particular note is that fine-tuning provides significant benefits across diverse amorphous compositions, even though the MLIP was optimized solely for a specific composition (In : Ga : Zn : O = 1 : 1 : 1 : 4).⁵³ For example, as shown in Fig. S10 and S11, the a-FT MLIP exhibits remarkable accuracy for various compositions with In : Ga : Zn : O ratios of 1 : 1 : 2 : 5, 2 : 2 : 1 : 7, and 3 : 1 : 1 : 7, including hydrogen-containing systems, outperforming the baseline 7net-0 model. The excellent transferability of the fine-tuned MLIP can be attributed to the similar local coordination environments that persist across different a-IGZO compositions.⁵⁴

3.2. H diffusion mechanism

We perform MD simulations using the fine-tuned MLIP models to investigate hydrogen diffusion in both a-IGZO and CAAC-IGZO systems. To ensure a sufficient number of diffusion events and meaningful MSD data within accessible simulation timescales, which is essential for reliably determining hydrogen diffusivity, the simulations are conducted at elevated temperatures above 650 K.

3.2.1. Amorphous IGZO. Hydrogen diffusion in a-IGZO is investigated over a wide temperature range from 650 to 1700 K. The Arrhenius plots of hydrogen diffusivity are shown in Fig. 4. Two distinct linear regions with different activation energies are observed, intersecting near 750 K. The low-temperature trend (below 750 K) extrapolated down to 300 K yields diffusivity values that are reasonably consistent with existing experimental data,^15,55 considering the associated measurement uncertainties. It is worth noting that two previous experiments reported distinct diffusivity values, which can be partly attributed to differences in the measurement methods. Although our results are closer to those of ref. 55, it is presently reasonable to consider such experimental data as supportive benchmarks rather than definitive quantitative references. Further experimental investigations are required to determine the hydrogen diffusivity more conclusively.

Fig. 4 Arrhenius plot of H diffusivity D_H as a function of 1000/T, where T is the temperature, in a- and CAAC-IGZO. Both systems exhibit two distinct linear regions with different activation energies. ₁¹H and ₁²D indicate Hydrogen and Deuterium, respectively.

We examine hydrogen diffusion processes in a-IGZO in the low-temperature regime, based on an analysis of MD trajectories. As shown in Fig. 5a and b, the a-IGZO structure features metal–oxygen ring networks, and the introduced hydrogen atoms primarily form O–H bonds. In this situation, hydrogen diffusion occurs mainly through two processes: intra-ring hopping (IRH) and ring-to-ring flipping (RTRF). IRH refers to the hopping of a hydrogen atom between two neighboring oxygen sites within the same ring, typically between sites facing each other (Fig. 5c). Since this process is confined within a single ring, long-range hydrogen diffusion requires an additional process to exit the ring structure. RTRF enables such an escape through a large-angle reorientation of the O–H bond (Fig. 5d). Following a successful RTRF transition, the hydrogen atom can further migrate to a new O site within the newly accessed ring via IRH. Repeated sequences of RTRF followed by IRH allow hydrogen to propagate across multiple ring structures, thereby enabling long-range diffusion.

Fig. 5 Diffusion processes in a-IGZO. (a) a-IGZO model, (b) ring configurations, (c) intra-ring hopping (IRH), and (d) ring-to-ring flipping (RTRF).

We calculate the energy barriers for IRH and RTRF using NEB calculations. To account for statistical variations inherent in amorphous systems, we sample multiple diffusion pathways by MLIP NEB—159 for IRH and 121 for RTRF. For IRH, most energy barriers are found to be below 0.5 eV, with an average of approximately 0.29 eV (Fig. 6a). In contrast, RTRF exhibits higher energy barriers, with an average of 0.76 eV (Fig. 6b), making it a much slower process and the rate-determining step for long-range hydrogen diffusion. The relatively high barrier of RTRF is attributed to the temporary breaking of an M–O bond during the transition, which is subsequently reformed (Fig. 6c).

Fig. 6 Energy barriers for (a) hydrogen hopping (159 samples) and (b) hydrogen flipping (121 samples). (c) NEB images of a representative flipping process, illustrating the sequential steps of M–O bond breaking, O–H bond rotation, and M–O bond reconstruction. Contributions of high- and low-frequency transition groups for (d) hopping and (e) flipping processes to hydrogen MSD at 650 K. High and low f_hop values are 847 ns⁻¹ and 122 ns⁻¹, respectively, and High and low f_flip values are 99 ns⁻¹ and 14 ns⁻¹, respectively.

As discussed above, RTRF processes are likely to determine the overall rate of long-range hydrogen diffusion, and thus, govern its diffusivity. The critical role of RTRF is further validated by quantifying the contributions of hopping and flipping to hydrogen transport. To this end, we analyze 48 H trajectories from a 60-ns MD simulation at 650 K. For each trajectory, the total number of hopping and flipping events is counted and converted into event frequencies (f_hop for hopping and f_flip for flipping). Here, flipping is considered to occur when a H atom rotates more than 90°. Fig. 6d illustrates averaged MSD as a function of Δt for two groups of hydrogen atoms: the top 50% with the highest f_hop s, and the remaining 50%. The two groups exhibit similar MSDs, indicating a weak correlation between MSD and f_hop. This suggests that, despite occurring more frequently than flipping (with average frequencies of f_hop = 485 ns⁻¹ and f_flip = 56 ns⁻¹), hopping processes alone are not critical in determining the overall diffusivity. This is because hopping typically takes place within a single ring and often occurs in a back-and-forth manner, contributing little to net displacement. In contrast, Fig. 6e shows the MSD for two groups of hydrogen atoms classified by f_flip, and evidently, hydrogen atoms with higher f_flip s display significantly larger MSDs. This observation highlights the importance of RTRF processes in enabling long-range hydrogen diffusion.

We analyze the local structural characteristics of O atoms that influence the RTRF energy barrier and identify two key factors. First, the RTRF barrier tends to increase with the metal coordination number of O atoms (CN_OM) bonded to H (Fig. S12a). This trend can be attributed to the restricted rotational motion of hydrogen, caused by enhanced electrostatic repulsion between H⁺ and nearby cations as CN_OM increases. Second, the barrier height increases when O atoms are coordinated with more Ga or In neighbors, whereas it decreases with a greater number of Zn neighbors (Fig. S12b). Because the RTRF mechanism involves partial breaking of M–O bonds, this behavior can be rationalized by the relative bond strengths: Ga–O bonds are the strongest, followed by In–O and Zn–O, consistent with previous reports.^21,56 In particular, when combined with the effect of CN_OM, the RTRF barrier becomes especially pronounced when O atoms bonded to H are coordinated with three Ga atoms.

Thus far, we have focused on hydrogen diffusion below 750 K. In this low-temperature regime, the M–O structural framework remains relatively rigid, with minimal displacement of metal and oxygen atoms over the simulation time. In contrast, at elevated temperatures above 750 K, the increased thermal energy and high structural flexibility of amorphous phase allow for greater atomic displacements of metal and oxygen atoms (Fig. S15). Under these conditions, the instantaneous breaking of M–O bonds occurs more frequently, facilitating O–H flipping processes that are inaccessible at lower temperatures. Furthermore, without breaking the O–H bond, oxygen atoms can migrate while carrying the O–H unit through the disordered atomic network, thereby expanding the available diffusion pathways for hydrogen. These effects collectively lead to a reduction in the effective activation energy for hydrogen diffusion at high temperatures.

3.2.2. c-axis aligned crystalline IGZO. Fig. 7a shows the crystal structure of CAAC-IGZO composed of In–O and Ga/Zn–O layers alternating along the c-axis. Within the Ga/Zn–O layers, Ga and Zn randomly occupy the same crystallographic sites. Hydrogen atoms incorporated into CAAC-IGZO primarily form O–H bonds and preferentially reside within the Ga/Zn–O layers or at the boundary between Ga/Zn–O and In–O layers, whereas incorporation within the In–O layers is energetically unfavorable.^57,58 Hydrogen diffusion in CAAC-IGZO is investigated at temperatures ranging from 650 to 1700 K, and the resulting Arrhenius plot reveals two distinct linear regimes, as shown in Fig. 4. At low-temperature regime, including room temperature, CAAC-IGZO exhibits a higher diffusivity than a-IGZO. However, unlike a-IGZO, the diffusion process in CAAC-IGZO is highly anisotropic, with ab-plane (in-plane) direction diffusion significantly exceeding that along the c-axis (Fig. 8).

Fig. 7 (a) CAAC-IGZO crystal structure. Identified H migration processes: (b) Ga/Zn–O nearest neighbor hopping (GZO-NNH), (c) Ga/Zn–O c-axis flipping (GZO-F), (d) In–O ab-plane hopping (IO-H), (e) In–O ab-plane flipping (IO-F), and (f) In–O c-axis penetrating (IO-P). (g) Averaged migration barrier and standard deviation from c-FT MLIP calculations for each process. The number of samples are 112 for GZO-NNH, 112 for GZO-F, 106 for IO-H, 24 for IO-F, and 32 for IO-P.

Fig. 8 Arrhenius plot of in-plane and out-of-plane D_H as a function of 1000/T, where T is the temperature, in CAAC-IGZO. It shows highly anisotropic diffusion process of CAAC-IGZO.

To uncover the origin of the anisotropic hydrogen transport, we analyze hydrogen trajectories from MLIP-MD simulations and identify five representative diffusion processes, as illustrated in Fig. 7b–f: (i) Ga/Zn–O nearest-neighbor hopping (GZO-NNH), (ii) Ga/Zn–O c-axis flipping (GZO-F), (iii) In–O ab-plane hopping (IO-H), (iv) In–O ab-plane flipping (IO-F), and (v) In–O c-axis penetrating (IO-P). The GZO-NNH process involves hydrogen hopping between oxygen atoms coordinated to a common Ga or Zn atom, exhibiting low energy barriers averaging around 0.28 eV. IO-F processes also exhibit low migration barriers—about 0.17 eV on average—and facilitate hydrogen escape from confined hopping sites through reorientation of a O–H bond. The combination of these low-barrier processes enables long-range hydrogen diffusion along the ab-plane through a zigzag hopping pathway. On the other hand, the GZO-F and IO-H processes can also contribute to in-plane hydrogen diffusion. However, they exhibit relatively higher energy barriers—approximately 0.50 eV for GZO-F and 0.48 eV for IO-H on average—which makes their contribution weaker than those of the low-barrier GZO-NNH and IO-F processes. At low temperatures, hydrogen diffusion is particularly dominated by the GZO-NNH and IO-F pathways, whose individual energy barriers closely correspond to the low-temperature activation energy (E_a) for in-plane diffusion (Fig. 8). However, as the temperature increases, higher-barrier processes such as GZO-F and IO-H also become activated. This increases the effective activation energy and simultaneously accelerates hydrogen diffusion along the ab-plane. This is the main reason for the transition of the diffusivity curve at ∼1000 K in Fig. 4 and 8. A detailed analysis of the respective contributions from low- and high-barrier processes to D₀ and E_a would provide deeper insight into hydrogen in-plane diffusion; however, this lies beyond the scope of the present study.

Notably, we observe that hydrogen diffusion across the In–O layer along the c-axis (out-of-plane)—a pathway critical for vertical long-range transport—is significantly suppressed (Fig. 7f). NEB calculations reveal that the energy barrier for IO-P pathway, which was the only vertical hydrogen migration process observed during MD simulations, is exceptionally high, approaching 1.2 eV (Fig. 7g), closely matching E_a = 1.11 eV obtained from the Arrhenius plot of out-of-plane diffusion (Fig. 8). Such high barriers explain the negligible diffusion events along the c direction in MD simulations, particularly at low-temperature region.

3.3. Kinetic stability of M–H bonds

Although hydrogen is known to preferentially form an O–H bond in a-IGZO systems, it can also bond with metal atoms by either substituting an oxygen site (H⁺_O) or occupying an interstitial site surrounded by metal atoms (H⁻_i).^52,59 These metal-coordinated hydrogen configurations, which involve M–H bonding, may be introduced during thin-film fabrication processes. In our simulations, H⁺_i predominantly forms owing to its entropic advantage, which arises from the large number of kinetically accessible O–H configurations, in addition to its higher thermodynamic stability.⁶⁰

Nonetheless, once M–H configurations form, due to their significant kinetic stability, such configurations are suggested to persist and remain immobile at typical device operating temperatures (300–400 K).⁵² To assess their kinetic stability, we evaluate the energy barriers for hydrogen escape from a metal-coordinated H_O site to an interstitial H_i site forming an O–H bond in a-IGZO (Fig. 9a). A total of 868 MLIP NEB calculations are performed for distinct diffusion pathways, yielding a relatively high average barrier of 1.42 eV (Fig. 9b). This barrier leads to a negligible Arrhenius factor of approximately 10⁻²¹ at 350 K, indicating that such transitions are effectively suppressed under typical operating conditions. The high energy barriers are attributed to the simultaneous breaking of M–H and M–O bonds during the transition, as illustrated in Fig. 9a.

Fig. 9 (a) NEB images of a representative H_O to H_i transition, illustrating the simultaneous breaking of M–H and M–O bonds involved in the process. (b) Distribution of hydrogen migration barriers obtained from NEB calculations over 868 distinct pathways.

3.4. Effects of hydrogen diffusion on device performance

As mentioned in the Introduction, the negative shift in V_th observed in oxide TFTs under NB(I)S conditions has been attributed to the diffusion of positively charged hydrogen—present in the form of O–H bonds in the a-IGZO channel—toward the gate insulator interface. This hypothesis implicitly assumes that hydrogen is sufficiently mobile to reach the interface within the stress duration; however, this has not been conclusively demonstrated. Based on the Arrhenius relation shown in Fig. 4, the hydrogen diffusivity is estimated to be approximately D_H = 5.9 × 10⁻¹⁷ cm² s⁻¹ at 350 K. This corresponds to a diffusion length, given by , ranging from approximately 1.1 to 11 nm for time scales between 10² and 10⁴ seconds. These diffusion lengths are comparable to the typical channel thickness, approximately 20–50 nm in display devices⁴ and below 20 nm in memory devices,^61,62 suggesting that hydrogen migration is likely to contribute to the observed V_th shift in experiments involving stress durations exceeding 10² seconds. Furthermore, since hydrogen diffusion is a thermally activated process—i.e., its diffusivity increases with temperature—its impact is expected to be more pronounced at elevated temperatures. Indeed, previous experiments reported that device instability is exacerbated with increasing temperature and longer bias application times.⁶³ Although quantum tunneling effects of hydrogen were neglected in this study,⁶⁴ their inclusion would slightly enhance the predicted diffusivity without altering our conclusions regarding the role of hydrogen diffusion in NB(I)S instability. As a note, under NB(I)S conditions, hydrogen diffusion may be accelerated by the electric field. However, this effect is likely confined to the vicinity of the interface, as the external bias is predominantly applied to the gate insulator and the channel region immediately adjacent to it.⁶⁵ Therefore, most hydrogen ions in the channel, which are located more than several nanometers away from the interface, are expected to have diffusion lengths comparable to our estimates.

In planar oxide TFTs used for display applications, hydrogen diffusion along the channel length direction during the TFT fabrication also plays a critical role in determining the carrier density—and thus the V_th. According to a previous work,²⁵ micrometer-scale a-IGZO channels in top-gate TFTs exhibit lower V_th values (i.e., higher carrier densities) as the channel length decreases. This experimental observation has been attributed to more substantial and uniform hydrogen incorporation in shorter channels during the plasma-enhanced chemical vapor deposition (PECVD) of the SiO₂ interlayer dielectric (ILD) on top of them. This explanation is based on the assumption that hydrogen generated during PECVD diffuses laterally into the oxide channel from the sidewalls in contact with the source and drain electrodes, eventually reaching the channel center in short channels. Given the ILD deposition duration of 30 minutes, the estimated diffusion lengths based on our calculated diffusivity are approximately 4.6 μm and 36.5 μm at process temperatures of 260 °C and 360 °C, respectively. While the actual amount of hydrogen incorporated into the channel depends on the diffusion process through the ILD and across the ILD/channel interface, our results clearly indicate that channels shorter than 10 μm can be significantly affected by hydrogen incorporation during ILD deposition—consistent with experimental observations. Therefore, to achieve a V_th near 0 V, precise control over hydrogen incorporation during the ILD deposition is essential and should be tailored to a given channel length.

On the other hand, compared to a-IGZO, CAAC-IGZO has been experimentally reported to exhibit significantly greater resistance to electrical instabilities under NB(I)S conditions.^66,67 Notably, as in a-IGZO, hydrogen concentrations exceeding 10¹⁹ cm⁻³ can be incorporated into CAAC-IGZO thin films during channel deposition.⁶⁸ One key factor contributing to the superior bias stability of CAAC-IGZO is the suppressed hydrogen diffusion along the c-axis, as demonstrated in our simulations. When we resolve the hydrogen diffusivity by crystallographic direction, the out-of-plane (i.e., c-axis) component is estimated to be only 1.7 × 10⁻¹⁸ cm² s⁻¹ at 350 K, corresponding to a diffusion length about one order of magnitude shorter than in a-IGZO. Consequently, hydrogen in the CAAC-IGZO channel is far less likely to reach the gate insulator interface under NB(I)S conditions compared to hydrogen in a-IGZO, contributing to the enhanced reliability of CAAC-based devices.

Note that the smaller c-axis diffusivity in CAAC-IGZO compared to a-IGZO, despite its lower E_a, arises from its lower diffusion pre-factor D₀ (Tables S4 and S7). We attribute the enhanced D₀ in a-IGZO, compared with CAAC-IGZO, to its higher migration entropy. Specifically, within transition state theory (TST), D₀ is expressed as:

where g is the geometrical factor, a is the hopping distance, f is the correlation factor, is the baseline attempt frequency, and ΔS_m is the migration entropy between the initial and transition states.^69,70 In TST, the term represents the effective attempt frequency. Given the similar local coordination environments, g is unlikely to differ significantly between a-IGZO and CAAC-IGZO. The hopping distances a in both amorphous and crystalline IGZO are comparable (1–3 Å), and thus cannot account for the large difference in D₀. The correlation factor f, which reflects the non-randomness of transport, typically ranges between 0 and 1; for interstitial diffusion it is usually close to unity, since the paths are less constrained. Likewise, the baseline attempt frequency is system-independent. Accordingly, these four terms are not sufficient to explain the three orders of magnitude difference in D₀ between a-IGZO and CAAC-IGZO (along the c-axis). Instead, the entropy contribution provides a plausible explanation. In a-IGZO systems with a more flexible atomic network than their crystalline counterparts, transitions are often associated with larger entropy changes, arising from more extended collective atomic motions that require the activation of a larger number of phonon modes (Fig. S16).^71,72 A similar trend was also reported for self-diffusion in Si, where crystalline and amorphous phases exhibit comparable activation energies but differ in pre-exponential factors by five orders of magnitude.⁷³

Regarding PBS instability, manifested as a positive shift in the threshold voltage, its origin is associated with the accumulation of negative charges at the interface. In this context, the contribution of H⁺_i diffusion to PBS instability is unlikely to be dominant. Moreover, as discussed above, hydrogen impurities in M–H configurations, which can exist as negatively charged species, are not sufficiently mobile to influence device instability. Instead, other mechanisms, such as electron trapping by interface states, including those associated with H⁺_i, are likely to play a more important role in PBS instability.

4. Conclusions

In this study, we conducted a comprehensive theoretical investigation of hydrogen diffusion in a- and CAAC-IGZO using MLIP-MD simulations. By fine-tuning the universal pretrained SevenNet-0 model, we developed phase-specific MLIP models that accurately capture the structural and dynamical behaviors of hydrogen in both disordered and crystalline environments. Our results reveal that in a-IGZO, long-range hydrogen transport is facilitated by a combination of intra-ring hopping and ring-to-ring flipping mechanisms. In contrast, hydrogen diffusion in CAAC-IGZO is highly anisotropic, dominated by in-plane motion, with vertical (c-axis) diffusion strongly suppressed.

Extrapolated diffusivities at operating temperatures (300–400 K) indicate that hydrogen in a-IGZO can reach the channel/insulator interface within typical stress durations, contributing to V_th instability under NB(I)S conditions. In contrast, the suppressed vertical diffusion in CAAC-IGZO suggests a reduced impact on bias-induced degradation. Furthermore, our study highlights that hydrogen incorporation during PECVD of the SiO₂ ILD can lead to enhanced n-type doping in shorter channels, thereby influencing V_th and contributing to the observed channel-length-dependent behavior in oxide TFTs. Overall, our study offers a detailed microscopic understanding of hydrogen transport mechanisms in IGZO-based semiconductors, establishing a foundation for designing oxide TFTs with high stability and performance.

Author contributions

H. Cho: conceptualization, methodology, investigation, validation, formal analysis, data curation, and writing – original draft; M. Moon and J. Kim: investigation and methodology support; E. Koh: project coordination, communication, and resource management; H.-D. Kim: project administration and resources; R. Kim and G. Park: project support and resources; S. Han: supervision and project involvement; and Y. Kang: conceptualization, supervision, methodology, data interpretation, and writing – review & editing;

Conflicts of interest

There are no conflicts to declare.

Data availability

Training data for the MLIPs will be available upon request, except for the open database. The SevenNet code used in this study is available in the GitHub repository: https://www.github.com/MDIL-SNU/SevenNet.

Supplementary data and additional information are provided in the supplementary information (SI) of this article. Supplementary information is available. See DOI: https://doi.org/10.1039/d5tc03190g.

Acknowledgements

This research was supported by Samsung Display Company Ltd and the National Research Foundation (NRF) funded by the Korean government (MSIT) (No. RS-2024-00407840). The computations were carried out at the Korea Institute of Science and Technology Information (KISTI) National Supercomputing Center (KSC-2025-CRE-0110) and at the Center for Advanced Computations (CAC) at Korea Institute for Advanced Study (KIAS).

References

1. K. Nomura, H. Ohta, A. Takagi, T. Kamiya, M. Hirano and H. Hosono, Nature, 2004, 432, 488–492.
2. J. Y. Kwon, K. S. Son, J. S. Jung, T. S. Kim, M. K. Ryu, K. B. Park, B. W. Yoo, J. W. Kim, Y. G. Lee, K. C. Park, S. Y. Lee and J. M. Kim, IEEE Electron Device Lett., 2008, 29, 1309–1311.
3. T. Kamiya, K. Nomura and H. Hosono, Sci. Technol. Adv. Mater., 2010, 11, 044305.
4. K. Ide, K. Nomura, H. Hosono and T. Kamiya, Phys. Status Solidi A, 2019, 216, 1800372.
5. S. Samanta, U. Chand, S. Xu, K. Han, Y. Wu, C. Wang, A. Kumar, H. Velluri, Y. Li, X. Fong, A. V.-Y. Thean and X. Gong, IEEE Electron Device Lett., 2020, 41, 856–859.
6. M. V. Lopez, B. Linares-Barranco, J. Lee, H. Erfanijazi, A. Patino-Saucedo, M. Sifalakis, F. Catthoor and K. Myny, Commun. Eng., 2024, 3, 102.
7. Q. Li, C. Gu, S. Zhu, Q. Hu, W. Zhao, X. Li, R. Huang and Y. Wu, Int. Electron Devices Meet., 2022, 2.7.1–2.7.4.
8. R. An, Y. Li, J. Tang, B. Gao, Y. Du, J. Yao, Y. Li, W. Sun, H. Zhao, J. Li, Q. Qin, Q. Zhang, S. Qiu, Q. Li, Z. Li, H. Qian and H. Wu, Int. Electron Devices Meet., 2022, 18.1.1–18.1.4.
9. F. Chen, M. Zhang, Y. Wan, X. Xu, M. Wong and H.-S. Kwok, J. Semicond., 2023, 44, 091602.
10. B. K. Yap, Z. Zhang, G. S. H. Thien, K.-Y. Chan and C. Y. Tan, Appl. Surf. Sci. Adv., 2023, 16, 100423.
11. G. Zhu, M. Zhang, L. Lu, M. Wong and H.-S. Kwok, IEEE Electron Device Lett., 2024, 45, 1602–1605.
12. M. Zhang, B. Wang, Y. Zhao, Z. Jiang, L. Lu, Y. Yan, L.-C. Chen, M. Wong and H.-S. Kwok, IEEE Electron Device Lett., 2025, 46, 432–435.
13. W.-T. Chen, S.-Y. Lo, S.-C. Kao, H.-W. Zan, C.-C. Tsai, J.-H. Lin, C.-H. Fang and C.-C. Lee, IEEE Electron Device Lett., 2011, 32, 1552–1554.
14. H. J. Kim, S. Y. Park, H. Y. Jung, B. G. Son, C.-K. Lee, C.-K. Lee, J. H. Jeong, Y.-G. Mo, K. S. Son, M. K. Ryu, S. Lee and J. K. Jeong, J. Phys. D: Appl. Phys., 2013, 46, 055104.
15. K. Nomura, T. Kamiya and H. Hosono, ECS J. Solid State Sci. Technol., 2013, 2, P5–P8.
16. H. Tang, K. Ishikawa, K. Ide, H. Hiramatsu, S. Ueda, N. Ohashi, H. Kumomi, H. Hosono and T. Kamiya, J. Appl. Phys., 2015, 118, 205703.
17. T. Toda, D. Wang, J. Jiang, M. P. Hung and M. Furuta, IEEE Trans. Electron Devices, 2014, 61, 3762–3767.
18. Y. Nam, H.-O. Kim, S. H. Cho and S.-H. K. Park, RSC Adv., 2018, 8, 5622–5628.
19. H. Y. Noh, W.-G. Lee, H. G. R., J.-H. Cha, J.-S. Kim, W. S. Yun, M.-J. Lee and H.-J. Lee, Sci. Rep., 2022, 12, 19816.
20. S.-I. Oh, G. Choi, H. Hwang, W. Lu and J.-H. Jang, IEEE Trans. Electron Devices, 2013, 60, 2537–2541.
21. T. Kamiya, K. Nomura and H. Hosono, Phys. Status Solidi A, 2010, 207, 1698–1703.
22. T. Kamiya, K. Nomura and H. Hosono, J. Disp. Technol., 2009, 5, 273–288.
23. H.-K. Noh, J.-S. Park and K. J. Chang, J. Appl. Phys., 2013, 113, 063712.
24. D.-G. Kim, T.-K. Lee, K.-S. Park, Y.-G. Chang, K.-J. Han and D.-K. Choi, Appl. Phys. Lett., 2020, 116, 013502.
25. H.-C. Chen, J.-J. Chen, K.-J. Zhou, G.-F. Chen, C.-W. Kuo, Y.-S. Shih, W.-C. Su, C.-C. Yang, H.-C. Huang, C.-C. Shih, W.-C. Lai and T.-C. Chang, IEEE Trans. Electron Devices, 2020, 67, 3123–3128.
26. J. K. Jeon, J. G. Um, S. Lee and J. Jang, AIP Adv., 2017, 7, 125110.
27. A. d J. d Meux, G. Pourtois, J. Genoe and P. Heremans, Phys. Rev. Appl., 2018, 9, 054039.
28. V. Shvachko, Surf. Sci., 1998, 411, L882–L887.
29. Y.-S. Chen, C. Huang, P.-Y. Liu, H.-W. Yen, R. Niu, P. Burr, K. L. Moore, E. Martínez-Pañeda, A. Atrens and J. M. Cairney, Int. J. Hydrogen Energy, 2025, 136, 789–821.
30. Y. Qin, P. Weiser, K. Villalta, M. Stavola, W. B. Fowler, I. Biaggio and L. Boatner, J. Appl. Phys., 2018, 123, 161506.
31. V. M. Reinertsen, P. M. Weiser, Y. K. Frodason, M. E. Bathen, L. Vines and K. M. Johansen, Appl. Phys. Lett., 2020, 117, 232106.
32. G. Kim, Comput. Mater. Sci., 2022, 203, 111109.
33. X. He, Y. Zhu, A. Epstein and Y. Mo, npj Comput. Mater., 2018, 4, 18.
34. M. L. D. Franckel, M. E. Turiansky, D. Waldhör and C. G. V. D. Walle, Phys. Rev. B, 2024, 110, L220101.
35. O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko and K.-R. Müller, Chem. Rev., 2021, 121, 10142–10186.
36. Y. Park, J. Kim, S. Hwang and S. Han, J. Chem. Theory Comput., 2024, 20, 4857–4868.
37. B. Deng, P. Zhong, K. Jun, J. Riebesell, K. Han, C. J. Bartel and G. Ceder, Nat. Mach. Intell., 2023, 5, 1031–1041.
38. S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt and B. Kozinsky, Nat. Commun., 2022, 13, 2453.
39. M. Grunert, M. Großmann, J. Hänseroth, A. Flötotto, J. Oumard, J. L. Wolf, E. Runge and C. Dreßler, J. Phys. Chem. C, 2025, 129, 9662–9669.
40. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. I. T. Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
41. S. Y. Lee, H. Kim, Y. Park, D. Jeong, S. Han, Y. Park and J. W. Lee, Forty-second International Conference on Machine Learning, 2025, https://openreview.net/forum?id=wiQe95BPaB.
42. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dulak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
43. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
44. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
45. A. Baldereschi, Phys. Rev. B, 1972, 7, 5212–5215.
46. K. Nomura, T. Kamiya, H. Ohta, T. Uruga, M. Hirano and H. Hosono, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 035212.
47. Y. Kang, H. Song, H.-H. Nahm, S. H. Jeon, Y. Cho and S. Han, APL Mater., 2014, 2, 032108.
48. A. Hino, S. Kosaka, S. Morita, S. Yasuno, T. Kishi, K. Hayashi and T. Kugimiya, ECS Solid State Lett., 2012, 1, Q51–Q53.
49. N. Kimizuka, M. Isobe and M. Nakamura, J. Solid State Chem., 1995, 116, 170–178.
50. H. Hosono, K. Nomura, Y. Ogo, T. Uruga and T. Kamiya, J. Non-Cryst. Solids, 2008, 354, 2796–2800.
51. D.-Y. Cho, J. Song, K. D. Na, C. S. Hwang, J. H. Jeong, J. K. Jeong and Y.-G. Mo, Appl. Phys. Lett., 2009, 94, 112112.
52. J. Bang, S. Matsuishi and H. Hosono, Appl. Phys. Lett., 2017, 110, 232105.
53. J. Kim, J. Lee, S. Oh, Y. Park, S. Hwang, S. Han, S. Kang and Y. Kang, arXiv, 2025, preprint, arxiv:2506.15223, https://arxiv.org/abs/2506.15223.
54. S. Kang, J. Chem. Phys., 2024, 161, 244102.
55. C.-Y. Chung, B. Zhu, R. G. Greene, M. O. Thompson and D. G. Ast, Appl. Phys. Lett., 2015, 107, 183503.
56. H.-K. Noh, K. J. Chang, B. Ryu and W.-J. Lee, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 115205.
57. H. Omura, H. Kumomi, K. Nomura, T. Kamiya, M. Hirano and H. Hosono, J. Appl. Phys., 2009, 105, 093712.
58. C. Chen, K.-C. Cheng, E. Chagarov and J. Kanicki, Jpn. J. Appl. Phys., 2011, 50, 091102.
59. Y. Kang, B. D. Ahn, J. H. Song, Y. G. Mo, H. Nahm, S. Han and J. K. Jeong, Adv. Electron. Mater., 2015, 1, 1400006.
60. M. Nakashima, M. Oota, N. Ishihara, Y. Nonaka, T. Hirohashi, M. Takahashi, S. Yamazaki, T. Obonai, Y. Hosaka and J. Koezuka, J. Appl. Phys., 2014, 116, 213703.
61. A. Belmonte, H. Oh, N. Rassoul, G. Donadio, J. Mitard, H. Dekkers, R. Delhougne, S. Subhechha, A. Chasin, M. J. V. Setten, L. Kljucar, M. Mao, H. Puliyalil, M. Pak, L. Teugels, D. Tsvetanova, K. Banerjee, L. Souriau, Z. Tokei, L. Goux and G. S. Kar, Int. Electron Devices Meet., 2020, 28.2.1–28.2.4.
62. B. Tang, Z. Fang, R. Wan, S. Hooda, J. F. Leong, Q. Wan, C.-K. Chen, E. Zamburg, J.-S. Kim and A. V.-Y. Thean, Int. Electron Devices Meet., 2024, 1–4.
63. G.-W. Chang, T.-C. Chang, J.-C. Jhu, T.-M. Tsai, K.-C. Chang, Y.-E. Syu, Y.-H. Tai, F.-Y. Jian and Y.-C. Hung, IEEE Trans. Electron Devices, 2014, 61, 2119–2124.
64. E. M. McIntosh, K. T. Wikfeldt, J. Ellis, A. Michaelides and W. Allison, J. Phys. Chem. Lett., 2013, 4, 1565–1569.
65. S. Li, M. Wang, D. Zhang, H. Wang and Q. Shan, IEEE J. Electron Devices Soc., 2019, 7, 1063–1071.
66. K. Park, H.-W. Park, H. S. Shin, J. Bae, K.-S. Park, I. Kang, K.-B. Chung and J.-Y. Kwon, IEEE Trans. Electron Devices, 2015, 62, 2900–2905.
67. Y. Kang, W. Lee, J. Kim, K. Keum, S.-H. Kang, J.-W. Jo, S. K. Park and Y.-H. Kim, Mater. Res. Bull., 2021, 139, 111252.
68. T. Ono, Y. Yanagisawa, Y. Komatsu, T. Aoki, Y. Jimbo, S. Ito, Y. Yamane, N. Okuno, H. Kunitake, H. Komagata, S. Sasagawa and S. Yamazaki, Int. Electron Devices Meet., 2020, 39.4.1–39.4.4.
69. C. Wert and C. Zener, Phys. Rev., 1949, 76, 1169–1175.
70. R. J. Friauf, J. Appl. Phys., 1962, 33, 494–505.
71. V. Naundorf, M.-P. Macht, A. Bakai and N. Lazarev, J. Non-Cryst. Solids, 1999, 250, 679–683.
72. A. Yelon, B. Movaghar and R. S. Crandall, Rep. Prog. Phys., 2006, 69, 1145.
73. F. Strauß, L. Dörrer, T. Geue, J. Stahn, A. Koutsioubas, S. Mattauch and H. Schmidt, Phys. Rev. Lett., 2016, 116, 025901.

---