Quench-rate and size-dependent behaviour in glassy Ge₂Sb₂Te₅ models simulated with a machine-learned Gaussian approximation potential

Current electronic memory devices are based on silicon complementary metal-oxide-semiconductor (CMOS) technology, but the roadmap for Moore's law continuing for such non-volatile (flash) memories, consisting of floating-gate (FG) MOS, field-effect transistor (FET) devices, is predicted to be nearing its end [1]. This is due to memory volatility occurring in such devices via electrons, trapped on the FG, tunnelling away through the very thin layers of dielectric, separating the FG and the semiconductor interface, associated with extreme down-scaling of device features. Thus, a new non-volatile memory (NVM) technology is required in order to bypass this technological-roadmap block, which is capable of further down-scaling. Phase-change random-access memory (PCRAM) is one such contender [2], and which is already being commercialized in the Micron-Intel Optane^TM solid-state memory device, based on 3D-Xpoint^TM cross-bar technology, and widely believed to use a PC material as the NVM element. PCRAM materials exhibit rather unique properties: they are electrically conducting in the crystalline (c-)state (logic state '1') but electrically resistive in the glassy (g-)state (logic state '0'); they can be switched repeatedly and reversibly between metastable c- and g-states (the phase change), with the rate-limiting crystallization times being of the order of a few nanoseconds; and the PC behaviour scales down to ∼ 2 nm.

The canonical PCRAM material is Ge₂Sb₂Te₅ (GST225), a composition along the pseudo-binary GeTe-Sb₂Te₃ tie-line of ternary Ge–Sb–Te semiconductor alloys. Although GST225 exhibits PC behaviour which is sufficient for PCRAM devices, the properties are sub-optimal in many regards [1]: the crystallization time is rather long (∼ 10 ns); the electrical-resistance contrast between g- and c-phases is insufficient for significant multilevel programming, in which multiple metastable resistance levels (corresponding to different g-/c-volume proportions) can be written into the same memory cell, thereby permitting multi-bit per cell operation; resistance 'drift', a slow time-dependent increase in electrical resistance, is observed in the g-state, which is deleterious to multilevel operation. Thus, it is imperative to search for new PC compositions with improved behaviour, in order to develop next-generation PCRAM devices. But, in order to do so effectively, it is necessary first to understand thoroughly, at the atomistic level, the origin of the behaviour of GST225, so that any deficiencies in its performance can be understood and ameliorated in next-generation PC materials.

If expensive random-structure-searching efforts to find new PC materials are to be avoided, and a more rational materials-design approach adopted, one needs to make use of the knowledge gained about the atomistic origin of PC behaviour in canonical PCRAM materials by constructing realistic structural models of the technologically relevant phases. Such an approach can help to direct compositional and structural searches for new materials in order to ameliorate their property limitations and is particularly important for glassy materials which exhibit atomic disorder. The current state of the art, for generating (electronically) accurate models of glasses, is ab initio molecular dynamics (AIMD), as driven by density-functional theory (DFT) approximations. A number of such simulations of PC materials have been published, on the canonical GST225 composition [3–7] as well as on several others [8–10]. However, because such calculations are so computationally intensive (scaling cubically (O(N³)) with the number of atoms, N, in the simulation box), this limits the maximum size of such models to ∼ 500 atoms, the maximum simulation time to ∼ 1 ns, and the number of simulations achievable in a given study to typically just one (thereby raising questions about the statistical significance of the calculated properties of such single, small models).

In order to circumvent these difficulties, it is necessary ideally to employ linear-scaling, O(N) interatomic potentials, which are able to approach DFT levels of energy and force accuracy. Such a potential was first constructed for the binary PC material, GeTe, by Sosso et al [11] using a neural-network, machine-learning method of fitting to the DFT-calculated electronic potential-energy surfaces of representative atomic configurations. We were the first to develop an O(N) potential for the ternary Ge-Sb-Te system, in particular Ge₂Sb₂Te₅ [12], using a different machine-learning approach (Gaussian process regression [13], rather than a neural network). The resulting Gaussian approximation potential (GAP) is DFT-accurate, linear-scaling and reasonably chemically transferable (for other GST-based compositions) [12]. The O(N) capability of the GAP allowed for the simulation of an ensemble of 30 × 315-atom MD models of GST225 [14], with this model size dictated by the models being subsequently used for DFT calculations (with hybrid exchange-correlation functionals) of their electronic properties, specifically for a statistical examination of the occurrence of gap states.

In this paper, we turn to the other two distinct advantages of using (DFT-accurate) O(N) potentials, such as the GAP for GST225, viz. that long timescales can be accessed through simulations, and that sufficiently large models can be simulated in order to overcome finite-size effects. In particular, the ability to run longer molecular-dynamics trajectories translates into smaller values of the quench rate that can be used in vitrification simulations, as carried out from the equilibrated liquid; and the capability of simulating ultra-large models means that size effects on computed properties can be examined systematically.

Glassy GST225 structures were generated using classical MD simulations in cubic simulation cells with periodic boundary conditions. All the MD simulations were performed using the open-source LAMMPS (large-scale atomic/molecular massively parallel simulator) package [15].

The interactions between the atoms were modelled using a machine-learned 'Gaussian approximation potential' (GAP) [13], recently developed for the Ge-Sb-Te system by our group [12, 16]. This potential model exhibits an accuracy close to that of the underlying density-functional theory (DFT) training set of data, while its linear-scaling nature provides the opportunity to model very large structures and investigate slow cooling rates. The validity of the force field has been discussed in our previous work, demonstrating that the potential can be used for rigorous modelling of the short- and medium-range order structure of glassy GST225 [12]. In addition, the complex distribution of local environments that coexist inside the glass can be successfully reproduced; hence, glassy models generated with the GAP potential appear not to be missing any of the atomic environments present in glassy GST225 [14].

System-size effects were studied by generating four glassy structures, with a quench rate of 15 K ps⁻¹ and varying in size by nearly two orders of magnitude with respect to the total number of atoms in each periodic cell, namely, 315, 900, 7200 atoms and 24 300 atoms, while the cell size in each case was calculated from the experimental density (5.88 g cm⁻³) [17], resulting in values of 21.65, 30.72, 61.44 and 92.27 Å, respectively (figure 1). In order to study the effect of the cooling rate, 900-atom glassy models were generated with six different quench rates, varying by two orders of magnitude, viz. 100, 50, 15, 10, 5 and 1 K ps⁻¹.

Zoom In Zoom Out Reset image size

It is noted that, in this study, we have consistently quenched with a linear cooling rate from the melt in order to systematically study the structural changes upon cooling. However, glassy models can be generated also with step-wise quenches [18–24], in which case the choice of cooling schedule can further influence the glassy structure . Ultimately, one would expect the quench rate in real devices to decrease with temperature according to Fourier's law. Stochastic boundary conditions can capture this phenomenology [25, 26], and have been applied as a robust tool for simulating the resulting non-equilibrium dynamics of irradiated amorphous GST225 [27, 28], but they were not considered in this work for simplicity.

All the glassy structures were generated following a melt-and-quench approach with the same protocol. For each glassy sample, the atoms were initially placed in a pseudo-random arrangement inside a cubic simulation box, with imposed constraints to avoid unphysically small interatomic distances. The canonical ensemble (constant number of particles, volume and temperature, or NVT) was applied to keep the density of the simulated amorphous model close to the experimental value. The Nosé–Hoover thermostat, with a relaxation constant of 40 fs, was chosen to control the temperature fluctuations. A time-step of 1.0 fs was used to integrate the equations of motion during the MD simulations. For each glassy model, the initial configuration was heated at 3000 K with a 30 ps MD run to ensure that the system was melted and homogenised. The molten structure was subsequently cooled down to 1200 K and equilibrated using first the NVT and then the NVE (microcanonical) ensembles, with 30 ps and 10 ps MD runs, respectively, in order to obtain a well-equilibrated liquid structure at this temperature. Glassy structures were obtained by gradually quenching the systems down to 300 K with a linear cooling rate in the NVT ensemble. A further 30 ps NVT MD run was carried out at 300 K in order to equilibrate the structures in each case.

Finally, structural data were collected from an additional 10 ps NVE MD run at 300 K, and all the properties related to the analysis of the short- and medium-range structure of the generated GST225 glass models were obtained by averaging over frames extracted from this run in each case.

We note that, for all the MD simulations presented in this work, the density of the glass models was fixed to the experimentally reported value for amorphous GST225. We would like to highlight that performing simulations with the NVT ensemble allows one to obtain structures fitting a pre-defined (in our case experimentally determined) density. In general, fixing the simulation-cell volume to match the target-glass density throughout the melt-and-quench MD simulations is a robust approach, employed and well tested in prior modelling studies, in order to generate glassy models representative of the experimental structures [6, 7, 14, 18–24]. In addition, from the calculations presented in this work, we did not aim to predict the density of the glass. We are aware that this fixed-volume approach will result in a glass structure with some residual stress in the cell; nevertheless, we believe that this does not affect the quality of the calculations, since the effect of the cell relaxation is expected to be small on the structural properties of the simulated systems [24]. Consequently, as the aim of the present study was to study the influence of the system size and cooling rate on the generated MD models, it appears appropriate to focus on glass structures obtained in an NVT ensemble.

3.1. Configurational energy

The configurational energy (U) of glass models obtained from melt-and-quench MD simulations was shown to have a logarithmic dependence on the quench rate (γ) [29, 30], as given by the equation below:

$$\begin{equation}U = {U_0} + {a_1}\left[ {1 + log{{\left( \gamma \right)}^{{a_2}}}} \right]\end{equation} \tag{ 1 }$$

where U₀ is an asymptotic limit for the configurational energy and a₁ and a₂ are fitting parameters. The same type of data has, with a similarly good agreement, been fitted to a power law [29, 30], such as the one given in the following equation:

$$\begin{equation}U = {U_0} + {b_1}{\gamma ^{{b_2}}}.\end{equation} \tag{ 2 }$$

The calculated average configurational energy of the different GST225 models was estimated from the production runs at room temperature carried out at the end of each melt-quench simulation. The resulting changes in configurational energy with quench rate for the GAP-generated amorphous GST225 models are given in figure 2. Two models are clear outliers, whether equation (1) or equation (2) was used for fitting the change in configurational energy with the cooling rate. These correspond to the cooling-rate extremes considered here. The model that was quenched with the smallest rate of 1 K ps⁻¹ has a much smaller configurational energy compared to the other amorphous models, and also smaller than the asymptotic limit for a glassy configuration, suggesting that it may be below the smallest critical quench rate needed for vitrification. When it comes to the model quenched at a relatively fast rate of 100 K ps⁻¹, it was found to show a very similar configurational energy to the one obtained from the 50 K ps⁻¹ model. Fast cooling rates of 100 K ps⁻¹ were found to deviate from the logarithmic dependence of the configurational energy in a previous study of multi-component oxide glasses [21]. This, together with a wide variability in the configurational energy of the fast quenched model, could explain the deviation.

Zoom In Zoom Out Reset image size

When it comes to finite-size effects, models of 900 and 24 300 atoms quenched at 15 K ps⁻¹ had very similar configurational energies to each other, whereas models with 315 and 7200 atoms, quenched at the same rate, had a slightly

higher energy. This reveals a significant variability in the model energy but no clear size effect.

3.2. Structure factor

The short-range order of the generated amorphous models was characterized by calculating the static structure factor for each glass. The total structure factors were calculated from the inverse Fourier transform of the radial distribution functions (RDFs) calculated for the models. The neutron-diffraction structure factors of glassy GST225 melt-and-quench GAP models of different sizes are shown in figure 3. The same maximum wavevector, Q_max = 25 Å⁻¹, and a smoothing factor of σ = 0.125 were used in the calculation of the structure factor from the RDFs of all models. We compare the simulated structure factors with an experimental measurement [31], which is displayed as the black solid curve in figures 3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows that GST225 GAP models are not significantly influenced by finite-size effects, at least beyond 900 atoms. A small difference can be seen in the low wavevector region of the smallest 315-atom model, when compared to the rest; nevertheless, the main features are largely the same.

Prompted by recent work on the system-size dependence of the fragile-to-strong transition in bulk metallic glasses [32, 33], we have carried out an Ornstein–Zernike analysis of the reduced pair-distribution function, G(r), along the lines of [34], for the largest glassy GST225 model containing 24 300 atoms. Details of this analysis can be found in the supporting information. We find that the sizes of all models considered in this work exceed the longest correlation length (5.6 Å for Te-Te correlations) by a factor of more than three, in agreement with the observation that any finite-size effects are limited to the medium and extended-range order of the glass. It is worth noting that the first reports of crystallization of the glass consisted of models that already exceeded this correlation length, although not by a large enough margin to rule out finite-size effects [5].

By contrast, the structure factors of the 900-atom models, quenched at different rates, show a sharpening of all the peaks with decreasing cooling rate, especially below 5 K ps⁻¹ (figure 3). This sharpening of the peaks in the structure factor is an indication of structural ordering and possible partial crystallization for the slowest quench rate, while the appearance of shoulders of the different peaks suggest potential phase segregation within the structure.

In order to visualize the resulting structure, the 900-atom glassy GST225 melt-and-quench GAP model quenched at 1 K ps⁻¹ is shown in figure 5. On visual inspection, the model seems to have segregated into GeTe- and SbTe-rich regions, and shows partial crystallization associated with the formation of ordered 4-membered square rings.

Zoom In Zoom Out Reset image size

We note that a significant dependence of the low-Q features of the structure factor of amorphous silica on the cooling rate has been reported in a previous simulation study [29], highlighting modifications in the medium-range structure of the glass when the cooling rate is decreased to values between 1–5 K ps⁻¹.

3.3. Chemical order

In order to capture chemical order in the amorphous structure, the binary order parameter of Cargill and Spaepen has been calculated [35], from the perspective of the quasi-binary (A/B) classification which is widely used in the phase-change materials community [36]. This classification groups Ge and Sb ions together as type 'A' and leaves Te as type 'B', highlighting the prevalent pattern of A-B alternation which is observed in most phase-change materials [37]. The quasi-binary chemical order parameter of a given structure is given by:

$$\begin{equation}\eta _{AB}^0 = \dfrac{{\dfrac{{{n_{AB}}n}}{{{x_B}{n_A}{n_B}}} - 1}}{{\eta _{AB}^{max}}}\end{equation} \tag{ 3 }$$

where x_A and x_B are the atomic fractions of each (generic) atom type. The quantities n_A and n_B are the average coordination numbers of A and B atoms. The quantities n_AB and n_BA are the average number of B atoms around A, and the average number of A atoms around B, respectively. Finally, n is the total average coordination number of the structure and the parameter η^max is given by:

$$\begin{equation}\eta _{AB}^{max} = \left\{\!\!\!\!{\begin{array}{*{20}{c}} {\dfrac{{{x_B}{n_B}}}{{{x_A}{n_A}}},{x_B}{n_B} \le {x_A}{n_A}\,} \\ {} \\ {\dfrac{{{x_A}{n_A}}}{{{x_B}{n_B}}},{x_B}{n_B} \gt {x_A}{n_A}}. \end{array}\,} \right.\end{equation} \tag{ 4 }$$

A related order parameter, denoted by α_x (where x = Ge, Sb, Te) can be obtained in order to highlight the difference in quasi-binary chemical order of an amorphous structure from the perspective of each species. This parameter has been used in a previous simulation study [4] to describe chemical order in GST225 and is given by:

$$\begin{equation}{a_x} = \begin{cases} {\begin{array}{*{20}{c}} {\!\!\frac{1 - \dfrac{n_{xB}}{{{x_A}{x_B}\left( {{n_A} + {n_B}} \right)}}}{1 - \dfrac{{n_x}}{{x_A}{x_B}\left( {{n_A} + {n_B}} \right)}},\,x = {\rm Ge,Sb}\,} \\ \\[-6pt] {\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\frac{{1 - \dfrac{{{n_{xA}}}}{{{x_A}{x_B}\left( {{n_A} + {n_B}} \right)}}}}{{1 - \dfrac{{{n_x}}}{{{x_A}{x_B}\left( {{n_A} + {n_B}} \right)}}}},\,x = {\rm Te}\,.} \end{array}\,} \end{cases}\end{equation} \tag{ 5 }$$

The two approaches share a common interpretation, namely that a value of 1 signifies a fully ordered structure, such as the ideal cubic phase, a value of 0 signifies a random arrangement and finally, a value of −1 signifies a tendency for clustering, or phase separation. The amorphous GAP structures generated in this work have chemical order parameters between 0.7 and 0.95 in all cases, which are very close to those reported from DFT simulations of melt-quenched amorphous models of GST225 [2–7]. We note, though, that GAP glassy models are slightly overstructured, with an increased proportion of chemically ordered (AB) bonds and a lower fraction of tetrahedral Ge, when compared to similarly generated DFT models, as mentioned in our previous work [12]. Analogously to the case of the structure factors, the finite size of the amorphous models has no clear influence on the chemical order, as can be seen in figure 6. By contrast, slower cooling rates result in an increase of the chemical order, as seen in figure 7, with two regimes, below and above 5 K ps⁻¹ (or 5 × 10¹² K s⁻¹).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.4. Bond and dihedral angles

Moving on to real-space structural indicators, the bond- and dihedral-angle distributions were calculated. Given the significant amount of chemical order present in all of the amorphous models, only the bond angles made from quasi-binary ordered triplets (Te–Ge–Te and Te–Sb–Te), and dihedral angles from quasi-binary ordered quartets (Te–Ge–Te–Ge and Te–Sb–Te–Sb), were analysed.

The distributions of Te–Ge–Te and Te–Sb–Te angles for models of different sizes, as well as for models generated with different cooling rates, are shown in figure 8. The usual distribution centred just below 90° corresponding to angles from defective octahedral environments, and a smaller peak in the region below 180° due to the angles formed between axial bonds, is observed in both cases [38]. The reason for the distortion away from the ideal octahedral angles of 90° and 180° is the presence of lone-pair electrons which occupy a larger volume than bonded electrons, as predicted by VSEPR theory [39]. The lone pairs 'push' the bonded atoms closer together and reduce these bond angles. The effect is larger for Sb-centred angles, due to the larger Sb-centred angles due to the large volume occupied by Sb lone pairs when compared to Ge lone pairs.

Zoom In Zoom Out Reset image size

A small number of tetrahedral Ge environments have been identified in all of the models, corresponding to the density of angles around 109.5°. As previously reported, the fraction of tetrahedra in the GAP glassy models (6%–11%) is lower when compared to DFT studies, leading to overstructured models [12].

Overall, the angular distributions become slightly wider in models with a larger size, and narrower, and hence more ordered, as the cooling rate is reduced. The Te–Ge–Te bond-angle distribution becomes more ordered in response to the slower quench rate, when compared to the behaviour of the Te–Sb–Te bond-angle distribution, which remains relatively unchanged.

Compared to the bond-angle distribution functions, the dihedral angular distribution functions are more sensitive to both the finite-size effects and the cooling rate. They nevertheless tell a similar story, as topological ordering is observed from an increase in the peaks around ϕ = 90° and 180°, and a decrease of the intensity around ϕ = 45° and 135°, is observed in smaller models, in figures 9(a) and (b). Models obtained with a slower cooling rate also become more ordered, as shown in figures 9(c) and (d), and torsions about Ge–Te bonds respond more strongly to the ordering effects compared to those around Sb–Te bonds.

Zoom In Zoom Out Reset image size

The tetrahedral angular order parameter of Errington and Debenedetti [40], for an atom i, is given by:

$$\begin{equation}{q_4} = 1 - \frac{3}{8}\mathop \sum \limits_{j = 1}^{n - 1} \mathop \sum \limits_{k = j + 1}^n {\left[ {\cos \left( {{\vartheta _{jik}}} \right) + \frac{1}{3}} \right]^2}\,\end{equation} \tag{ 6 }$$

where the j and k indices run over nearest neighbours of atom i. This parameter is 0 for an ideal gas or an octahedral environment and 1 for a tetrahedral local environment. Other important values are ⅞ for a 3-coordinated pyramidal environment, ⅝ for a 4-fold coordinated defective octahedron (see-saw structure) and ⅓ for a 5-fold defective octahedron (square pyramidal environment). Figure 10 shows the change in the distribution of the q₄ parameter with the quench rate. From 100 K ps⁻¹ to 10 K ps⁻¹, the distribution moves towards more pyramidal and tetrahedral environments. Below 10 K ps⁻¹, however, the structure starts to become more octahedral with a dramatic change below 5 K ps⁻¹. For the 1 K ps⁻¹ model, we obtain a wide distribution corresponding largely to octahedral and defective octahedral environments.

Zoom In Zoom Out Reset image size

3.5. Rings

Shortest-path ring statistics were calculated using the Franzblau algorithm [41] as implemented in the QUIP code [42]. The resulting ring distributions are shown for models of different sizes in figure 11, while those of models quenched at different rates can be seen in figure 12. The 4-fold ring is the most probable one in all cases, and this is expected since it is an important structural motif in phase-change memory materials. ABAB square rings are also building blocks of the cubic crystalline phase and play an important role during the fast crystallization of the amorphous phase [4, 37].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Models of a smaller size, as well as those that have been quenched more slowly, show an increased amount of topological order. The 315-atom model has a slightly higher proportion of 4-fold shortest-path rings compared to larger models. Cooling rates have a more significant effect, especially below 5 K ps⁻¹, when the proportion of 4-, 6- and 8-fold shortest-path rings increases sharply, indicating structural ordering and, in particular, possible partial crystallization.

3.6. Voids

The GST225 material has intrinsic vacancies and exhibits cation disorder even as a cubic crystal. Metallic or degenerate semiconductor compositions that are fully-filled 'rocksalt-like' structures (such as GST224 or even GeTe) can lower their energy by forming cation vacancies, as well as by undergoing structural distortions, with the GST225 composition being close to the energetically optimal amount of cation vacancies [43]. The distribution of these vacancies within the crystalline lattice has been shown to be correlated with properties of the material, and to govern metal-insulator and cubic-to-hexagonal transitions in GST225 [44]. While vacancies cannot be clearly defined in the absence of a crystalline lattice, voids play a similar role in the glassy GST225 material where they were shown to enhance ionic diffusion [45]. Furthermore, different amorphous phases have also been identified as a function of pressure [46].

In order to further describe the changes in the topology of the glassy network, the statistics of non-interstitial voids within the structure were evaluated for the different glassy GST225 models. These voids have been identified on a fine mesh, with a maximal spacing of 0.3 Å, using a Voronoi construction with a cutoff of 2.8 Å [45]. The probability density of voids in the amorphous structure, assumed to be spherical with a radius r, was fitted using a log-normal distribution, as prompted by previous studies on vitreous silica [47, 48].

The resulting probability distributions of voids, as a function of the radius, for amorphous GST225 models of different sizes are shown in figure 13. The distribution of voids becomes sharper with increasing number of atoms, while the total void volume remains roughly the same. Therefore, it is more likely to find smaller voids in larger models.

Zoom In Zoom Out Reset image size

The void distribution as a function of the quench rate is shown in figure 14. For quench rates higher than 5 K ps⁻¹, there is a gradual sharpening of the distribution as the quench rate is lowered from 100 K ps⁻¹ down to about 10 K ps⁻¹. The distribution of voids for a quench rate of 5 K ps⁻¹ is broader than for 10 K ps⁻¹ but still sharper than for 50 or 100 K ps⁻¹. The void distribution changes dramatically below a quench rate of 5 K ps⁻¹. Given the significant density difference between the amorphous and crystalline phases of GST225, the much wider void-size distribution obtained for a quench rate of 1 K ps⁻¹ is another indicator of partial crystallization in this case. The cubic crystalline phase of GST225 is known to have a significant number of 'cation' vacancies, but it nevertheless has a higher density (6.30 g cm⁻³) than the corresponding amorphous phase (5.88 g cm⁻³). Given that the volume of the cells was kept constant during the quench, to maintain the experimental amorphous density, it is not surprising that large voids appear when the system is cooled slowly enough, as it tries to form a denser crystalline phase, instead of the glass.

Zoom In Zoom Out Reset image size

The above-mentioned results from a study of model-size and quench-rate effects on the structure of glassy models of GST225, simulated using the linear-scaling GAP potential, generally match expectations. For instance, it is under-standable that the contribution of larger rings to the ring-size distribution should be diminished in smaller-size models (especially in the smallest 315-atom model), and revealed in the largest models (figure 11). The ring-size distributions for all the simulated glassy models show a clear even–odd alternation, indicative of the medium-range order present in these glassy structures and resulting from the prevalence of 4-membered (ABAB-like) rings, a structural memory of the metastable rocksalt crystal structure of this material. Equally, the secondary maximum in the distribution for 12-rings, evident in figure 11, is significantly reduced with increasing model size and illustrates the changes in the medium-range order of the material.

However, there are much more pronounced quenching-rate effects, especially for the slowest quench rates. For instance, the numbers of 4- and 6-rings, and to a lesser extent of 8-rings, are greatly enhanced for the lowest quench rate, 1 K ps⁻¹ (figure 12). The dihedral-angle distribution also sharpens appreciably for the slowest quench rate, particularly for Te–Ge–Te–Ge dihedral angles (figure 9(c)). The Ge-centred (but not the Sb-centred, interestingly) bond-angle distribution shows a similar effect for near-90° bond angles (figure 8). Furthermore, the total structure factor exhibits fine structure, reminiscent of the presence of crystallinity, for the slowest quench rate GST225 model (figure 4). The second largest component of the Fourier transform (FFT) of the atomic density in the case of the 1 K ps⁻¹ model, commonly used as an order parameter for crystallization [5, 12], is low (0.25) when compared to crystalline models at a higher density (> 0.8), but is more than double that of the other glassy models at the experimental amorphous density (∼ 0.11). Being sensitive to the periodicity of the structure, the FFT order parameter might be low due to the formation of large voids in this model. The number of ABAB square rings, also used as an order parameter for the crystallization [12, 37], is nearly 70% that of an ideal cubic crystal. Considering the low configurational energy, the increased proportion of ABAB square rings, and the local densification of the structure, we believe the 1 K ps⁻¹ model to be partially crystalline.

The pattern in the chemical order and tetrahedral angular order parameter q₄ corroborate three distinct regimes in quenched GST225. For very fast quenches, faster than 100 K ps⁻¹, the structure resembles that of the frozen liquid. Below this value, we see significant AB-ordering and the q₄ parameter indicates the formation of pyramidal and tetrahedral environments. The familiar structure of glassy GST225, similar to that previously reported from DFT computational studies, is obtained around 10 K ps⁻¹. This ordering process, together with the formation of tetrahedral environments, is in agreement with the idea that an 'ideal' amorphous phase can be distinguished from both the liquid and the crystalline phases due to the formation of tetrahedral local environments. Finally, below 5 K ps⁻¹, the structure undergoes another dramatic ordering process and becomes more connected and almost entirely octahedral, resembling the metastable cubic crystal structure instead of the glass.

Thus, from these quench-rate-dependent simulation results, it can be concluded that, for this (GAP) interatomic potential for GST at least, the critical quench rate, below which liquid GST225 cannot be vitrified, lies in the range 1–5 K ps⁻¹ (1–5 × 10¹² K s⁻¹), and probably nearer the lower value than the upper one. This estimate is in very good agreement with recent DFT-based MD simulations of GST225, using the PBE exchange-correlation functional, where the critical quench rate was estimated to be 4.5 K ps⁻¹ [49]. However, it should be noted that these DFT-based simulations were for rather small model sizes (109 atoms, rather than the 900-atom models used in the corresponding present GAP-MD simulations), and so their estimate for the critical vitrification quench rate might be somewhat too high, since it is generally easier to crystallize small models than larger models (of a similar density), due to finite-size effects.

We report the results of a series of molecular-dynamics simulations of the canonical phase-change non-volatile memory material, GST225, made using a new, linear-scaling, machine-learned interatomic (GAP) potential. We have explored the effects of model size, varying the number of atoms in the simulation cells by nearly two orders of magnitude (up to 24 300 atoms), and of quench rate, again varied by two orders of magnitude, down to 1 K ps⁻¹ (10¹² K s⁻¹). It was found that, for the slowest quench rate, the models were partially crystalline rather than glassy, and so this value of quench rate gives an estimate for the critical minimum quench rate from the liquid, below which this phase-change material cannot be vitrified.

F C M acknowledges financial support from the UK Engineering and Physical Sciences Research Council (EPSRC) Centre for Doctoral Training in Computational Methods for Materials Science under grant number EP/L015552/1, and resources provided by the 'Cambridge Service for Data Driven Discovery' (CSD3, http://csd3.cam.ac.uk) system operated by the University of Cambridge Research Computing Service (http://www.hpc.cam.ac.uk), and funded by EPSRC Tier-2 capital grant number EP/P020259/1.

Via our membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202, EP/R029431), this work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk).

Data supporting this publication can be found in an online repository at https://doi.org/10.17863/CAM.49880.