Ab Initio Study of Novel Phase‐Change Heterostructures

Abstract Neuromorphic devices constitute a novel approach to computing that takes inspiration from the brain to unify the processing and storage units. Memories based on phase‐change materials (PCMs) are potential candidates for such devices due to their non‐volatility and excellent scalability, however their use is hindered by their conductance variability and temporal drift in resistance. Recently, it has been shown that the utilization of phase‐change heterostructures consisting of nanolayers of the Sb2Te3 PCM interleaved with a transition‐metal dichalcogenide, acting as a confinement material, strongly mitigates these problems. In this work, superlattice heterostructures made of TiTe2 and two prototypical PCMs, respectively GeTe and Ge2Sb2Te5 are considered. By performing ab initio molecular dynamics simulations, it is shown that it is possible to switch the PCMs without destroying the superlattice structure and without diffusion of the atoms of the PCM across the TiTe2 nanolayers. In particular, the model containing Ge2Sb2Te5 shows weak coupling between the two materials during the switching process, which, combined with the high stability of the amorphous state of Ge2Sb2Te5, makes it a very promising candidate for neuromorphic computing applications.


Cell Optimization
We optimize the PCHs' cells using the QUICKSTEP code included in the CP2K package [1].We use a mixed Gaussian and plane-wave basis set to expand the Kohn-Sham orbitals (KSOs) and the charge density.Regarding the KSOs, we employ a Gaussian-type basis set of triple zeta (double zeta for Ti atoms) plus polarization quality.As for the charge density, we use a plane wave basis with an energy cutoff of 320 Ry.We model the valence configurations of Ge, Sb, Te, and Ti (4s 2 4p 2 , 5s 2 5p 3 , 5s 2 5p 4 , and 4s 2 3d 2 , respectively) using scalar-relativistic Goedecker-Teter-Hutter pseudopotentials [2].We also incorporate the Van der Waals (VdW) correction DFT-Grimme D3 [3] including the calculation of the C 9 coefficient, representing a 3-body interaction term.We minimize the cells through a BFGS minimizer that relies on the diagonalization of the entire Hessian matrix.We compute the stress tensor at each step.We set strict convergence criteria: the maximum geometry change (MAX DR) is 10 −5 Bohr, the maximum root mean square (RMS DR) geometry change is 5 • 10 −6 Bohr, the maximum force of each geometrical configuration (MAX FORCE) is 1.5•10 −6 Hartree/Bohr and finally the maximum root mean square force (RMS FORCE) of each geometrical configuration is 10 −6 Hartree/Bohr.Finally, we ensure that the optimization process preserves the cells' symmetry.In Table S1, we compare the relaxed lattice parameters of bulk GeTe, GST-225 and TiTe 2 at T = 0 K obtained with pure PBE, PBE with D2 corrections [4] and PBE with D3 corrections.

GST-225
GeTe TiTe Table S1: Comparison of a and c parameters of hexagonal bulk models of GST-225, GeTe, and TiTe 2 at T = 0 K.For GeTe the conventional hexagonal cell with 6 atoms is used.The wording "no corr" refers to the case of pure PBE calculations, while the others refer to D2 and D3 corrections added to the PBE functional.All the parameters are expressed in Å units.

Ab-Initio Molecular Dynamics Simulations
To determine the structural and dynamic characteristics of the confined PCMs in their crystalline (600 K), liquid (1300 K), and amorphous (300 K after quenching) phases we conduct ab initio molecular dynamics (AIMD) simulations using the CP2K package.In this context, we expand the charge density on a larger basis than the cell optimization setting, raising the energy cutoff to 500 Ry.The program employs periodic boundary conditions, with the Brillouin zone only sampled at the Γ point.We use the Born-Oppenheimer molecular dynamics (BOMD) and the second-generation Car-Parrinello molecular dynamics (SGCPMD) to ensure efficiency and accuracy.The latter is a technique developed by Kühne et al. [5], which is faster than BOMD and standard CPMD but intrinsically dissipative.The dissipation is described by the intrinsic damping coefficient "noisy-gamma" γ D , which depends on the system under consideration.γ D is extracted from preliminary simulations.A second, extrinsic damping coefficient, γ L , can also be added.After careful analysis, we set γ L = 4 • 10 −3 f s −1 and γ D = 3.5 • 10 −5 f s −1 for both PCH models.
We perform BOMD simulations with the Bussi thermostat [6] for warming up the system (first to 600 K, then to 1700 K).The Bussi thermostat allows for sampling the canonical ensemble through a stochastic rescaling of the velocities.In our simulations, the rescaling occurs every 10 timesteps (5 fs).This approach ensures a rapid temperature increase, while also generating optimized wavefunctions at the appropriate temperature.These wavefunctions are subsequently employed to initiate the SGCPMD.The timesteps in our simulations are 0.5 fs for BOMD-Bussi and 2 fs for SGCPMD.
We generate the crystalline phases at 600K through the following protocol: initially we conduct BOMD-Bussi simulations to warm up the crystalline PCHs from 0 K to 600 K covering trajectories that lasted 1 picosecond.After that, we equilibrate the systems at 600K for 20 picoseconds using SGCPMD.Finally, we generate trajectories spanning 30 picoseconds with SGCPMD, during which we calculate and analyze various observables.
We implement the two methods mentioned above to simulate also the PCH partial liquid phase at 1300K with some variations.When implementing the BOMD-Bussi simulation, we increase the temperature from 600 K to 1700 K, encompassing a trajectory that lasted 1 picosecond.Next, we apply the SGCPMD method holding the system at 1700 K for 10 picoseconds to accelerate the melting of the PCMs while ensuring that the CM retains its crystalline structure throughout the simulation.Following this, we reduce the temperature to 1300 K and equilibrate the system for 30 picoseconds.Finally, we generate a trajectory spanning 50 picoseconds, during which we compute and analyze all relevant observables.
To generate the amorphous phases, we employ SGCPMD only.We conduct an initial quenching with a rate of 10 15 K/s to bring the system's temperature from 1300K to 900K.After equilibrating the system at 900K for a few picoseconds, we perform a second quenching with a rate of 10 13 K/s until reaching a temperature of 300K.Then we generate trajectories of approximately 12 picoseconds and calculate and analyze the relevant observables (excluding the initial 2 ps for the analysis).

Analysis of the VdW gaps
Prior to cell optimization, the two PCHs exhibit a layered structure in their crystalline phase, with the total atomic density profile ρ(z) appearing as a sequence of Dirac delta functions, each with an amplitude corresponding to the number of atoms present in the layer.The thickness of the VdW gap in each PCH is determined by the peak-to-peak distance between the ρ(z) components relative to the planes of the PCM and the CM facing each other.Thus, in the GST-225/TiTe 2 heterostructure, the gap thickness is defined by the distance between the Te planes of both GST-225 and TiTe 2 facing each other.For GeTe/TiTe 2 , it is calculated from the distance between the Te planes of TiTe 2 and the adjacent Ge or Te plane of GeTe.However, the system's symmetry is reduced after cell optimization, and the layered structure may become distorted.This effect is more prominent at higher temperatures, as thermal motion further disrupts the crystalline order.Obviously, this effect reaches its maximum after partial melting of the PCM (as well as in the partly amorphous models at lower temperature).To calculate the VdW gap thickness in such instances, we determine the average position of the PCM and CM interfaces by identifying the local maxima of the ρ(z) histogram.The VdW gap is then defined as the spatial region, depicted in cyan color, between the PCM and CM ρ(z) nearest maxima, which are marked by dashed vertical lines in the bottom plots of Figures (S1) and (S2).We use a Gaussian filter to increase the precision of the estimate of the peak values.The coordination numbers (CNs) are a useful metric for studying the spatial correlations between the individual particles of a system.They are particularly helpful for characterizing the structure around an atom or molecule.Liquids and amorphous systems can be described as complex bond networks with welldefined local structures (coordination shells).However, characterizing these coordination shells in terms of a fixed set of distances between particles is ineffective, as in a liquid at finite temperature, the individual atoms are in constant motion, and distances constantly fluctuate.Additionally, quenched amorphous systems suffer from statistical variability of the coordination pattern.
A more effective approach to characterizing disordered structures is using a distribution function of distances in the coordination shell to analyze the average structure near the particles of the system.This distribution function, called radial distribution function (RDF), measures the probability of finding two particles a distance r apart and displays peaks at particular r values, which are characteristic of the structure, with peak widths primarily determined by the temperature and density of the system.To calculate the CNs, we first determine the RDF using AIMD simulations and then identify its first peak.We then integrate the RDF up to the position r cut of the first minimum to determine the number of atoms enclosed in a sphere of radius r cut centered on each particle (see Eq.3 in the main text).This protocol is also valid for partial RDFs.It follows that CNs depends sensitively on the choice of r cut and even a small variation can lead to significant changes in the CNs.This is why the RDF remains the primary benchmark for comparing our results with those indicated in [7,8].In fact, in the main text, we make it clear how many of our results align with those of these studies, despite having different CNs due to the use of different r cut .In this section, we present the CNs at T = 1300 K that we would obtain if we used the r cut values defined in the referenced works.As specified in Tables S2 and S3, the third and fourth columns report CN values that are quite compatible.Furthermore, we report in Figures(S3)    Table S3: Comparison of CNs concerning confined l-GST-225.Column 2 outlines the cutoff radii utilized in the study by Schumacher et al. [8], while column 3 presents the overall CNs reported by the same authors in their work.In column 4, we present the total CNs we obtained when applying the cutoffs specified in column 2.

Partial Angular Distribution Functions
Partial angle distribution functions (pADFs) serve as pivotal tools for characterizing the angular relationships among atomic triplets.Specifically, pADFs capture the statistical distribution of angles formed by atoms of species α, β and γ.The pADF quantifies the likelihood of observing a particular angle θ between the specified αβγ atomic triplet.This statistical metric enables a nuanced exploration of angular distributions, revealing insights into the local geometric arrangements and correlations inherent to the atomic structure.To compute the pADF we analyzed the trajectories obtained from AIMD simulations.We report in Figures (S5) and (S6) the pADFs relative to confined l-GST-225 and confined l-GeTe at T = 1300 K.

Coordination and q Order Parameter
The Errington-DeBenedetti q order parameter is a crucial metric in the study of phase transitions and structural properties of fluids.We report in Figures(S11) and (S12) the coordination of confined l-GST-225 and confined l-GeTe.We also report in Figures(S13) and (S14) the q-distributions relative to the prominent environments selected by the coordination in the previous two figures.

Diffusion Coefficient
The Einstein formula let us compute the diffusion coefficient D z along the vertical (z) direction: The average in Eq.( 1) designates an average over time and an average over several AIMD trajectories.To calculate MSD z (t) from a trajectory, it is essential to divide such trajectory into several displacements.To that aim, we compute MSD z (t) at equally spaced time points ∆t, 2∆t, • • • , T run ∆t, where ∆t represents the time step and T run designates the length of the entire AIMD simulation.In this way, we derive T run (T run − 1)/2 non-trivial forward displacements z(j) − z(i) that correspond to time lags (j − i)∆t, where 1 ≤ i < j ≤ T .The quantity z(i) denotes the z atomic coordinate at time i.While this definition involves many distinct forward displacements at small time lags, it entails very few such displacements at large time lags.However, since the time step is fixed, we can disregard this and refer to the time lag t when selecting j = (i + t).As a result, we align the definition of time lag with the time index and obtain the following formula for MSD z , which can be applied to each trajectory: ( To compute Eq.( 2), we wrote a Python code designed to handle different materials (GST-225, GeTe) with specific simulation parameters.The main functionalities include loading trajectory coordinates, dividing the simulation into sub-runs, determining the slab position for each atom, and calculating the time-lag MSD z for each slab.We begin by assigning each atom to a slab at time t = 0, relying on its z coordinate.Second, for all the atoms in a slab, we calculate the average in Eq.( 1) in the interval [0; t surv ] where the surviving time is defined in the main text.This operation is repeated every 1 ps (500 MD steps).Therefore, the second time window is [1, t surv + 1] and so on until the last, N w -th one, which is [N w , t surv + N w ], where t surv + N w is shorter but very close to the total duration of the simulation.The trajectories of atoms within each slab during a specific sub-run are stored in a trajectory matrix organized such that each row corresponds to a trajectory of a single particle, and each column corresponds to a time step within the specified sub-run.After the trajectories are organized in matrices, Eq.( 2) is calculated for each trajectory at each time lag over all pairs of time steps within the sub-run.This process is parallelized using multiprocessing to enhance computational efficiency.Then we perform the ensemble average marginalizing over the rows of the trajectory matrix.The result is a time-ensemble-average MSD z of a sub-run of a specific slab.Repeating these process for each sub-run gives the set of red curves depicted in Figures(S15) and (S16).From the average of these curves we obtain the black curves which are then fitted with a linear function in the region relative to diffusive regime.The result of the fit are listed in Table S4 and refers to Fig. (6) S5.We observe a slight increase of N T eT e with respect to [9,10].We observe some discrepancies between our ADF and the one of Ref.
[40] in the main text: in particular, we do not observe clear features at 90 • and 120 • degrees, contrary to their models.These discrepancies appear to be due to confinement effects.Calculation of the Te-Te-Te ADF for an independent a-GeTe bulk model at T = 300 K obtained with a machine-learning potential leads to more satisfactory agreement with the findings of Ref.
[40], as the two features at 90 • and 120 • degrees are better resolved.

Figure S1 :
Figure S1: Atomic density profile ρ(z) of the GST-225/TiTe 2 heterostructure.From left to right the T = 300 (amorphous), 600 and 1300 K cases are depicted.In each plot, from top to bottom, the Ge, Sb, Te, Ti and total atomic density profiles are shown.The black dashed lines in the bottom plots (PCH) depict the ρ(z) maxima defining the thickness of the VdW gaps.The VdW gaps are the region of space shaded in cyan and enclosed between two ρ(z) maxima.

Figure S2 :
Figure S2: Atomic density profile ρ(z) of the GeTe-225/TiTe 2 heterostructure.From left to right the T = 300 (amorphous), 600 and 1300 K cases are depicted.In each plot, from top to bottom, the Ge, Te, Ti and total atomic density profiles are shown.The black dashed lines in the bottom plots (PCH) depict the ρ(z) maxima defining the thickness of the VdW gaps.The VdW gaps are the region of space shaded in cyan and enclosed between two ρ(z) maxima.
and (S4) the CNs N αβ as a function of r cut of confined liquid GST-225 (l-GST-225) and confined liquid GeTe (l-GeTe) at T = 1300 K, where α and β denote atomic species.

Figure S3 :
Figure S3: CN N αβ of confined l-GST-225 as a function of r cut at T = 1300 K. Different colors refer to different αβ atomic pairs.Crosses represent the CN values that we obtained by means of the cutoff radii specified inTable 4 of the main text.

Figure S4 :
Figure S4: CN N αβ of confined l-GeTe as a function of r cut at T = 1300 K. Different colors refer to different αβ atomic pairs.Crosses represent the CN values that we obtained by means of the cutoff radii specified in Table 4 of the main text.

Figure S5 :
Figure S5: pADFs of confined l-GST-225 at T = 1300 K.The upper plot refer to all the Ge-centered triplets; the mid plot refer to all the Sb-centered triplets; the lower plot shows all the Te-centered triplets.Each pADF is normalized to the relative number of atomic triplets.The black dashed line is centered at 90 • while the grey dash line is centered at 109.5 • .

Figure S6 :
Figure S6: pADFs of confined l-GeTe at T = 1300 K.The upper plot shows all the Ge-centered triplets; the lower plot shows all the Te-centered triplets.Each pADF is normalized to the relative number of atomic triplets.The black dashed line is centered at 90 • while the grey dash line is centered at 109.5 • .

Figure S7 :
Figure S7: Sub-ALTBCs of confined l-GeTe at T = 1300 K.The left plot shows the Ge-centered ALTBC and the right plot shows the Te-centered ALTBC.

Figure S 8 :
Figure S 8: Sub-ALTBCs of confined l-GST-225 at T = 1300 K.The upper plots shows the Ge-centered (left) and Tecentered (right) ALTBCs.The lower plot shows the Sb-centered ALTBC.

Figure S9 :
Figure S9: Total and sub-ALTBCs of bulk l-GST-225 at T = 1300 K.The upper plots show the total (left) and Ge-centered ALTBCs.The lower plots show the Sb-centered (left) and Te-centered (right) ALTBCs.

Figure S10 :
Figure S10: Total and sub-ALTBCs of bulk l-GeTe at T = 1300 K.The upper left plot shows the total ALTBC and the upper right plot shows the Ge-centered ALTBC.The lower plot shows the Te-centered ALTBC.

Figure S11 :
Figure S11: Coordination of confined l-GST-225 at T = 1300 K.The upper plots show Ge (left) and Sb (right) coordination; The lower plot shows Te coordination.

Figure S12 :
Figure S12: Coordination of confined l-GeTe at T = 1300 K.The left plot shows Ge coordination; the right plot shows Te coordination.

Figure S13 :
Figure S13: Errington-De Benedetti local order parameter q of confined l-GST-225 at T = 1300 K.The upper plot shows the 4, 5 and 6-fold coordination of Ge atoms; the mid plot shows the 4, 5 and 6-fold coordination of Sb atoms.The lower plot shows the 3, 4, and 5-fold coordination of Te atoms.

Figure S14 :
Figure S14: Errington-De Benedetti local order parameter q of confined l-GeTe at T = 1300 K.The upper plot shows the 4, 5 and 6-fold coordination of Ge atoms; the lower plot shows the 3, 4, and 5-fold coordination of Te atoms.

Figure S15 :
Figure S15: Average MSD along z direction (MSD z ) of confined l-GST-225 in each slab as a function of the timestep at T = 1300 K. Red curves relate to MSD z of different sub-runs of the AIMD trajectory.The blue curve is the Average MSD z and the dashed black line is a linear fit of the MSD z in the diffusion regime.

Figure S 16 :
Figure S 16: Average MSD along z direction (MSD z ) of confined l-GeTe in each slab as a function of the timestep at T = 1300 K. Red curves relate to MSD of different sub-runs of the AIMD trajectory.The blue curve is the Average MSD z and the dashed black line is a linear fit of the MSD z in the diffusion regime.

Figure S17 :
Figure S17: Inset a) D z self-diffusion coefficient profiles in confined GST-225 at T = 1300 K. Ge and Sb atoms exhibit comparable trends, whereas the mobility of Te atoms is generally suppressed.Insets b-d) show various Te, Ge and Sb density profiles in the PCH.Each curves is generated by sampling the MD trajectory at different times: starting from a pre-melted crystalline phase (purple curves), passing through an overheated phase (Melting onset -melting complete on the colorbar) and finally reaching a thermalized liquid phase at 1300 K (Thermalization onset on the colorbar).These curves indicate that the decrease in mobility of the Te atoms is due to a layering effect induced in the PCM by the interface.Indeed, the two peaks corresponding to the Te layers of GST-225 next to the TiTe 2 trilayer disappear slowly (contrary to the peaks corresponding to the Ge and Sb layers) and non-negligible peaks are present even in the fully liquid state.This is in line with the fact that the presence of the Te layers of the confinement material favors the presence of Te-rich regions near the interfaces, resulting in lower mobility for this atomic species.

Figure S18 :
Figure S18: Snapshots of a-GST-225 (left) and a-GeTe (right) at T = 300 K. We depict Te, Ti, Ge and Sb atoms in yellow, brown, red and cyan respectively.

Figure S20 :
Figure S20: Coordination of a-GeTe at T = 300 K.The left plot shows Ge coordination; the right plot shows Te coordination

Figure S21 :
Figure S21: CN N αβ of confined a-GST-225 as a function of r cut at T = 300 K. Different colors refer to different αβ atomic pairs.Crosses represent the CN values that we obtained by means of the cutoff radii specified in Table9of the main text.

Figure S22 :
Figure S22: CN N αβ of confined a-GeTe as a function of r cut at T = 300 K. Different colors refer to different αβ atomic pairs.Crosses represent the CN values that we obtained by means of the cutoff radii specified in Table9of the main text.

Figure S23 :
Figure S23: Partial ADFs (pADFs) of a-GST-225 at T = 300 K.The upper plot shows all the Ge-centered triplets; the mid plot shows all the Sb-centered triplets; the lower plot shows all the Te-centered triplets.Each pADF is normalized to the relative number of atomic triplets.The grey dashed line is centered at 90 • while the black dashed line is centered at 109.5 • .

Figure S24 :
Figure S24: sub ADFs of a-GeTe at T = 300 K.Each sub-ADF is normalized to the relative number of atomic triplets.The black dashed line is centered at 90 • .

Figure S25 :
Figure S25: Partial ADFs (pADFs) of a-GeTe.The upper plot shows all the Ge-centered triplets; the lower plot shows all the Te-centered pADFs.Each pADF is normalized to the relative number of atomic triplets.The grey dashed line is centered at 90 • while the black dashed line is centered at 109.5 • .

Figure S26 :
FigureS26: Te-Te-Te ADF of confined a-GeTe at T = 300 K upon inclusion of the second Te-Te coordination shell at r cut ≈ 5.22 Å.We observe some discrepancies between our ADF and the one of Ref.[40]  in the main text: in particular, we do not observe clear features at 90 • and 120 • degrees, contrary to their models.These discrepancies appear to be due to confinement effects.Calculation of the Te-Te-Te ADF for an independent a-GeTe bulk model at T = 300 K obtained with a machine-learning potential leads to more satisfactory agreement with the findings of Ref.[40], as the two features at 90 • and 120 • degrees are better resolved.

Figure S27 :
Figure S27: Nearest neighbors of the Ge atoms with tetrahedral coordination, 6.3% are Ge, 71% are Te, and 22.6% are Sb.

Figure S28 :
Figure S28: Errington-De Benedetti local order parameter q of a-GeTe at T = 300 K.The upper plot shows the 4 and 5-fold coordination of Ge atoms; the lower plot shows the 3, 4, and 5-fold coordination of Te atoms.

6. 5
Distribution of Primitive Rings

Figure S29 :
Figure S29: Distribution of primitive rings in a-GeTe at T = 300 K.The vertical black bars in the ring histogram are error bars computed with the RINGS code.

Table 4
of the main text.

Table S2 :
[7]parison of CNs concerning confined l-GeTe.Column 2 outlines the cutoff radii utilized in the study by Weber et al.[7], while column 3 presents the overall CNs reported by the same authors in their work.In column 4, we present the total CNs we obtained when applying the cutoffs specified in column 2.

Table S4 :
of the main text.D z profile of the two confined l-PCMs at T = 1300 K as the slab number written in the first column varies.The region ascribed to confined l-GST-225 (l-GeTe) is divided into six (five) slabs.

Table S5 :
[9,10]ison of CNs concerning a-GST-225.Column 2 outlines the cutoff radii utilized in the works of Caravati et al.[9,10], while column 3 presents the overall CNs reported by the same authors in their work.In column 4, we present the total CNs we obtained when applying the cutoffs specified in column 2.