Device‐to‐Materials Pathway for Electron Traps Detection in Amorphous GeSe‐Based Selectors

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 a synergistic approach between techniques able to correlate specific material properties with device characteristics. Using a material‐to‐device multiscale technique, a reliable protocol for an efficient characterization of the active traps in amorphous GeSe chalcogenide is proposed. The resulting trap maps trace back the specific features of materials responsible for the measured findings, and connect them to an atomistic description of the sample. The metrological approach can be straightforwardly extended to other materials and devices, which is very beneficial for an efficient material‐device codesign and the optimization of novel technologies.


Device-to-Materials Approach
The characterization of the material parameters (e.g., composition, thickness, morphology, stoichiometry) that mostly affect the device response is a tremendous task, especially when it deals with the structural and electronic properties of amorphous systems. Even though the computational materials discovery is a major task for the design of new technology, this top-down approach is often infeasible due to economic constraints of the sample's fabrication process. Thus, an alternative solution is to follow the inverse device-to-materials pathway that extracts the intrinsic properties of the materials from the experimental data. Since the electrical characteristics of chalcogenide selectors and their temperature dependence are well explained by a trap-assisted carrier model, the identification of defects/traps becomes crucial in the process.
Here, we developed the three-step approach, summarized in Figure 1: Step 1: Experimental setting involves the sample preparation, the device fabrication, and the measurement of the electrical characteristics (I-V, G-V).
Step 2: The experimental output is used as the input of a specific defect spectroscopy (DS) functionality, implemented in the device-oriented modeler Ginestra (R) , [22] which is specifically developed to extract intrinsic material properties and characteristics of the defects (thermal ionization and relaxation energies, distribution within the bandgap) from the electrical measurements.
Step 3: The resulting defect mapping is interpreted in terms of ab initio simulations, based on density functional theory.
The core of the process is the DS functionality that relies on the accurate modeling and simulation of the carrier-defect interactions that are at the bases of the defect-assisted trapping and transport (trap-assisted tunneling, TAT), [23] as described in the framework of the multiphonon theory which accounts for carriers-phonons coupling and lattice relaxation processes [23][24][25] (see the Experimental Section). The DS approach comprises four main steps: i) The measured device/structure is recreated within the Ginestra (R) simulation platform, starting from the known device and material parameters (layers thickness, bandgap, dielectric constant, etc.). ii) The measured experimental data are loaded into the DS tool, as a function of the temperature, voltage, and frequency (the latter for G-V data only). In order to guarantee high accuracy and uniqueness of the results, the DS methodology has to be applied to multiple electrical data (e.g., I-V and G-V as in this paper) measured under multiple conditions (voltages, temperatures, and frequencies). iii) Definition of the material and trap parameters to be extracted as well as their variation range to be considered during the optimization procedure. This information is defined in the input file. iv) An automatic optimization procedure accurately reproduces the experimental data by varying the selected www.advelectronicmat.de parameters within the specified intervals by minimizing the overall error between the experimental and the simulated electrical characteristics. A modified root mean square error is used for this purpose. At the end of the procedure, the DS module provides the energy-versus-position (E, z) map of the defects responsible for the observed voltage, temperature, and frequency dependence of the experimental data (uniform spatial distributions in the whole GeSe thickness are considered in this work for simplicity). It also provides selected material/ device parameters (e.g., work function and energy bandgap) as well as the defects relaxation energy (E REL ), which represents the energy required for the lattice relaxation process to occur when the defect state changes its charge state. [23][24][25] E REL is estimated from the temperature dependence and determines the charge capture/emission time constants and their temperature dependence and it is linked to the defects' atomic structure. Therefore, it can be used, together with the extracted thermal ionization energy band, to identify the nature of the defects responsible for the observed electrical data. [23] 3. Results

Step 1: Device Fabrication and Experimental Data Generation
Three different Ge x Se 100−x OTS layers were prepared and integrated in two-terminal MIM selectors, by combining materials with different thicknesses (10 and 20 nm) and Ge concentration (x = 50 and 60) as reported in Table 1. The amorphous Ge x Se 100−x films were prepared by room temperature physical vapor deposition (PVD). The layers were integrated into TiN/ GeSe/TiN two-terminal structures in a 300 mm process flow and passivated with a low-temperature back end of line process. The details of fabrication and measurement sets are reported in Section S1 (Supporting Information). Each current-voltage (I-V) or conductance-voltage (G-V) measurement was performed on a fresh device, and multiple (at least three) samesized devices were measured at each measurement condition.

Step 2: Device Modeling and Trap Mapping
The proposed DS technique is applied to the I-V and G-V data measured on the GeSe samples D21 and D22, described above; results for the D16 system are reported in Section S1 (Supporting Information). Since the conductance is found to be independent of the frequency, only the data measured at 2 kHz are considered for the DS procedure.
GeSe has a small bandgap from ≈0.8 to ≈1.5 eV depending on the stoichiometry, which is indicative of a small Schottkylike energy barrier at the TiN/GeSe interface (barrier larger than the bandgap would correspond to an Ohmic-like contact). The resulting electrical characteristics for devices D21 and D22 are shown in Figure 2a,b and c,d, respectively. An excellent agreement between experiments and simulations is obtained for both stacks in the whole range of voltages and temperatures considered. The consistency of the applied model and of corresponding results has been checked against the thickness change, by comparing the data for devices D22 and D16 (see Figures S1 and S2, Supporting Information).
Starting from the I-V and G-V characteristics, we extracted the material and trap parameters, as reported in Table 2. The Ge-rich system has a lower bandgap (1.1 eV) than the stoichiometric one (1.45 eV), in agreement with previous works [26,27] and consistent with the reduction of the experimental gap observed in the D21 device. The DS defect maps reported in 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 the sample. The trap density (N T,i ) is similar for electrons and holes, and it is higher for the Ge-rich case than stoichiometric one.
However, important differences can be detected in the two samples. In the case of Ge 50 Se 50 trap bands are slightly narrower and located on average at ≈0.4 and ≈0.2 eV from the valence and conduction band edge, respectively. The Ge 60 Se 40 compound is characterized by broader trap bands located on average at ≈0.1 and ≈0.3 eV from valence band top and the conduction band bottom. At zero bias traps are either fully occupied or empty, with the Fermi level lying in the middle of the gap. Once the external bias is applied, both carrier types contribute to conduction, with a net prevalence of electron conduction through the upper trap band (n-type). However, conduction type may depend on structural parameters of the device, such as thickness, effective bandgap, electrode work-function. Indeed, the analysis for the D16 sample exhibits an equal contribution of both hole and electron carriers.
In order to investigate the role of traps in the transport mechanisms within the GeSe layer, we recalculated the I-V characteristics of the D22 sample, by selectively deactivating specific scattering terms (Figure 4a). In particular, we considered the electrical characteristics resulting i) without the inclusion of the defect contributions, ii) without the trap-to-trap transitions, and iii) without the trap-to-band transitions. The results at different temperatures for the D22 sample are summarized in Figure 4b contributions has smaller effects on the I-V characteristics: trap-to-trap seems 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 that the most relevant contribution for current stems from the presence of the trap bands within the bandgap, while the other mechanisms are responsible for the specific details of the resulting electrical characteristics. This can be justified in terms of the simultaneous small bandgap and high trap density of both systems that make trap-to-trap and trap-to-band energetically similar. Another important parameter that can be extracted through the DS tool is the relaxation energy (E REL ), which is directly connected to the atomic nature of the trap and is typically derived from the temperature dependence of the charge transport (i.e., of the leakage current, as shown in Section S1, Supporting Information). However, an accurate extraction of this key parameter is possible in principle only if the observed current-temperature dependence originates from purely trap-assisted transport. In the case of the thick, small bandgap GeSe layers considered here, this could be not the case and other mechanisms, such as tunneling from/to the electrodes may play a role: simulations show that the charge transport occurs by means 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 dependent on the temperature, the E REL value of 0.1-0.5 eV for both electrons and holes that is consistently extracted for the different GeSe samples has to be verified also with other methods (e.g., ab initio calculations).
The extracted trap density is consistent among samples and within the range typically reported in the literature. Besides characterizing the properties and distribution of GeSe traps, the DS tool also allows the extraction of the GeSe effective thickness (t GeSe ). As reported in Table 2, the extracted value of t GeSe is very consistent for both stoichiometries and in very good agreement with the nominal thickness.
This analysis demonstrates the DS capability to characterize the trap distribution of samples, and to discern among different compositions, temperatures, and geometrical properties of the device. However, the DS tool is not sufficient to unravel the structural properties and chemical origin of those traps. For crystalline systems, the analysis of defect states is routinely done by considering the effects of point defects (e.g., vacancies, interstitial, impurities) that usually cause the appearance of localized states in the semiconductor band gap. In the present  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 temperature, on the different OTS devices. E G is the (mobility) bandgap, t GeSe is the film thickness, E T,i and ΔE T,i are average and spread of the traps thermal ionization energy, E REL,i is the traps relaxation energy, and N T,i is the traps density for electrons (i = el) and holes (i = h), respectively. case, the amorphous structure of GeSe prevents from using this intuitive approach and demands for an extensive microscopic investigation of realistic disordered structures.

Step 3: First Principles Calculations
In order to interpret the DS results from a microscopic point of view, we performed a high throughput study of the structural and electronic properties of GeSe compounds based on first principles calculations. Given the intrinsic disordered nature of the amorphous materials, the understanding of the interplay between local-order structures and electronic trap states requires the capability of treating large systems at the atomistic level and of having access to the corresponding electronic structure. In the absence of direct experimental results about the atomic structure of the glassy phase at room temperature, the identification of a reliable geometry is a major challenge for theory. In order to obtain realistic models of the amorphous structure without resorting to experimental information or adjusted interatomic potentials, we carried out an extensive set of independent simulations performed with different multiscale techniques and different codes, including molecular dynamics (MD) using interactions from both classical potentials (LAMMPS code [28] ) and density functional theory (SIESTA [29,30] and QuantumEspresso [31] ). As a starting point, we considered the atomic model of Ge 50 Se 50 resulting from classical MD simulations on a large simulation cell including 4480 atoms (labeled big cell). The structure is obtained through nanosecond-long melt-andquench cycles (5 K ps −1 cooling rate), following the procedure described in ref. [32]. This approach allows the system to rearrange at the local-and medium-order structures avoiding unphysical "frozen-liquid" structures, and it has been demonstrated to generate a realistic model for the amorphous GeSe compound. [32,33] The increase of Ge content from 50% to 60% in Ge x Se 100−x causes a systematic change of the average coordination number of Ge and Se species that moves from 3:3 (as  www.advelectronicmat.de in the crystal) to a phase where Ge has a mixed fourfold (as in Ge-bulk) and threefold coordination. The different short-range structures and their spatial arrangement in the medium range are expected to have a critical role in the formation/appearance of the electrical trap states (see below). We assumed the big cell structure as the reference testbed for the further comparisons.
We first exploited the capability of the SIESTA code to treat large systems by simulating one 4480 atoms structure extracted from the MD simulation at 300 K (i.e., one frame from MD trajectory). The structure was fully relaxed using the DFT forces computed in SIESTA, through a total-energy-and-force optimization approach to allow the system to readjust the atomic distances and angles on the basis of the DFT potential. The comparison between the original MD and the DFT results is summarized in Figure 5a-c, which reports the partial radial distribution function g(r) of the (a) GeGe, (b) GeSe, and (c) SeSe pair distances. Except for the expected minor differences due to bond-length optimization, the optimized DFT results (labeled DFT big cell) well reproduce the classical ones (labeled MD big cell).
However, a meaningful analysis of an amorphous system requires the investigation of several frames to gain a statistical average of the structural properties. This makes the use of such a large cell practically unfeasible for extensive DFT simulations. This opens the problem on how to generate smaller atomic models that are simultaneously structurally reliable and numerically affordable. We generated two independent sets of structures, by exploiting two different approaches and DFT codes.
In the first one (labeled m&q) we performed a sequence of melting-and-quenching cycles on model structures, by using the SIESTA code. Initial Markov chain Monte Carlo simulations have been used to generate the disordered starting structures of Ge x Se 100−x (x = 50, 60), with N = 336 atoms in a cubic box of side L = 21.0 Å. Starting from these initial structures, the amorphous models have been obtained through meltingand-quenching cycles, carried on by using ab initio MD. The quenching step is particularly critical to obtain trustable amorphous results. Preliminary results (see Section S2, Supporting Information) indicate that i) the relaxation of volume under constant pressure, ii) the use of a relatively large basis set (in particular, double-ζ plus polarization [DZP] [34] ) during the quenching procedure, and iii) a low quenching rate are necessary steps for the building of realistic GeSe amorphous structures. The details of the quenching-and-melting approach are summarized in Table 3. Further computational details and accuracy tests are reported in the Experimental Section and in the Supporting Information.
In the second approach (labeled cut), we extracted from the original MD trajectories on the large (4480 atoms) Ge x Se 100−x (x = 50, 60) models, 12 smaller cubic boxes containing 500 atoms each, with the same stoichiometries and densities of the parent MD systems. These smaller cells (cut from the larger box and imposing periodic boundary conditions) were relaxed using classical MD at a constant temperature of 300 K for 10 ns, in order to allow the broken bonds at the borders to rearrange and match the new periodic boundary conditions. Each system has been further optimized at the DFT level, by using the QuantumEspresso code [31] (see the Experimental Section and Supporting Information for further details).
The averaged radial distribution functions of model structures resulting from the m&q and cut approaches are shown in Figure 5d-f. A direct comparison clearly indicates that the final configurations resulting from the two approaches resemble the original big one, in particular the radial distribution functions are very similar to the original results ( Figure 5, top panels). This confirms that we generated atomically different www.advelectronicmat.de but statistically equivalent systems, whose structural characteristics are independent from the procedure and the numerical details used to realize the model. Similar considerations hold for the corresponding electronic structures, even though in this case the specific choice of the numerical details (e.g., exchangecorrelation functional) may affect the quantitative value of parameters such as the bandgap. However, as evident from the comparison of the electronic density of states (DOS) reported in Figure S8 (Supporting Information), both DFT approaches give the same qualitative results. This assures that the trap states analysis is robust and independent from the specific choice of the used DFT code and functional.
On the basis of these results, we performed a set of calculations to generate model structures for the two different stoichiometries corresponding to the experimental samples D21 and D22, i.e., Ge 50 Se 50 and Ge 60 Se 40 . The Se-rich counterpart Ge 40 Se 60 has also been simulated for comparison (see the Supporting Information). For all GeSe compounds, several (12) atomic structures have been generated to have a statistical meaningful insight of the amorphous characteristics. The complete set of DOS plots for each simulated structure is summarized in Figures S9-S11 (Supporting Information). The electronic properties of the three stoichiometries have similar spectral features: the valence DOS plots are characterized by three distinctive peaks (labeled A, B, and C in Figure S12, Supporting Information) in agreement with recent X-ray spectroscopy. [35] Peak C corresponds to the energy region going from −14 to −12 eV and has mostly a Se(4s) character. The central region, from −10 to −6 eV (peak B), is dominated by the Ge(4s) contribution, while the region from −5 to 0 eV (peak A) has a mixed Se(4p)-Ge(4p) character, with a major contribution from Se states, which varies with the Ge/Se ratio. The lowest energy part of the conduction band has a mixed GeSe behavior as well, but here the major contributions derive from the Ge atoms, as a function of the composition. This characteristic peak distribution is typical of GeSe compounds, such as crystalline GeSe and GeSe 2 . [36][37][38] Despite these overall common aspects, the gap regions, which are the most relevant for the transport properties, exhibit different characteristics depending on the stoichiometry. Figure 6 shows the DOS plots resulting from the ensemble average on the 12 structures of Ge 50 Se 50 and Ge 60 Se 40 , respectively, as resulting from the cut approach. A Hubbardlike potential has been included to both Ge and Se atoms, in order to correct the severe underestimation of the bandgap due to standard DFT approach (see Section S2, Supporting Information, for further details). The zero-energy reference is set to the corresponding Fermi level. Red diamonds indicate the inverse participation ratio (IPR) of single particle orbitals, [39,40] which is a measure of the spatial localization of the electronic states: the higher the IPR the more localized the states.
The mobility gap E g is determined as the energy-level border separating the localized and extended bands, i.e., as the energy region where the IPR increases with respect to the deep valence/conduction bands. Increasing the Ge amount in the material, the mobility gap E g decreases from 1.4 eV (x = 50%) to 0.7 eV (x = 60%) in agreement with the experimental data and the DS results. The density, the localization and the energy distribution of defect states in the gap strongly depend on the stoichiometry. The nature of the mobility-gap states is found to be different from the Valence-Alternating-Pairs (VAP) model for Se-rich GeSe. [41][42][43] In the case of Ge 50 Se 50 [panel (a)], the gap is almost empty with a few defect states that accumulate close to the edges of the bandgap, forming two trap bands (dashed cyan lines) for holes and for electrons, in agreement with the DS results of Figure 3. The states forming these defect bands are spatially localized (Figure 7a,b) but involve different local coordination depending on n-type or p-type character of the trap. A similar picture also holds for the Se-rich Ge 40 Se 60 compound ( Figure S13, Supporting Information), Table 3. Melt-and-quench (m&q) protocol for generating the amorphous structures. SZ and DZP indicate single-ζ basis and doublepolarized-ζ basis sets; NPT indicates isothermal-isobaric ensemble.
Initial structure Generate initial structure with MCMC (336 atoms). Then relax structure (cell and atom coordinates)  Figure 3b. However, in this case the in-gap states are spatially delocalized all over the cell (Figure 7c,d), so they cannot properly be considered as trap states, but rather as conductive states. This goes in parallel with the enhancement of the current measured in sample D22 ( Figure 2). Figure 7 shows the single-particle charge density plots of a few selected defect states of the GeSe systems for both occupied (a, c) and empty (b, d) states. Gray polyhedrons indicate ordered threefold and fourfold local-order structures centered on Ge atoms, that are among the most frequent local arrangements found in these systems. [33] From the analysis of the spatial distribution of the in-gap states, we draw a few important observations: i) at odd with the usual crystalline case, the localized states of these amorphous structures are not localized around the single atomic site (e.g., the point defect), but involve a few bonds. ii) The defect states often appear at interconnections among a few threefold and fourfold polyhedrons, indicating a key role of Ge coordination on the formation of trap states, in agreement with previous theory reports. [38] iii) In the stoichiometric case [panels (a) and (b)], the abundance of Se favors the formation of GeSe bonds arranged in ordered local structures (local order). This justifies the presence of the reduced number of trap states in the Ge 50 Se 50 system. On the contrary, the deficiency of Se in the Ge-rich case causes the formation of homopolar GeGe bonds, which tend to be interconnected across the cell instead of organizing in folded structures (i.e., higher delocalization). Note that the number of polyhedrons decreases as the Ge content increases. iv) Since the top of the valence band (peak A, Figure S12, Supporting Information) has a prevalent Se character, for all stoichiometries (including Ge 60 Se 40 ) 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 heuristically be extended to the case of the stoichiometric compound (x = 50%) in view of the comparison with the crystalline counterpart, the same concepts are completely ill defined in the case of off-stoichiometric systems. Indeed, the 20% of compositional unbalance as in Ge 40 Se 60 or Ge 60 Se 40 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 crucial 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 common feature of the latter. All electron defect states involve Gehomopolar bonds. However, not all the GeGe bonds give rise to a state in the gap. The question is: what is the origin of these states? The answer stems from the structural analysis of the medium range order of the system, and in particular from the formation of GeSe complexes with different coordination.
To illustrate this point, we considered the trap states relative to a single snapshot of Ge 50 Se 50 , assumed as a testbed of the system. Similar arguments can be extended to the other GeSe cases. We focused on one occupied (T h ) and one empty (T el ) trap state, which decorate the band edges of the system (Figure 8), and we analyzed the coordination of the atoms closely involved in these trap states. The occupied defect state (T h ) is localized around three close Ge atoms that are not forming homopolar bonds, but that are simultaneously bonded to Se atoms having different coordination, namely, threefold (III) and fourfold (IV) Se atoms. This Se-coordination mix results in the formation of an unbalanced and localized charge distribution around the Ge atoms that causes the appearance of the trap state in the gap. The T el state is localized across a triple GeGeGe chain, where one Ge (the one on the right) is bonded to three inequivalent Se atoms, having II, III, and IV coordination, respectively. Again, the trap state originates from the off-coordination of the bonded Se atoms, and not from the Ge itself or from the formation of the Gehomopolar bonds. To confirm this statement we analyzed an extended state (labeled No trap in Figure 8), where we identified the formation of a GeGe bond. However, both Ge are bonded to homo-coordinated Se atoms (Se III ). This prevents the localization of charge and the formation of a trap state, in agreement with extended machine learning investigations on GeSe alloys. [44] We can conclude that trap states originate from local electron under/over-coordination, in agreement with previous general statements. [18,38,44] Yet, our results indicate that off-coordination is not associated to a single atomic site as usually suggested, but rather to nm-extended complexes (medium range) that include both chemical species. The multielement/multibond complexity of the trap origin justifies the different interpretations appeared in the literature, often based www.advelectronicmat.de on simplified structural models, such as homopolar bonds or folding of single chemical elements (e.g., Ge-fold structures).
Since trap-assisted transport implies the hopping of charge to/from localized states and the 0/−1 (0/+) charge-state transition for electron (hole) carriers, we investigated the effect of including an extra charge to the system. For each initial snapshot of the Ge 50 Se 50 system, we simulated the case of an extra hole or an extra electron in the unit cell. We fully relaxed the charged structures and we averaged the resulting electronic structures, as summarized in Figure 9.
Even though the two characteristic trap bands of Figure 3 are still visible, the inclusion of the extra charge in the system imparted a redistribution of the defect energy within the bandgap that results more occupied with respect to the neutral case. We observe also an increase of the average IPR of the defect states, which corresponds to an increase of the spatial charge localization favored by the local structural relaxation (i.e., polaron-like effect). The direct evaluation of the chargestate transition energies for amorphous systems is an open challenge within the solid-state community, and goes beyond the aim of this work.

Conclusion
We presented a multiscale device-to-materials approach for the identification and the analysis of the trap distribution in electronic devices. Starting from the experimental electrical characteristic of two-terminal TiN/a-GeSe/TiN selector devices, we fully investigated the properties of the active material at the microscopic level and its transport response. This approach allowed us to connect the performances and reliability of the experimental devices to the composition and the stoichiometry of the system. In particular, we unravel the role of charge traps in amorphous chalcogenides and the relation of these defect states with local order of the system. The proposed D2M protocol is general and it can be applied to any CMOS, logic or memory device, independently from the  www.advelectronicmat.de structure, the composition, and the chemical species of the active materials.
Our learning approach combines characterization and metrology aspects that may accelerate the optimization of material-design capabilities by exploring new materials and new compounds, starting from the desired electrical performances of the device.

Experimental Section
Fabrication and Characterization: Amorphous Ge x Se 100−x films were prepared by room-temperature PVD. TiN/GeSe/TiN MIM capacitors were integrated in a 300 mm CMOS process flow. Electrical characterization was performed with a Keithley 2536 SMU and an Agilent E4980A LCR meter. For further details see Section S1 (Supporting Information).
Device Modeling: Device simulations were performed with the Ginestra software package (Ginestra), a multiscale simulation platform that selfconsistently describes all the mechanisms relevant for logic, [23,45,46] and memory devices. [47,48] Non-defect-related conduction mechanisms such as direct/Fowler-Nordheim tunneling, thermionic emission, and drift were considered consistently with defect-assisted charge trapping and transport (TAT) model that controls the conduction in a wide variety of materials. [23,45,49,50] TAT is described in the framework of the multiphonon theory accounting for carriers-phonons coupling and lattice relaxation processes, which depends on the atomistic structure of the defect and is described by the relaxation energy parameter (EREL). Hydrodynamic transport, current-induced power dissipation and temperature increase, distortion and breakage of atomic bonds promoted by field, temperature and electron injection, [46] diffusion of atomic species, [47] ferroelectric switching, [48] and phase change were also self-consistently considered, thus enabling the modeling of material modifications associated to the electrical stress and specific material properties. Calculations were performed considering the local potential given by the applied bias and the defect charge state and occupation. Such a comprehensive and selfconsistent description of carriers transport and trapping mechanisms is at the basis of the DS technique [21] that was developed for extracting defects and properties of material used in modern and novel logic and memory devices.
Materials Simulations: Classical MD Simulations: The Ge 50 Se 50 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 employed for the description at the atomistic level of GeSe 2 or InP systems. [51] The classical MD simulations were carried out by using 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 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 repulsive interaction between atoms was included to make their distribution more homogeneous in the cell, and potential which depends on the Ge and Se coordination, to impose a Ge:Se coordination of 4:2. The resulting geometry was then relaxed prior to the meltingand-quenching procedure by using a conjugate gradient algorithm. Ab initio molecular dynamics (AIMD) were then carried out within the isobaric-isothermal (NPT) ensemble at zero external pressure with an annealing algorithm, by rescaling the velocities and the pressure at each time step of Δt = 1 fs and the characteristic relaxation time of τ relax . The simulations started by equilibrating each initial configuration at 1200 K for 12 ps, at which the system melts and loses the memory of the initial configuration. After equilibration, each system was quenched to 800 K over a period of 3 ps; this corresponds to an average cooling rate of 133 K ps −1 . During the melting and the quenching to 800 K procedure, the τ relax was set to 1000 fs. The obtained configuration at 800 K was quenched again to 400 K in 3 ps. At this temperature, the structure was equilibrated additionally for 3 ps using now the DZP and a τ relax of 750 ps.
The relaxation and AIMD calculations were carried out with the SIESTA code. [29,30] SIESTA employs local basis sets, and norm-conserving pseudopotentials from the pseudo-dojo dataset, [52] in which the Se and Ge 4s and 4p electrons are considered as valence states. The generalized gradient approximation (GGA)-Perdew-Burke-Ernzerhof (PBE) [53] 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 DZP basis functions were used. The basis set was generated with an orbital-confining scheme with a cutoff radius corresponding to the energy shift of 0.001 Ry. [34] Only the Γ-point was used for the Brillouin zone integration. The convergence criterion for the geometry optimization was set to a maximum force of 0.04 eV Å −1 on each atom. Extended accuracy test on the choice of the basis set is reported in Section S2 (Supporting Information).
Cell Cutting (Cut) Approach: For each stoichiometry, 12 atomic structures of Ge x Se 100−x each including 500 atoms were generated. The initial geometries were extracted from the large (4480 atoms) structures generated with the classical interatomic potentials, discussed above. Each structure was then optimized by using a total energy and forces approach based on DFT calculations, as implemented in the QuantumEspresso package. [31] The electron-ion interactions were treated within the projector augmented wave (PAW) framework. Single particle wave functions and charge were expanded in planewaves up to energy cutoff of 30 and 300 Ry, respectively. The exchange-correlation functional was described by the vdW-DF2-b86r exchange-correlation functional [55] that takes into account the van der Waals contributions which are known to play a crucial role in amorphous chalcogenides. [56] All the atomic positions were allowed to move until the maximum ionic forces were less than 0.002 eV Å −1 . Integrations over the Brillouin zone were sampled at the Γ-point only. In order to overcome the well-known bandgap underestimation of DFT, a DFT+U correction was included to the electronic structure, within the pseudohybrid Hubbard density functional approach (namely ACBN0 [57] ) at single point scf-calculations. The resulting Hubbard-like values are: U Ge = 0.16 eV and U Se = 2.28 eV. The detailed analysis of the effect of the inclusion of the Hubbard-U potential in comparison with the hybrid HSE06 approach is reported in the Supporting Information.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.