Research Article  www.advelectronicmat.de Device-to-Materials Pathway for Electron Traps Detection in Amorphous GeSe-Based Selectors Amine Slassi, Linda-Sheila Medondjio, Andrea Padovani, Francesco Tavanti, Xu He, Sergiu Clima, Daniele Garbin, Ben Kaczer, Luca Larcher, Pablo Ordejón, and Arrigo Calzolari* 1. Introduction The choice of the ideal material employed in selector devices is a tough task both from the theoretical and experimental side, especially due to the lack of Chalcogenides are the materials of choice for nonvolatile phase-change memory a synergistic approach between techniques able to correlate specific material (PCM)[1–3] and ovonic threshold switching properties with device characteristics. Using a material-to-device multiscale (OTS)[4,5] selector devices.[6,7] While for technique, a reliable protocol for an efficient characterization of the active the case of PCM memories good can- traps in amorphous GeSe chalcogenide is proposed. The resulting trap maps didate materials have been identified trace back the specific features of materials responsible for the measured mainly within the family of germanium- antimony-tellurium (GST) compounds findings, and connect them to an atomistic description of the sample. The (with different stoichiometry[8,9] and metrological approach can be straightforwardly extended to other materials doping[10,11]), the ideal material for OTS and devices, which is very beneficial for an efficient material-device codesign has still to be identified notwithstanding and the optimization of novel technologies. the many candidates proposed so far, such as GeSe, GeS, AsTe, AlTe, TeAsSiGe, SiGeInAsTe, GeTeSbS, and TeAsSiGeP.[12] Particularly relevant is the case of amor- A. Slassi, F. Tavanti, A. Calzolari phous GeSe alloys[5,13–15]: While a lot of effort has been devoted CNR-NANO Istituto Nanoscienze to GeTe as a prototypical phase change material, and to GeSe2 Centro S3, Modena I-41125, Italy as an efficient optoelectronic system, much less is known E-mail: arrigo.calzolari@nano.cnr.it about the quasi-stoichiometric GeSe compounds. Experimental L.-S. Medondjio, X. He, P. Ordejón results indicate that GeSe exhibits ovonic electrical switching Catalan Institute of Nanoscience and Nanotechnology (ICN2) in the amorphous phase[13,16] and that the modulation of the CSIC and BIST Campus Bellaterra, Barcelona 8193, Spain stoichiometry Ge/Se ratio affects the electrical response of the A. Padovani[+] device, in terms of I–V characteristics, threshold potential, University of Modena and Reggio Emilia switching response, endurance, and power dissipation.[5] Selec- DIEF tors based on Se-rich GeSe materials exhibit reduced leakage Via P. Vivarelli, 10, Modena, MO 41125, Italy current, improved thermal stability but higher threshold voltage S. Clima, D. Garbin, B. Kaczer with respect to the Ge-rich counterparts[17] that, at odd, show imec Kapeldreef 75, Leuven B-3001, Belgium superior operation endurance. Despite the huge effort spent on L. Larcher this subject, there is no ubiquitous consensus on the origin of AMAT – Applied Materials Italia this behavior and on the role of the Ge concentration on the Via Meuccio Ruini 74L, Reggio Emilia, Italy characteristics of selectors. On general grounds, the research of the optimal elec- The ORCID identification number(s) for the author(s) of this article tronic material requires a deep understanding of its micro- can be found under https://doi.org/10.1002/aelm.202201224. scopic properties (especially atomic defects, composition, and doping) and how they determine the electrical characteristics © 2023 The Authors. Advanced Electronic Materials published by of the final device. While device characteristics such as power Wiley-VCH GmbH. This is an open access article under the terms of the consumption, endurance, switching time, or variability are Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly certainly correlated to the properties of the constituent mate- cited. rials, univocal structure–property relations have not yet been unraveled. The origin of this mismatch relies on the fact that [+] Present address: Formerly Applied Materials Italia, Via Meuccio Ruini experimental measurements, such as I–V and C–V charac- 74L, Reggio Emilia, Italy teristics, are indirect and extensive, which stem from a set of conditions (e.g., temperature, bias, composition, doping) DOI: 10.1002/aelm.202201224 that average and/or compensate in the final results. For these Adv. Electron. Mater. 2023, 9, 2201224 2201224 (1 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de reasons, a coherent protocol able to treat at the same level of 2. Device-to-Materials Approach accuracy both the materials and the device would be of great value for the comprehension and the further design of the The characterization of the material parameters (e.g., compo- final device. sition, thickness, morphology, stoichiometry) that mostly affect Prior works exploited density functional theory (DFT) cal- the device response is a tremendous task, especially when it culations and quasi-Fermi level models[18,19] to evaluate the deals with the structural and electronic properties of amor- electronic structure of the materials, as well as simplified phous systems. Even though the computational materials dis- analytical models to estimate the material atomic structure, covery is a major task for the design of new technology, this glass transition temperature, and subthreshold current.[20] top-down approach is often infeasible due to economic con- Although providing useful insights, these approaches do straints of the sample’s fabrication process. Thus, an alterna- not give a comprehensive and self-consistent picture of tive solution is to follow the inverse device-to-materials pathway the device electrical characteristics, since they are not able that extracts the intrinsic properties of the materials from the to represent the material impact on the device electrical per- experimental data. Since the electrical characteristics of chal- formances, i.e., the response of the material in the specific cogenide selectors and their temperature dependence are well device configuration. Indeed, GeSe structure, composition, explained by a trap-assisted carrier model, the identification of and trap impact on the device characteristics of OTS selec- defects/traps becomes crucial in the process. tors remain elusive. Here, we developed the three-step approach, summarized In this work, we propose a multiscale/multiphysics approach in Figure 1: Step 1: Experimental setting involves the sample that combines experiments, device modeling, and atomistic preparation, the device fabrication, and the measurement of material simulations to identify the traps in amorphous mate- the electrical characteristics (I–V, G–V). Step 2: The experi- rials affecting the observed subthreshold conduction of selec- mental output is used as the input of a specific defect spectros- tors. In particular, we adopt a device-to-materials approach that copy (DS) functionality, implemented in the device-oriented extracts traps/defects characteristics from the experimental modeler Ginestra(R),[22] which is specifically developed to extract device electrical data and connects them to the microscopic intrinsic material properties and characteristics of the defects properties of the material, in a learning process that correlates (thermal ionization and relaxation energies, distribution within the electrical properties of the specific device to the electronic the bandgap) from the electrical measurements. Step 3: The properties of the constituent material. While a few attempts of resulting defect mapping is interpreted in terms of ab initio using this device-to-materials approach have been proposed in simulations, based on density functional theory. the case of crystalline systems,[21] similar approaches have not The core of the process is the DS functionality that relies been used to unveil the role of defect/traps in amorphous mate- on the accurate modeling and simulation of the carrier-defect rials, where the lack of long-range order opens to conceptual interactions that are at the bases of the defect-assisted trapping problems about the structural models and the defect definition, and transport (trap-assisted tunneling, TAT),[23] as described in especially for off-stoichiometric systems that do not have a crys- the framework of the multiphonon theory which accounts for talline counterpart. carriers–phonons coupling and lattice relaxation processes[23–25] This methodology is applied to two-terminal metal/insulator/ (see the Experimental Section). The DS approach comprises metal (MIM) selectors, based on amorphous GexSe100−x films four main steps: i) The measured device/structure is recre- with different thickness and stoichiometry, where x = {50, 60} ated within the Ginestra(R) simulation platform, starting from is the percentage of Ge in the sample. While Se-rich GeSe com- the known device and material parameters (layers thickness, pounds have been largely investigated, even in comparison with bandgap, dielectric constant, etc.). ii) The measured experi- the GeSe2 crystal, here we focus on the case of Ge-rich systems mental data are loaded into the DS tool, as a function of the (e.g., Ge60Se40) for which a crystalline counterpart does not temperature, voltage, and frequency (the latter for G–V data exist. Our results self-consistently explain the experimental cur- only). In order to guarantee high accuracy and uniqueness of rent at different temperatures, the conductance data at different the results, the DS methodology has to be applied to multiple frequencies, the changes of material properties as a function of electrical data (e.g., I–V and G–V as in this paper) measured the Ge content and the atomic defects dominating trap-assisted under multiple conditions (voltages, temperatures, and fre- transport. The microscopic analysis of electron trap distribution quencies). iii) Definition of the material and trap parameters to is eventually carried out in terms of short and medium range be extracted as well as their variation range to be considered order specific of the selected stoichiometry, without invoking during the optimization procedure. This information is defined ill-defined comparisons with defects (e.g., vacancies) in crystal- in the input file. iv) An automatic optimization procedure accu- line systems. rately reproduces the experimental data by varying the selected Figure 1. Sketch of the device-to-materials workflow. Adv. Electron. Mater. 2023, 9, 2201224 2201224 (2 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de parameters within the specified intervals by minimizing the GeSe has a small bandgap from ≈0.8 to ≈1.5 eV depending overall error between the experimental and the simulated on the stoichiometry, which is indicative of a small Schottky- electrical characteristics. A modified root mean square error like energy barrier at the TiN/GeSe interface (barrier larger is used for this purpose. At the end of the procedure, the DS than the bandgap would correspond to an Ohmic-like contact). module provides the energy-versus-position (E, z) map of the The resulting electrical characteristics for devices D21 and D22 defects responsible for the observed voltage, temperature, and are shown in Figure 2a,b and c,d, respectively. An excellent frequency dependence of the experimental data (uniform spa- agreement between experiments and simulations is obtained tial distributions in the whole GeSe thickness are considered for both stacks in the whole range of voltages and tempera- in this work for simplicity). It also provides selected material/ tures considered. The consistency of the applied model and of device parameters (e.g., work function and energy bandgap) as corresponding results has been checked against the thickness well as the defects relaxation energy (EREL), which represents change, by comparing the data for devices D22 and D16 (see the energy required for the lattice relaxation process to occur Figures S1 and S2, Supporting Information). when the defect state changes its charge state.[23–25] EREL is esti- Starting from the I–V and G–V characteristics, we extracted mated from the temperature dependence and determines the the material and trap parameters, as reported in Table 2. The charge capture/emission time constants and their temperature Ge-rich system has a lower bandgap (1.1 eV) than the stoichi- dependence and it is linked to the defects’ atomic structure. ometric one (1.45 eV), in agreement with previous works[26,27] Therefore, it can be used, together with the extracted thermal and consistent with the reduction of the experimental gap ionization energy band, to identify the nature of the defects observed in the D21 device. The DS defect maps reported in responsible for the observed electrical data.[23] Figure 3 show for both compositions the presence of two trap bands: one for holes and one for electrons. In both cases, trap distributions remain spatially uniform through the thickness of 3. Results the sample. The trap density (NT,i) is similar for electrons and holes, and it is higher for the Ge-rich case than stoichiometric 3.1. Step 1: Device Fabrication and Experimental Data Generation one. However, important differences can be detected in the two Three different GexSe100−x OTS layers were prepared and inte- samples. In the case of Ge50Se50 trap bands are slightly nar- grated in two-terminal MIM selectors, by combining materials rower and located on average at ≈0.4 and ≈0.2 eV from the with different thicknesses (10 and 20 nm) and Ge concentra- valence and conduction band edge, respectively. The Ge60Se40 tion (x = 50 and 60) as reported in Table 1. The amorphous compound is characterized by broader trap bands located on GexSe100−x films were prepared by room temperature physical average at ≈0.1 and ≈0.3 eV from valence band top and the con- vapor deposition (PVD). The layers were integrated into TiN/ duction band bottom. At zero bias traps are either fully occu- GeSe/TiN two-terminal structures in a 300 mm process flow pied or empty, with the Fermi level lying in the middle of the and passivated with a low-temperature back end of line process. gap. Once the external bias is applied, both carrier types con- The details of fabrication and measurement sets are reported tribute to conduction, with a net prevalence of electron conduc- in Section S1 (Supporting Information). Each current–voltage tion through the upper trap band (n-type). However, conduc- (I–V) or conductance–voltage (G–V) measurement was per- tion type may depend on structural parameters of the device, formed on a fresh device, and multiple (at least three) same- such as thickness, effective bandgap, electrode work-function. sized devices were measured at each measurement condition. Indeed, the analysis for the D16 sample exhibits an equal con- tribution of both hole and electron carriers. In order to investigate the role of traps in the transport 3.2. Step 2: Device Modeling and Trap Mapping mechanisms within the GeSe layer, we recalculated the I–V characteristics of the D22 sample, by selectively deactivating The proposed DS technique is applied to the I–V and G–V specific scattering terms (Figure 4a). In particular, we consid- data measured on the GeSe samples D21 and D22, described ered the electrical characteristics resulting i) without the inclu- above; results for the D16 system are reported in Section S1 sion of the defect contributions, ii) without the trap-to-trap tran- (Supporting Information). Since the conductance is found to be sitions, and iii) without the trap-to-band transitions. The results independent of the frequency, only the data measured at 2 kHz at different temperatures for the D22 sample are summarized are considered for the DS procedure. in Figure 4b–d. The absence of trap contributions [panel (b)] imparts a reduction of the simulated current (dashed lines) by one order of magnitude over the entire voltage range [−2 V, Table 1. MIM capacitors based on GexSe100−x with varying Ge content +2 V], and a remarkable worsening of the agreement with the x and varying chalcogenide thickness. Data about D16 sample are experimental data (circles), only partially compensated by the reported in the Supporting Information. increase of the temperature. The suppression of the trap-to-trap [panel (c)] or of the trap-to-band [panel (d)] contributions has Sample ID M/I/M capacitor smaller effects on the I–V characteristics: trap-to-trap seems D16 TiN/10 nm Ge60Se40/TiN more relevant at low bias (V < 1.0 eV), while trap-to-band has a major effect for higher bias (V > 1.5 eV). We can conclude D21 TiN/20 nm Ge50Se50/TiN that the most relevant contribution for current stems from D22 TiN/20 nm Ge60Se40/TiN the presence of the trap bands within the bandgap, while the Adv. Electron. Mater. 2023, 9, 2201224 2201224 (3 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de Figure 2. a,c) I–V and b,d) C–V characteristics of 20-nm-thick TiN/Ge50Se50/TiN (D21, top row) and TiN/Ge60Se40/TiN (D22, bottom row) capacitors, at different temperatures. Dots (lines) correspond to measured (simulated) data. other mechanisms are responsible for the specific details of of the leakage current, as shown in Section S1, Supporting the resulting electrical characteristics. This can be justified in Information). However, an accurate extraction of this key para- terms of the simultaneous small bandgap and high trap density meter is possible in principle only if the observed current–tem- of both systems that make trap-to-trap and trap-to-band ener- perature dependence originates from purely trap-assisted trans- getically similar. port. In the case of the thick, small bandgap GeSe layers con- Another important parameter that can be extracted through sidered here, this could be not the case and other mechanisms, the DS tool is the relaxation energy (EREL), which is directly con- such as tunneling from/to the electrodes may play a role: nected to the atomic nature of the trap and is typically derived simulations show that the charge transport occurs by means from the temperature dependence of the charge transport (i.e., of carrier excitation at the metal TiN/GeSe interface and then through the defect band by means of trap-to-trap transitions. Since both processes (excitation and trap-to-trap) are strongly Table 2. Material and trap parameters extracted by using the DS utility to reproduce I–V and G–V data measured, as a function of tempera- dependent on the temperature, the EREL value of 0.1–0.5 eV ture, on the different OTS devices. EG is the (mobility) bandgap, tGeSe for both electrons and holes that is consistently extracted for is the film thickness, ET,i and ΔET,i are average and spread of the traps the different GeSe samples has to be verified also with other thermal ionization energy, EREL,i is the traps relaxation energy, and NT,i methods (e.g., ab initio calculations). is the traps density for electrons (i = el) and holes (i = h), respectively. The extracted trap density is consistent among samples and within the range typically reported in the literature. Besides Parameter Ge50Se50 Ge60Se40 characterizing the properties and distribution of GeSe traps, the D21 D22 DS tool also allows the extraction of the GeSe effective thick- EG 1.45 eV 1.10 eV ness (tGeSe). As reported in Table 2, the extracted value of tGeSe is ET,el 0.23 eV 0.34 eV very consistent for both stoichiometries and in very good agree- ΔET,el 0.15 eV 0.2 eV ment with the nominal thickness. EREL,el 0.15 eV 0.4 eV This analysis demonstrates the DS capability to characterize the trap distribution of samples, and to discern among dif- NT,el 5× 1018 cm−3 2.48 × 1019 cm−3 ferent compositions, temperatures, and geometrical properties ET,h 1.1 eV 0.9 eV of the device. However, the DS tool is not sufficient to unravel ΔET,h 0.15 eV 0.2 eV the structural properties and chemical origin of those traps. EREL,h 0.15 eV 0.3 eV For crystalline systems, the analysis of defect states is routinely NT,h 5 × 1018 cm−3 2.48 × 1019 cm−3 done by considering the effects of point defects (e.g., vacancies, interstitial, impurities) that usually cause the appearance of tGeSe 21.3 nm 20.5 nm localized states in the semiconductor band gap. In the present Adv. Electron. Mater. 2023, 9, 2201224 2201224 (4 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de Figure 3. Defect band extracted by applying DS on the data in Figure 2 for (a) D21 and (b) D22 device, at V = −1.0 V. Colored dots represent point defects, the side color scale indicates the current driven by each trap. Horizontal gray lines mark the Fermi level of the two external leads (i.e., the applied bias), the red and blue lines across the sample are the calculated quasi-Fermi level for holes and electrons, respectively. case, the amorphous structure of GeSe prevents from using structure without resorting to experimental information or this intuitive approach and demands for an extensive micro- adjusted interatomic potentials, we carried out an extensive scopic investigation of realistic disordered structures. set of independent simulations performed with different mul- tiscale techniques and different codes, including molecular dynamics (MD) using interactions from both classical potentials 3.3. Step 3: First Principles Calculations (LAMMPS code[28]) and density functional theory (SIESTA[29,30] and QuantumEspresso[31]). In order to interpret the DS results from a microscopic point of As a starting point, we considered the atomic model of view, we performed a high throughput study of the structural Ge50Se50 resulting from classical MD simulations on a large and electronic properties of GeSe compounds based on first simulation cell including 4480 atoms (labeled big cell). The principles calculations. Given the intrinsic disordered nature structure is obtained through nanosecond-long melt-and- of the amorphous materials, the understanding of the inter- quench cycles (5 K ps−1 cooling rate), following the procedure play between local-order structures and electronic trap states described in ref. [32]. This approach allows the system to requires the capability of treating large systems at the atom- rearrange at the local- and medium-order structures avoiding istic level and of having access to the corresponding electronic unphysical “frozen-liquid” structures, and it has been demon- structure. In the absence of direct experimental results about strated to generate a realistic model for the amorphous GeSe the atomic structure of the glassy phase at room temperature, compound.[32,33] The increase of Ge content from 50% to 60% the identification of a reliable geometry is a major challenge for in GexSe100−x causes a systematic change of the average coordi- theory. In order to obtain realistic models of the amorphous nation number of Ge and Se species that moves from 3:3 (as Figure 4. a) Graphical sketch of relevant transport mechanisms assumed in the DS modeler. Simulation of I–V characteristics of D22 sample at different temperatures, by selecting specific transport mechanism: b) without traps, c) without trap-to-trap transitions, and d) without trap-to-band transitions. Dots (lines) correspond to measured (simulated) data. Adv. Electron. Mater. 2023, 9, 2201224 2201224 (5 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de in the crystal) to a phase where Ge has a mixed fourfold (as in of side L = 21.0 Å. Starting from these initial structures, the Ge-bulk) and threefold coordination. The different short-range amorphous models have been obtained through melting- structures and their spatial arrangement in the medium range and-quenching cycles, carried on by using ab initio MD. The are expected to have a critical role in the formation/appearance quenching step is particularly critical to obtain trustable of the electrical trap states (see below). We assumed the big cell amorphous results. Preliminary results (see Section S2, Sup- structure as the reference testbed for the further comparisons. porting Information) indicate that i) the relaxation of volume We first exploited the capability of the SIESTA code to treat under constant pressure, ii) the use of a relatively large basis large systems by simulating one 4480 atoms structure extracted set (in particular, double-ζ plus polarization [DZP][34]) during from the MD simulation at 300 K (i.e., one frame from MD tra- the quenching procedure, and iii) a low quenching rate are jectory). The structure was fully relaxed using the DFT forces necessary steps for the building of realistic GeSe amorphous computed in SIESTA, through a total-energy-and-force opti- structures. The details of the quenching-and-melting approach mization approach to allow the system to readjust the atomic are summarized in Table 3. Further computational details and distances and angles on the basis of the DFT potential. The accuracy tests are reported in the Experimental Section and in comparison between the original MD and the DFT results is the Supporting Information. summarized in Figure 5a–c, which reports the partial radial In the second approach (labeled cut), we extracted from the distribution function g(r) of the (a) GeGe, (b) GeSe, and (c) original MD trajectories on the large (4480 atoms) GexSe100−x SeSe pair distances. Except for the expected minor differences (x = 50, 60) models, 12 smaller cubic boxes containing due to bond-length optimization, the optimized DFT results 500 atoms each, with the same stoichiometries and densities (labeled DFT big cell) well reproduce the classical ones (labeled of the parent MD systems. These smaller cells (cut from the MD big cell). larger box and imposing periodic boundary conditions) were However, a meaningful analysis of an amorphous system relaxed using classical MD at a constant temperature of 300 K requires the investigation of several frames to gain a statistical for 10 ns, in order to allow the broken bonds at the borders to average of the structural properties. This makes the use of such rearrange and match the new periodic boundary conditions. a large cell practically unfeasible for extensive DFT simulations. Each system has been further optimized at the DFT level, by This opens the problem on how to generate smaller atomic using the QuantumEspresso code[31] (see the Experimental Sec- models that are simultaneously structurally reliable and numer- tion and Supporting Information for further details). ically affordable. We generated two independent sets of struc- The averaged radial distribution functions of model struc- tures, by exploiting two different approaches and DFT codes. tures resulting from the m&q and cut approaches are shown in In the first one (labeled m&q) we performed a sequence of Figure 5d–f. A direct comparison clearly indicates that the final melting-and-quenching cycles on model structures, by using configurations resulting from the two approaches resemble the SIESTA code. Initial Markov chain Monte Carlo simulations the original big one, in particular the radial distribution func- have been used to generate the disordered starting structures tions are very similar to the original results (Figure 5, top of GexSe100−x (x = 50, 60), with N = 336 atoms in a cubic box panels). This confirms that we generated atomically different Figure 5. Radial distribution functions g(r) for a,d) GeGe, b,e) GeSe, c,f) SeSe pair distances of Ge50Se50, simulated through different techniques. Results reported in the top layer are relative to the 4480-atom system (big cell); results in the bottom layer correspond to smaller cells (m&q and cut). Adv. Electron. Mater. 2023, 9, 2201224 2201224 (6 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de Table 3. Melt-and-quench (m&q) protocol for generating the amor- region, from −10 to −6 eV (peak B), is dominated by the Ge(4s) phous structures. SZ and DZP indicate single-ζ basis and double- contribution, while the region from −5 to 0 eV (peak A) has a polarized-ζ basis sets; NPT indicates isothermal–isobaric ensemble. mixed Se(4p)-Ge(4p) character, with a major contribution from Se states, which varies with the Ge/Se ratio. The lowest energy Initial structure Generate initial structure with MCMC (336 atoms). Then relax part of the conduction band has a mixed GeSe behavior as structure (cell and atom coordinates) well, but here the major contributions derive from the Ge Melt 6 ps at 1200 K, SZ, NPT + 6 ps at 1200 K NPT atoms, as a function of the composition. This characteristic Quench 1 1200–800 K, 3 ps in total, 1 fs per step, peak distribution is typical of GeSe compounds, such as crys- τ relax = 1000 fs: SZ, NPT talline GeSe and GeSe2.[36–38] Despite these overall common Quench 2 800–400 K, 3 ps in total, 1 fs per step, aspects, the gap regions, which are the most relevant for the τ relax = 750 fs; SZ, NPT transport properties, exhibit different characteristics depending Quench 3 400 K, 3 ps in total, 1 fs per step, on the stoichiometry. τ relax = 750 fs; DZP, NPT Figure 6 shows the DOS plots resulting from the ensemble average on the 12 structures of Ge50Se50 and Ge60Se40, respectively, as resulting from the cut approach. A Hubbard- but statistically equivalent systems, whose structural character- like potential has been included to both Ge and Se atoms, in istics are independent from the procedure and the numerical order to correct the severe underestimation of the bandgap details used to realize the model. Similar considerations hold due to standard DFT approach (see Section S2, Supporting for the corresponding electronic structures, even though in this Information, for further details). The zero-energy reference case the specific choice of the numerical details (e.g., exchange- is set to the corresponding Fermi level. Red diamonds indi- correlation functional) may affect the quantitative value of cate the inverse participation ratio (IPR) of single particle parameters such as the bandgap. However, as evident from the orbitals,[39,40] which is a measure of the spatial localization of comparison of the electronic density of states (DOS) reported in the electronic states: the higher the IPR the more localized Figure S8 (Supporting Information), both DFT approaches give the states. the same qualitative results. This assures that the trap states The mobility gap Eg is determined as the energy-level analysis is robust and independent from the specific choice of border separating the localized and extended bands, i.e., as the used DFT code and functional. the energy region where the IPR increases with respect to the On the basis of these results, we performed a set of cal- deep valence/conduction bands. Increasing the Ge amount culations to generate model structures for the two different in the material, the mobility gap Eg decreases from 1.4 eV stoichiometries corresponding to the experimental samples (x = 50%) to 0.7 eV (x = 60%) in agreement with the experi- D21 and D22, i.e., Ge50Se50 and Ge60Se40. The Se-rich coun- mental data and the DS results. The density, the localization terpart Ge40Se60 has also been simulated for comparison (see and the energy distribution of defect states in the gap strongly the Supporting Information). For all GeSe compounds, sev- depend on the stoichiometry. The nature of the mobility-gap eral (12) atomic structures have been generated to have a sta- states is found to be different from the Valence-Alternating- tistical meaningful insight of the amorphous characteristics. Pairs (VAP) model for Se-rich GeSe.[41–43] In the case of Ge50Se50 The complete set of DOS plots for each simulated structure is [panel (a)], the gap is almost empty with a few defect states that summarized in Figures S9–S11 (Supporting Information). The accumulate close to the edges of the bandgap, forming two trap electronic properties of the three stoichiometries have similar bands (dashed cyan lines) for holes and for electrons, in agree- spectral features: the valence DOS plots are characterized by ment with the DS results of Figure 3. The states forming these three distinctive peaks (labeled A, B, and C in Figure S12, Sup- defect bands are spatially localized (Figure 7a,b) but involve dif- porting Information) in agreement with recent X-ray spectros- ferent local coordination depending on n-type or p-type char- copy.[35] Peak C corresponds to the energy region going from acter of the trap. A similar picture also holds for the Se-rich −14 to −12 eV and has mostly a Se(4s) character. The central Ge40Se60 compound (Figure S13, Supporting Information), Figure 6. Averaged density of states (DOS) plot for amorphous GeSe compounds at different stoichiometries. Red diamonds indicate the inverse participation ratio (IPR) of single-particle states. Vertical gray (cyan) dashed lines mark the mobility gap Eg (trap bands) of the system. The zero-energy reference is set to the Fermi level of each system. Adv. Electron. Mater. 2023, 9, 2201224 2201224 (7 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de valence band (peak A, Figure S12, Supporting Information) has a prevalent Se character, for all stoichiometries (including Ge60Se40) hole defects have a mixed GeSe character, which cannot be simply ascribed to the presence or the absence of a single element, as for the case of point defects in crystals. While the concept of vacancy or species-deficiency may heuris- tically be extended to the case of the stoichiometric compound (x = 50%) in view of the comparison with the crystalline coun- terpart, the same concepts are completely ill defined in the case of off-stoichiometric systems. Indeed, the 20% of compositional unbalance as in Ge40Se60 or Ge60Se40 compounds causes an overall modification of the interatomic bonding (bond atoms, bond length, bond angles) and not simply a local perturbation next neighbor to the defect. From this analysis we conclude that: 1) defect states have an interbond spatial distribution of >1 nm; 2) they always have a partial Se contribution; and 3) Ge coordination plays a cru- cial role in the formation of defect states. Hole and electron defects have different characteristics: the former do not exhibit the presence of homopolar GeGe bonds that are, instead, a Figure 7. Single-particle charge-density plots of selected defect states in the mobility gap of a,b) Ge50Se50 and c,d) Ge60Se40 for occupied (left) and common feature of the latter. All electron defect states involve empty (right) states. Gray polyhedrons indicate ordered threefold and Gehomopolar bonds. However, not all the GeGe bonds fourfold local-order structures centered on Ge atoms. Ge (Se) atoms are give rise to a state in the gap. The question is: what is the colored in purple (green), respectively. Red (blue) isosurfaces indicate the origin of these states? The answer stems from the structural positive (negative) sign of the single-particle wavefunction. analysis of the medium range order of the system, and in par- ticular from the formation of GeSe complexes with different even though the defect states are more distributed across the coordination. bandgap, giving rise to broader trap bands. More complex is To illustrate this point, we considered the trap states relative the case of Ge60Se40 (Figure 6b): even though two trap bands to a single snapshot of Ge50Se50, assumed as a testbed of the are still visible, the defect states are largely spread within the system. Similar arguments can be extended to the other GeSe energy gap. This agrees with the appearance of the broader cases. We focused on one occupied (Th) and one empty (Tel) trap defect bands in Figure 3b. However, in this case the in-gap state, which decorate the band edges of the system (Figure 8), states are spatially delocalized all over the cell (Figure 7c,d), so and we analyzed the coordination of the atoms closely involved they cannot properly be considered as trap states, but rather as in these trap states. The occupied defect state (Th) is localized conductive states. This goes in parallel with the enhancement around three close Ge atoms that are not forming homopolar of the current measured in sample D22 (Figure 2). bonds, but that are simultaneously bonded to Se atoms having Figure 7 shows the single-particle charge density plots of a different coordination, namely, threefold (III) and fourfold (IV) few selected defect states of the GeSe systems for both occupied Se atoms. This Se-coordination mix results in the formation (a, c) and empty (b, d) states. Gray polyhedrons indicate ordered of an unbalanced and localized charge distribution around the threefold and fourfold local-order structures centered on Ge Ge atoms that causes the appearance of the trap state in the atoms, that are among the most frequent local arrangements gap. The Tel state is localized across a triple GeGeGe chain, found in these systems.[33] From the analysis of the spatial dis- where one Ge (the one on the right) is bonded to three inequiv- tribution of the in-gap states, we draw a few important obser- alent Se atoms, having II, III, and IV coordination, respectively. vations: i) at odd with the usual crystalline case, the localized Again, the trap state originates from the off-coordination of the states of these amorphous structures are not localized around bonded Se atoms, and not from the Ge itself or from the for- the single atomic site (e.g., the point defect), but involve a few mation of the Gehomopolar bonds. To confirm this statement bonds. ii) The defect states often appear at interconnections we analyzed an extended state (labeled No trap in Figure 8), among a few threefold and fourfold polyhedrons, indicating where we identified the formation of a GeGe bond. However, a key role of Ge coordination on the formation of trap states, both Ge are bonded to homo-coordinated Se atoms (SeIII). This in agreement with previous theory reports.[38] iii) In the stoi- prevents the localization of charge and the formation of a trap chiometric case [panels (a) and (b)], the abundance of Se favors state, in agreement with extended machine learning investiga- the formation of GeSe bonds arranged in ordered local struc- tions on GeSe alloys.[44] We can conclude that trap states origi- tures (local order). This justifies the presence of the reduced nate from local electron under/over-coordination, in agreement number of trap states in the Ge50Se50 system. On the contrary, with previous general statements.[18,38,44] Yet, our results indi- the deficiency of Se in the Ge-rich case causes the formation cate that off-coordination is not associated to a single atomic of homopolar GeGe bonds, which tend to be interconnected site as usually suggested, but rather to nm-extended complexes across the cell instead of organizing in folded structures (i.e., (medium range) that include both chemical species. The mul- higher delocalization). Note that the number of polyhedrons tielement/multibond complexity of the trap origin justifies the decreases as the Ge content increases. iv) Since the top of the different interpretations appeared in the literature, often based Adv. Electron. Mater. 2023, 9, 2201224 2201224 (8 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de Figure 8. a) Density of states and IPR of a representative snapshot of the Ge50Se50 system. Labels mark the energy positions of the selected states of panels (b)–(d). b,c) The atomic coordination pattern of the hole (Th) and electron (Tel) trap states, respectively. d) An extended state (No trap) is included for comparison. on simplified structural models, such as homopolar bonds or challenge within the solid-state community, and goes beyond folding of single chemical elements (e.g., Ge-fold structures). the aim of this work. Since trap-assisted transport implies the hopping of charge to/from localized states and the 0/−1 (0/+) charge-state transi- tion for electron (hole) carriers, we investigated the effect of 4. Conclusion including an extra charge to the system. For each initial snap- shot of the Ge50Se50 system, we simulated the case of an extra We presented a multiscale device-to-materials approach for hole or an extra electron in the unit cell. We fully relaxed the the identification and the analysis of the trap distribution charged structures and we averaged the resulting electronic in electronic devices. Starting from the experimental elec- structures, as summarized in Figure 9. trical characteristic of two-terminal TiN/a-GeSe/TiN selector Even though the two characteristic trap bands of Figure 3 devices, we fully investigated the properties of the active are still visible, the inclusion of the extra charge in the system material at the microscopic level and its transport response. imparted a redistribution of the defect energy within the This approach allowed us to connect the performances and bandgap that results more occupied with respect to the neu- reliability of the experimental devices to the composition and tral case. We observe also an increase of the average IPR of the the stoichiometry of the system. In particular, we unravel the defect states, which corresponds to an increase of the spatial role of charge traps in amorphous chalcogenides and the rela- charge localization favored by the local structural relaxation tion of these defect states with local order of the system. The (i.e., polaron-like effect). The direct evaluation of the charge- proposed D2M protocol is general and it can be applied to state transition energies for amorphous systems is an open any CMOS, logic or memory device, independently from the Figure 9. Averaged DOS and IPR plot for Ge50Se50 in the inclusion of a) an extra hole (Q = +1) and b) an extra electron (Q = −1) in the simulation cell, where Q is the total charge of the system. Vertical dashed line marks the Fermi level of the system, assumed as the zero-energy reference. Adv. Electron. Mater. 2023, 9, 2201224 2201224 (9 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de structure, the composition, and the chemical species of the again to 400 K in 3 ps. At this temperature, the structure was equilibrated active materials. additionally for 3 ps using now the DZP and a τrelax of 750 ps. Our learning approach combines characterization and The relaxation and AIMD calculations were carried out with the SIESTA code.[29,30] SIESTA employs local basis sets, and norm-conserving metrology aspects that may accelerate the optimization of mate- pseudopotentials from the pseudo-dojo dataset,[52] in which the Se and rial-design capabilities by exploring new materials and new Ge 4s and 4p electrons are considered as valence states. The generalized compounds, starting from the desired electrical performances gradient approximation (GGA)-Perdew–Burke–Ernzerhof (PBE)[53] of the device. functional was used to describe the exchange-correlation. DRSLL van der Waals correction was used in all calculations.[54] Single-ζ (SZ) basis functions[34] were used for the initial relaxation and finite temperature molecular-dynamics simulation, except in the last AIMD run in which 5. Experimental Section DZP basis functions were used. The basis set was generated with Fabrication and Characterization: Amorphous GexSe100−x films an orbital-confining scheme with a cutoff radius corresponding to were prepared by room-temperature PVD. TiN/GeSe/TiN MIM the energy shift of 0.001 Ry.[34] Only the Γ-point was used for the capacitors were integrated in a 300 mm CMOS process flow. Electrical Brillouin zone integration. The convergence criterion for the geometry characterization was performed with a Keithley 2536 SMU and an optimization was set to a maximum force of 0.04 eV Å−1 on each atom. Agilent E4980A LCR meter. For further details see Section S1 (Supporting Extended accuracy test on the choice of the basis set is reported in Information). Section S2 (Supporting Information). Device Modeling: Device simulations were performed with the Ginestra Cell Cutting (Cut) Approach: For each stoichiometry, 12 atomic software package (Ginestra), a multiscale simulation platform that self- structures of GexSe100−x each including 500 atoms were generated. consistently describes all the mechanisms relevant for logic,[23,45,46] and The initial geometries were extracted from the large (4480 atoms) memory devices.[47,48] Non-defect-related conduction mechanisms such structures generated with the classical interatomic potentials, discussed as direct/Fowler–Nordheim tunneling, thermionic emission, and drift above. Each structure was then optimized by using a total energy and were considered consistently with defect-assisted charge trapping and forces approach based on DFT calculations, as implemented in the transport (TAT) model that controls the conduction in a wide variety of QuantumEspresso package.[31] The electron–ion interactions were materials.[23,45,49,50] TAT is described in the framework of the multiphonon treated within the projector augmented wave (PAW) framework. Single theory accounting for carriers–phonons coupling and lattice relaxation particle wave functions and charge were expanded in planewaves up to processes, which depends on the atomistic structure of the defect and energy cutoff of 30 and 300 Ry, respectively. The exchange-correlation is described by the relaxation energy parameter (EREL). Hydrodynamic functional was described by the vdW-DF2-b86r exchange-correlation transport, current-induced power dissipation and temperature increase, functional[55] that takes into account the van der Waals contributions distortion and breakage of atomic bonds promoted by field, temperature which are known to play a crucial role in amorphous chalcogenides.[56] and electron injection,[46] diffusion of atomic species,[47] ferroelectric All the atomic positions were allowed to move until the maximum ionic switching,[48] and phase change were also self-consistently considered, forces were less than 0.002 eV Å−1. Integrations over the Brillouin zone thus enabling the modeling of material modifications associated to were sampled at the Γ-point only. In order to overcome the well-known the electrical stress and specific material properties. Calculations were bandgap underestimation of DFT, a DFT+U correction was included performed considering the local potential given by the applied bias and to the electronic structure, within the pseudohybrid Hubbard density the defect charge state and occupation. Such a comprehensive and self- functional approach (namely ACBN0[57]) at single point scf-calculations. consistent description of carriers transport and trapping mechanisms The resulting Hubbard-like values are: UGe = 0.16 eV and USe = 2.28 eV. is at the basis of the DS technique[21] that was developed for extracting The detailed analysis of the effect of the inclusion of the Hubbard-U defects and properties of material used in modern and novel logic and potential in comparison with the hybrid HSE06 approach is reported in memory devices. the Supporting Information. Materials Simulations: Classical MD Simulations: The Ge50Se50 system was simulated by using a classical molecular dynamics approach, according to Tavanti et al.[32] The system contains 4480 atoms and the force field used is based on the well-established Vashishta potential Supporting Information employed for the description at the atomistic level of GeSe2 or InP Supporting Information is available from the Wiley Online Library or systems.[51] The classical MD simulations were carried out by using from the author. the LAMMPS package.[28] The complete description of the adopted model and the analysis of the classical MD results were reported elsewhere.[32,33] Melting-and-Quenching Approach: Various random supercells with Acknowledgements Ge:Se coordination number 4:2 and 3:3 were constructed by a Markov Chain Monte Carlo method. In the generation of the structures, a A.S. and L.-S.M. contributed equally to this work. This work was repulsive interaction between atoms was included to make their supported in part by EC through H2020-NMBP-TO-IND project GA no. distribution more homogeneous in the cell, and potential which 814487 (INTERSECT) and H2020- NMBP-TO-IND-2018-2020 project GA depends on the Ge and Se coordination, to impose a Ge:Se coordination no. 953167 (OpenModel). ICN2 was supported by the CERCA Program of 4:2. The resulting geometry was then relaxed prior to the melting- of Generalitat de Catalunya, and by grant PGC2018- 096955-B-C43 and-quenching procedure by using a conjugate gradient algorithm. funded by Spain’s MCIN/AEI/10.13039/501100011033 and by ERDF “A Ab initio molecular dynamics (AIMD) were then carried out within the way of making Europe”. L.-S.M. received funding from the European isobaric–isothermal (NPT) ensemble at zero external pressure with an Union’s Horizon 2020 research and innovation programme under annealing algorithm, by rescaling the velocities and the pressure at each the Marie Skłodowska-Curie Grant Agreement No. 754558 (PREBIST time step of Δt = 1 fs and the characteristic relaxation time of τrelax. The – COFUND). L.-S.M.’s work has been performed in the context of the simulations started by equilibrating each initial configuration at 1200 K Physics Ph.D. program of the Universitat Autonòma de Barcelona. for 12 ps, at which the system melts and loses the memory of the initial Computer resources were provided in part by the TACC supercomputing configuration. After equilibration, each system was quenched to 800 K center (TX, USA) through QUANTRANS project, and by the Spanish over a period of 3 ps; this corresponds to an average cooling rate of 133 K Red Nacional de Supercomputación through Grant Nos. FI-2020-2- ps−1. During the melting and the quenching to 800 K procedure, the τrelax 0037, FI-2020-3-0034, FI-2021-1-0016, FI-2021-2-0017, FI-2021-3-0031, FI- was set to 1000 fs. The obtained configuration at 800 K was quenched 2022-1-0035, and FI-2022-2-0042 at MareNostrum 4, BSC-CNS. Adv. Electron. Mater. 2023, 9, 2201224 2201224 (10 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de Conflict of Interest [21] A. Padovani, B. Kaczer, M. Pesic, A. Belmonte, M. Popovici, L. Nyns, D. Linten, V. v. Afanas’ev, I. Shlyakhov, Y. Lee, H. Park, The authors declare no conflict of interest. L. Larcher, IEEE Trans. Electron. Device 2019, 66, 1892. [22] Applied Materials Ginestra, https://www.appliedmaterials.com/ products/applied-mdlx-ginestrasimulation-software. [23] L. Vandelli, A. Padovani, L. Larcher, R. G. Southwick, W. B. Knowlton, Data Availability Statement G. Bersuker, IEEE Trans. Electron. Devices 2011, 58, 2878. The data that support the findings of this study are available from the [24] K. Huang, A. Rhys, Proc. R. Soc. London 1950, 204, 406. corresponding author upon reasonable request. [25] M. Zhang, Z. Huo, Z. Yu, J. Liu, M. Liu, J. Appl. Phys. 2011, 110, 114108. [26] D. I. Bletskan, J. J. Madyar, V. N. Kabaciy, Semiconductors 2006, 40, 137. Keywords [27] Z. Jahangirli, Turk. J. Phys. 2011. devices, DFT, electrical characteristics, Ginestra, materials, ovonic [28] S. Plimpton, J. Comput. Phys. 1995, 117, 1. threshold switching, selectors, trap defects [29] J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, D. Sánchez-Portal, J. Phys.: Condens. Matter 2002, 14, 2745. Received: November 10, 2022 [30] A. García, N. Papior, A. Akhtar, E. Artacho, V. Blum, E. Bosoni, Revised: January 7, 2023 P. Brandimarte, M. Brandbyge, J. I. Cerdá, F. Corsetti, R. Cuadrado, Published online: February 7, 2023 V. Dikan, J. Ferrer, J. Gale, P. García-Fernández, V. M. García-Suárez, S. García, G. Huhs, S. Illera, R. Korytár, P. Koval, I. Lebedeva, L. Lin, P. López-Tarifa, S. G. Mayo, S. Mohr, P. Ordejón, A. Postnikov, Y. Pouillon, M. Pruneda, et al., J. Chem. Phys. 2020, 152, 204108. [1] S. Raoux, W. Wełnic, D. Ielmini, Chem. Rev. 2010, 110, 240. [31] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, [2] H.-S. P. Wong, S. Raoux, S. Kim, J. Liang, J. P. Reifenberg, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari B. Rajendran, M. Asheghi, K. E. Goodson, Proc. IEEE 2010, 98, 2201. Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, S. Baroni, [3] D. Loke, T. H. Lee, W. J. Wang, L. P. Shi, R. Zhao, Y. C. Yeo, J. Chem. Phys. 2020, 152, 154105. T. C. Chong, S. R. Elliott, Science 2012, 336, 1566. [32] F. Tavanti, B. Dianat, A. Catellani, A. Calzolari, ACS Appl. Electron. [4] P. Noé, A. Verdy, F. d’Acapito, J.-B. Dory, M. Bernard, G. Navarro, Mater. 2020, 2, 2961. J.-B. Jager, J. Gaudin, J.-Y. Raty, Sci. Adv. 2020, 6, 1. [33] F. Tavanti, A. Calzolari, ACS Omega 2022, 7, 23255. [5] S.-D. Kim, H.-W. Ahn, S. y. Shin, D. S. Jeong, S. H. Son, H. Lee, [34] E. Artacho, D. Sanchez-Portal, P. Ordejon, A. Garcia, J. M. Soler, B.-k. Cheong, D. W. Shin, S. Lee, ECS Solid State Lett. 2013, 2, Phys. Status Solidi B 1999, 215, 809. Q75. [35] S. Jia, N. Shi, J. Shen, R. Wu, Q. Liu, Z. Song, M. Zhu, Phys. Status [6] J. Song, J. Woo, A. Prakash, D. Lee, H. Hwang, IEEE Electron Device Solidi B 2021, 258, 2100038. Lett. 2015, 36, 681. [36] M. Cobb, D. A. Drabold, R. L. Cappelletti, Phys. Rev. B 1996, 54, [7] S. Kim, J. Zhou, W. D. Lu, IEEE Trans. Electron Devices 2014, 61, 12162. 2820. [37] J. Y. Raty, V. v. Godlevsky, J. P. Gaspard, C. Bichara, M. Bionducci, [8] R. E. Simpson, P. Fons, A. v. Kolobov, T. Fukaya, M. Krbal, T. Yagi, R. Bellissent, R. Céolin, J. R. Chelikowsky, Ph. Ghosez, Phys. Rev. B J. Tominaga, Nat. Nanotechnol. 2011, 6, 501. 2001, 64, 235209. [9] O. Abou El Kheir, M. Bernasconi, Nanomaterials 2021, 11, 2382. [38] S. Clima, D. Garbin, K. Opsomer, N. S. Avasarala, W. Devulder, [10] Y.-H. Lee, P. J. Liao, V. Hou, D. Heh, C.-H. Nien, W.-H. Kuo, I. Shlyakhov, J. Keukelier, G. L. Donadio, T. Witters, S. Kundu, G. T. Chen, S.-M. Yu, Y.-S. Chen, J.-Y. Wu, X. Bao, C. H. Diaz, in 2021 B. Govoreanu, L. Goux, C. Detavernier, V. Afanas’ev, G. S. Kar, IEEE Int. Reliability Physics Symp. (IRPS), IEEE, Piscataway 2021, p. 1. G. Pourtois, Phys. Status Solidi RRL 2020, 14, 1900672. [11] P. Noé, C. Vallée, F. Hippert, F. Fillot, J.-Y. Raty, Semicond. Sci. [39] N. C. Murphy, R. Wortis, W. A. Atkinson, Phys. Rev. B 2011, 83, Technol. 2018, 33, 013002. 184206. [12] N. Saxena, A. Manivannan, Phys. Status Solidi RRL 2022, 16, [40] S. Caravati, M. Bernasconi, T. D. Kühne, M. Krack, M. Parrinello, 2200101. J. Phys.: Condens. Matter 2010, 22, 399801. [13] Z. Chai, W. Zhang, R. Degraeve, S. Clima, F. Hatem, J. F. Zhang, [41] W. Zhou, M. Paesler, D. E. Sayers, Phys. Rev. B 1991, 43, 2315. P. Freitas, J. Marsland, A. Fantini, D. Garbin, L. Goux, G. S. Kar, [42] E. Gulbrandsen, H. B. Johnsen, M. Endregaard, T. Grande, IEEE Electron. Device Lett. 2019, 40, 1269. S. Stølen, J. Solid State Chem. 1999, 145, 253. [14] B. Song, H. Xu, S. Liu, H. Liu, Q. Liu, Q. Li, Appl. Phys. A 2019, 125, [43] M. T. M. Shatnawi, C. L. Farrow, P. Chen, P. Boolchand, A. Sartbaeva, 772. M. F. Thorpe, S. J. L. Billinge, Phys. Rev. B 2008, 77, 094134. [15] C. Yoo, W. Kim, J. W. Jeon, E.-S. Park, M. Ha, Y. K. Lee, C. S. Hwang, [44] M. Xu, M. Xu, X. Miao, InfoMat 2022, 4. ACS Appl. Mater. Interfaces 2020, 12, 23110. [45] A. Padovani, L. Larcher, G. Bersuker, P. Pavan, IEEE Electron. Device [16] G. Liu, T. Li, L. Wu, Y. Chen, B. Liu, Z. Ma, S. Song, Z. Song, ACS Lett. 2013, 34, 680. Appl. Nano Mater. 2019, 2, 5373. [46] A. Padovani, D. Z. Gao, A. L. Shluger, L. Larcher, J. Appl. Phys. 2017, [17] N. S. Avasarala, G. L. Donadio, T. Witters, K. Opsomer, 121, 155101. B. Govoreanu, A. Fantini, S. Clima, H. Oh, S. Kundu, W. Devulder, [47] L. Larcher, F. M. Puglisi, A. Padovani, L. Vandelli, P. Pavan, in 2016 M. H. van der Veen, J. van Houdt, M. Heyns, L. Goux, G. S. Kar, in 26th Int. Workshop on Power and Timing Modeling, Optimization and 2018 IEEE Symp. VLSI Technology, IEEE, New Jersey 2018, p. 209. Simulation (PATMOS), IEEE, USA 2016, p. 128. [18] H. Li, J. Robertson, Sci. Rep. 2019, 9, 1867. [48] M. Pesiu, V. di Lecce, D. Pramanik, L. Larcher, in 2018 Int. Conf. [19] G. Liu, L. Wu, X. Chen, T. Li, Y. Wang, T. Guo, Z. Ma, M. Zhu, Simulation of Semiconductor Processes and Devices (SISPAD), S. Song, Z. Song, J. Alloys Compd. 2019, 792, 510. IEEE, USA 2018, p. 111. [20] A. Velea, K. Opsomer, W. Devulder, J. Dumortier, J. Fan, [49] K. A. Nasyrov, V. A. Gritsenko, J. Appl. Phys. 2011, 109, 093705. C. Detavernier, M. Jurczak, B. Govoreanu, Sci. Rep. 2017, 7, [50] G. Jegert, A. Kersch, W. Weinreich, P. Lugli, IEEE Trans. Electron. 8103. Devices 2011, 58, 327. Adv. Electron. Mater. 2023, 9, 2201224 2201224 (11 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH 2199160x, 2023, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/aelm.202201224 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License www.advancedsciencenews.com www.advelectronicmat.de [51] P. Vashishta, R. Kalia, I. Ebbsjo, Solid State Ionics 1989, 32, 872. [57] L. A. Agapito, S. Curtarolo, M. B. Nardelli, Phys. Rev. X 2015, 5, 011006. [52] M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, [58] T. D. Kühne, M. Iannuzzi, M. del Ben, V. v. Rybkin, P. Seewald, D. R. Hamann, X. Gonze, G.-M. Rignanese, Comput. Phys. F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, Commun. 2018, 226, 39. D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, [53] J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865. U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, [54] G. Román-Pérez, J. M. Soler, Phys. Rev. Lett. 2009, 103, 096102. T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, [55] T. Thonhauser, V. R. Cooper, S. Li, A. Puzder, P. Hyldgaard, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, D. C. Langreth, Phys. Rev. B 2007, 76, 125112. et al., J. Chem. Phys. 2020, 152, 194103. [56] A. Bouzid, C. Massobrio, M. Boero, G. Ori, K. Sykina, E. Furet, Phys. [59] J. Heyd, G. E. Scuseria, M. Ernzerhof, J. Chem. Phys. 2003, 118, Rev. B 2015, 92, 134208. 8207. Adv. Electron. Mater. 2023, 9, 2201224 2201224 (12 of 12) © 2023 The Authors. Advanced Electronic Materials published by Wiley-VCH GmbH