Nucleus-driven crystallization of amorphous Ge 2 Sb 2 Te 5 : A density functional study

Early stages of nucleus-driven crystallization of the prototype phase change material Ge 2 Sb 2 Te 5 have been studied by density functional/molecular dynamics simulations for amorphous samples (460 and 648 atoms) at 500, 600, and 700 K. All systems assumed a ﬁxed cubic seed of 58 atoms and 6 vacancies. Crystallization occurs within 600 ps for the 460-atom system at 600 and 700 K, and signs of crystallization (nucleus growth, percolation) are present in the others. Crystallization is accompanied by an increase in the number of “ ABAB squares” ( A : Ge, Sb, B : Te), and atoms of all elements move signiﬁcantly. There is no evidence of cavity movement to the crystal-glass interface, as suggested recently, and the existence of Te-Te, Ge-Ge, Ge-Sb, and Sb-Sb (“wrong”) bonds is an inevitable consequence of rapid crystallization.


I. INTRODUCTION
Phase change (PC) materials play crucial roles in rewritable media and are promising candidates for wider applications in nonvolatile computer memory such as PC-random access memory (PC-RAM). 1,2 Their function is based on the rapid and reversible change between nanoscale amorphous (a) and crystalline (c) spots in a thin polycrystalline film, and the phases are characterized by different resistivities and optical properties. Two families of chalcogenide alloys dominate practical applications: (GeTe) 1−x (Sb 2 Te 3 ) x pseudobinary alloys, and doped alloys of Sb and Te near the eutectic composition Sb 70 Te 30 . Digital versatile disk (DVD)-RAM and Blu-ray Disc are examples of the first family, where Ge 2 Sb 2 Te 5 (GST-225, x = 1/3, referred to below as GST) is often taken as a prototype, and DVD-rewritable disks provide examples of the second. 3 Amorphous bits of materials in both groups are stable at room temperature for some decades, and the recrystallization of these bits can be accelerated greatly by laser irradiation or ohmic heating. The striking characteristic of PC materials is the extremely rapid (order of nanoseconds) crystallization of the amorphous phase. An understanding of this process requires a detailed knowledge of the structures of both ordered and disordered phases, but information about the amorphous phase, in particular, was lacking until very recently. The interplay between the theory and experiment has now identified "ABAB squares" (A: Ge, Sb; B: Te) as the dominant motif in amorphous GST and other members of this family. [4][5][6][7] This pattern is also prevalent in the metastable (rock salt) structure of GST, and it is plausible that the rapid amorphous-tocrystalline transition be viewed as a reorientation (nucleation) of disordered ABAB squares supported by the space provided by cavities. The crystallization of GST has been simulated in samples of less than 200 atoms with interesting results. 6,[8][9][10] The most striking results of the recent simulations (180 atoms, 600 K) were (a) the critical crystalline nucleus comprises 5-10 "ABAB cubes," (b) the process results in a cavity-free crystal, and all cavities segregate to the glass-crystal interface.
The role of finite-size effects in 180-atom samples is not yet clear.
We study here the crystallization of GST using density functional (DF) calculations combined with molecular dynamics (MD). A sample of a-GST with 460 atoms was studied at 500, 600, and 700 K, and a second sample of 648 atoms was simulated at 600 K. In all cases we used a fixed crystalline seed (58 atoms, 6 vacancies), and we do not discuss the onset of nucleation. However, our larger samples reduce finite-size effects, and we show the effect of choosing different annealing temperatures. Simulations of this scale (up to 648 atoms over 1 ns) are near the limit of present day DF/MD calculations. Still larger samples (4096 atoms) of GeTe have been studied using a classical force field derived from a neural networks representation of the potential energy surface derived from many DF calculations. 11 In Sec. II we outline the methods we use, and the results are presented in Sec. III. We discuss the results and their interpretation in Sec. IV.

A. DF/MD simulations
The simulations follow the pattern adopted in our earlier work. The CPMD program 12,13 is used with Born-Oppenheimer MD, a predictor-corrector algorithm, 14 and a time step of 3.0236 fs (125 a.u.). We employ an NVT ensemble (constant particle number N, volume V, temperature T) with cubic simulation cell, periodic boundary conditions, and a single point (k = 0) in the Brillouin zone; the PBEsol approximation 15 for the exchange-correlation energy; and norm-conserving and scalar-relativistic Troullier-Martins pseudopotentials. 16 We use a plane-wave basis with a kinetic energy cutoff of 20 Ry, and the temperature is controlled by a Nosé-Hoover thermostat (chain length 4, frequency 800 cm −1 ). 17,18 Simulations have been performed on two samples of a-GST (Fig. 1), with 460 and 648 atoms, respectively. We embedded in both a crystalline "seed" of 4 × 4 × 4 sites (13 Ge, 13 Sb, and 32 Te atoms, 6 vacancies) in a rock salt structure with lattice constant 3.0Å. The structure of the crystallite adopted the model of Yamada,19,20 which assumed that one sublattice of the rock salt structure comprises Te atoms, the other a random arrangement of Ge, Sb, and vacancies. We removed the atoms of the amorphous structure inside this volume, fixed the coordinates of the seed, and optimized the resulting structure.

460-atom structure at 500, 600, and 700 K
The initial geometry used the structure of a-GST that gave excellent agreement between DF simulations and experimental x-ray diffraction (XRD) and x-ray photoemission spectroscopy (XPS) measurements. 7 With fixed coordinates of the seed, we performed DF/MD simulations at 500, 600, and 700 K over 600 ps (500/600 K) and 350 ps (700 K). The path to crystallization at 600 K is shown in Figs. 1(a)-1(d). The amorphous and crystalline densities of GST differ, and the size of the cubic simulation cell was changed from 24.629Å (amorphous density of 0.0308 atoms/Å 3 ) to 24.060Å (crystalline density of 0.0330 atoms/Å 3 ) in five steps of 0.114Å.

648-atom structure at 600 K
The initial coordinates of the 648-atom sample were obtained from the a recent computer-aided deposition simulation. 21 The size of the cubic box was reduced in two steps from 27.736Å (amorphous density) to 27.150Å (density 0.0324 atoms/Å 3 ), which is incommensurate with the rock salt structure with this number of atoms. Complete crystallization was not expected and was not found at 600 K [see Fig. 1(f)].

B. Structure analysis
The analysis of the structures found focuses on several aspects. The radial (spherically averaged) properties are studied using pair distribution functions, and the bond orientational order (and the definition of "crystalline" atoms) are based on an order parameter appropriate for a rock salt structure. Cavities (vacancies) play an essential role in the phase change. Structural phase transitions of this nature are often discussed in terms of percolation, and we study the percolation of crystalline atoms as the structures evolve.

Pair distribution functions
The pair distribution function (PDF) g(r) is the spherically averaged distribution of interatomic vectors where ρ is the number density. In a system with components α, β,. . ., the PDF g αα , g αβ , . . . , can be calculated if ρ in Eq. (1) is replaced by ρ αβ = ρ √ c α c β , where c α and c β are the concentrations of elements α and β. The PDF for crystalline structures show sharp Bragg peaks that broaden progressively with increasing disorder. X-ray or neutron scattering measurements usually provide the total structure factor S(Q). If required, this and the partial structure factors S αβ (Q) can be found by Fourier transformation of the corresponding PDF. 4 Cavities play an essential role in determining the properties of GST and related systems, and α and/or β can be taken to be a cavity. The "center" of a cavity is determined as described below.
The local structure can also be characterized by the average coordination numbers, which are found by integrating g αβ to its first minimum R min N αβ = R min 0 dr 4πr 2 ρ αβ (r)g αβ (r). (2)

Bond orientational order
Crystallization is the ordering of a disordered system, and it is essential to quantify the directional order of a system of bonds and to define "crystalline" atoms. Angular correlations at an atomic scale can be achieved by projecting interatomic vectors r ij onto a basis of spherical harmonics Y lm ( r ij ), and the order parameter of Steinhardt, Nelson, and Ronchetti 22 has proved to be valuable in several contexts: 23 whereQ and N b (i) includes the atom i and its neighbors, thereby averaging over an atom and its first coordination shell and including intermediate range order. The single atom quantity where N(i) is the number of neighbors for atom i. For cubic structures, the first nonzero value ofQ is for = 4, and we define a "crystalline atom" as one for whichQ 4 0.6.

Percolation
The connection between statistical mechanics, particularly phase transitions, and percolation has a vast literature. 24 In the case of GST, simultaneous measurements of electrical resistivity and optical reflectivity showed a significant influence of percolation on electrical properties, but a negligible influence on optical properties. 25 We have analyzed percolation of the simulation trajectories by locating "crystalline" atoms as defined above and determined whether there was a continuous path of such atoms (maximum bond length 3.20Å) from an atom i to the same atom in a neighboring cell.

Cavities
Cavities (vacancies, voids) are defined and monitored as in our previous work. 4,5,26 We used a mesh of 0.057Å to determine the grid points (cavity domains) at least 2.8Å from all atoms. Cavities were determined for these domains using the Voronoi construction for disordered systems (the analog of the Wigner-Seitz cell in crystals). All points on a domain surface were used in the construction. The cavity center is the center of the largest sphere that can be placed inside the cavity without overlapping an atom. As noted above, the correlations of these centers with atoms and other cavities were monitored throughout, as were the distribution of cavity volumes and the total volume. The cavities have been visualized using the PYMOLDYN program. 27

III. RESULTS
The 460-atom simulations in 600 and 700 K temperatures showed full crystallization, and the former is visualized in In the simulation at 500 K and in the 648-atom simulation at 600 K [Figs. 1(e) and 1(f)] there was crystal growth around the seed, but crystallization was not complete within 600 ps.

A. Pair distribution functions, total energy, and crystalline atoms
The evolution of the partial PDF at 600 K is shown in Fig. 2. The crystalline seed gives rise to Bragg peaks at 3.00, 4.24, 5.20, 6.00, 6.71, and 7.35Å, which are enhanced by crystallization. The increase is most rapid between 280-300 and 480-500 ps at 600 K (Fig. 2). The partial PDF at 500 and 700 K are given in the Supplementary Information (SI). 28 The rapid increase at 700 K occurs between 80-100 and 180-200 ps. The total coordination numbers at the beginning and end of the 460-atom simulations (500, 600, 700 K) were calculated from the PDF for Ge, Sb, and Te atoms and cavities and are also given in the SI. 28 They show the expected change from amorphous to crystalline values for the samples that order and smaller changes otherwise. The PDF for cavities are discussed in detail in Sec. III E. Figure 3 shows the DF energy, fraction of crystalline atoms (via order parameterQ 4 ), and number of ABAB squares in all simulations. Crystallization is more rapid at 700 K than at other temperatures, with the crystalline fraction increasing most rapidly between 150-205 ps, but it is also complete within 600 ps at 600 K. The number of ABAB squares and the strength of the Bragg peaks correlate well with the number of crystalline atoms. At 600 K, the number of crystalline atoms increases by 40% between 330 and 435 ps, and it changes by 50% between 150 and 205 ps at 700 K. The more rapid crystallization at higher temperatures is consistent with the higher atomic mobility and with the results of recent ultrafast heating calorimetry measurements. 29 The fraction of crystalline atoms for each element is shown in Fig. 4. The fraction is lowest in Ge until near the end of crystallization at both 600 and 700 K, which is related to the coexistence of tetrahedral (minority) and octahedral Ge in the amorphous state. 4,30 As crystallization proceeds, the tetrahedral component of tetrahedral Ge becomes negligible, and the Ge fraction becomes comparable to the other two.
More information about the crystal growth mechanism can be found by studying the first two shells around the (4 × 4 × 4) seed: a 6 × 6 × 6 and an 8 × 8 × 8 cubes. The distance from the nearest fixed atom determines the shell of a particular atom.  to the second shell. The rapid crystallization (∼350-435 ps for 600 K, ∼150-200 ps for 700 K) is accompanied by an increase in the number of crystalline atoms in the first shell by ∼40 (600 K) and ∼60 (700 K), while the corresponding numbers in the next shells are ∼130 and ∼170. The first shells in these simulations are Ge-rich at 600 K and Sb-rich at 700 K.

B. Wrong bonds
"Wrong bonds" (Ge-Ge, Ge-Sb, Te-Te, and Sb-Sb) are those that do not occur in a widely accepted model of the metastable cubic structure, 19,20 and their existence is evident at 600 K in maxima in the partial PDF (Fig. 2) of the fully crystallized samples and in Fig. 6. The numbers decrease during the simulations, most rapidly during fast crystallization, and are largest in the systems (460 atom at 500 K, 648 atom at 600 K) that do not crystallize completely. Nevertheless, they are present throughout all simulations. The mobility and the ability to find energetically favorable configurations are higher at 700 K than at 600 K, resulting in fewer wrong bonds in the fully crystallized structure. The concentration of Te causes the large number of Te-Te bonds, and the Ge-Ge wrong  bond maxima are weak. There are only four Ge-Ge bonds near the end, but the first peak shows an unusually short bond less than 2.5Å. The closest Ge-Ge atom pair is adjacent to two cavities and separated by only 2.37Å (Fig. 7), and they have tetrahedral local configurations. This is exceptional for Ge, where the other atoms form the (defective) octahedral rock salt structure. There are more Ge-Sb bonds than Sb-Sb bonds after crystallization at 600 K, while the reverse is true at 700 K.

C. Percolation and crystal growth
Percolation of the crystalline cluster was analyzed from the trajectories, and the results are shown in Fig. 8. The simulations began with the same structure, but no preferred direction for crystal growth was found. Percolation occurred at 600 K in the order z, y, and x directions, and x, z, and y directions for 700 K. Percolation is present in the 500 K simulation in two directions, although the system did not crystallize fully (as in the 648-atom simulation). Percolation can be achieved with as few as ∼20% of the atoms being "crystalline," and its onset occurs well before the critical stage of crystallization. This is illustrated in Figs. 1(b) and 1(c) and in Supplementary Fig. 4. 28 In the last figure, percolation results in a network of periodic images of the simulation cell prior to rapid crystallization (∼ 350-435 ps) in the 460-atom simulation at 600 K.

D. Atomic motion, mean square displacement
Crystallization is also reflected in the mean-square displacement (MSD, Fig. 9) of the atomic coordinates from the starting structure. Changes in MSD slow down and eventually stop as crystallization occurs. Sb is slightly more mobile than Ge at 500 K until ∼500 ps, where the MSD of Ge increases rapidly. At 600 K, Ge has a higher MSD then Sb, while the reverse is true at 700 K. We note that Sb has a higher mobility in the liquid. 4 The evolution of σ (MSD)/MSD, where σ (MSD) is the standard deviation of the MSD of individual atoms, shows that there is a wide range of atomic mobilities in individual elements (e.g., for Ge atoms at 600 K and Sb atoms at 700 K). On the other hand, Sb atoms at 600 K and Ge atoms at 700 K have less variance than Te atoms. In general, Ge and Sb atoms are more mobile than Te, atoms crystallizing in the first shell outside the seed move less than atoms crystallizing in the second shell, and the crystallite grows faster in the direction parallel to the crystal lattice alignment of the seed.
The evolution of σ (MSD)/MSD (Fig. 9) shows that some atoms move very little, while others diffuse far during the simulations. The motion of individual atoms is an interesting aspect of the crystallization process, and examples of mobile Ge/Sb with almost stationary Te, and the reverse, can be found. Two particular cases are discussed in the SI. 28 Figure 10 shows the evolution of the total cavity volume in the 460-atom simulations at 600 and 700 K. The slight downward trend prior to crystallization is associated with the smaller cell, but the variation for constant cells (e.g., 0-100 ps, 200-300 ps, etc.) is not large at 600 K. The total cavity volume is rather stable during the rapid stage of crystallization, but there is a small increase towards the end (at 600 K from 400-450 ps, at 700 K from 150-200 ps). On the other hand, the overall volume reduction from the amorphous to crystalline state is significant, of the order of 1/3, which is consistent with the 7% higher density and smaller fraction of cavities/vacancies (10% of total volume) in c-GST. 7 The radial distribution of cavities from the center of the seed is shown in Fig. 11. The first two maxima arise from the seed and do not fluctuate. The finite size of the simulation cell means that the distribution tails off after ∼12Å, but there are cavities at all shorter distances down to the cutoff of the cavity domains (2.8Å). The cavities show definite order after crystallization  at 600 K [ Fig. 11(a)], but order is less pronounced at 700 K [ Fig. 11(b)].

E. Cavities
The size distribution of cavities is shown in Fig. 12 for the crystallized structures at 600 and 700 K. The distributions are similar, with periodic peaks at multiples of ∼35Å 3 , approximately the size of a single vacancy in the GST rock salt lattice with the above definition of a cavity. The 700 K distribution resembles a liquid before crystallization. The unusual peak at ∼80Å 3 before crystallization at 600 K is due in part to the merged (multi)cavities in the fixed seed and its interface.
The PDF involving cavities (Fig. 13) provide us with important information. There is a clear tendency for Ge and Sb atoms to move during crystallization away from cavities, although the final (rock salt) structure shows nearest-neighbor Sb-cavity "wrong bonds." Te atoms make up more than half of the total, and the number of Te-cavity bonds remains large throughout the simulation (note the different scales). The crystalline PDF show prominent peaks for cavities with Ge, Sb, and "cav" at distances corresponding to the opposite edges of ABAB squares, and the maximum at 5.2Å in Fig. 13(c) corresponds to Te atoms and cavities at opposite corners of ABAB cubes. There are peaks in the Ge-cav, Sb-cav, and cav-cav PDF near 6.0Å, twice the AB distance in the rock salt structure. The overall picture of crystallization is the motion of Ge and Sb atoms away from cavities to occupy sites on one sublattice of the rock salt structure. The phase transition is, however, so rapid that "wrong bonds" are inevitable.

F. Electronic densities of states
The projected electronic densities of states (DOS) (Supplementary Fig. 7) show small changes upon crystallization. There is a small (<0.5 eV) semiconducting gap at the Fermi energy with localization at the band edges, the DOS at the upper end of the Te s band increases significantly, and the gap between the two s bands at −10.5 eV increases slightly. The gap between the s and p bands around −6 eV decreases towards zero, which is typical for the electronic band structure of chalcogenides with the rock salt structure. 31

IV. DISCUSSION
The present series of calculations has focused on the crystallization of 460-and 648-atom samples of amorphous Ge 2 Sb 2 Te 5 (GST) from a fixed crystalline nucleus (58 atoms, 144113-7 6 vacancies) up to 600 ps. The 460-atom sample was simulated at 500, 600, and 700 K, the 648-atom sample at 600 K. The onset of crystallization from a completely disordered sample will be addressed in future work, but the present results go beyond earlier work in several ways. Our simulations are the largest performed to date for crystallization in GST and have studied the process at three temperatures. We focus on the changes in the numbers of "ABAB squares" (A: Ge, Sb, B: Te) and "wrong bonds," and the difference in densities between the amorphous and crystalline phases (0.0308 and 0.0330 atoms/Å 3 , respectively) is taken into account during crystallization by adjustment of the cell size. Crystallization is defined in terms of bond orientational order parameters, and we study percolation in all samples.
Crystal growth of the nucleus is evident in all samples, but the process is complete only in the 460-atom sample at 600 K (∼400 ps) and 700 K (∼200 ps). The 648-atom sample is larger, and the final unit cell was incommensurate with the crystalline structure. Although these features appear to hinder the process, the crystallite also grows in this simulation. Each process involves a preliminary stage where the nucleus grows gradually and forms percolating connections (narrow necks) with replicas in the neighboring cells. For example, Fig. 1(c) and Supplementary Fig. 4 show that the system has paths with many crystalline atoms going through the cell boundaries in the x, y, and z directions before the rapid stage of crystallization occurs. Changes in the atomic structure prior to crystallization (thermal restructuring, pre-annealing) were found in the recent 180-atom DF/MD simulations of GST at 420 K and used to explain the more rapid crystallization observed experimentally. 10 Our simulation at 500 K does not crystallize fully, but the number of crystalline atoms increases, and percolation with the neighboring replicas begins.
Growth is clearly favored by the mobility brought by higher temperatures. Recent ultrafast differential scanning calorimetry measurements by Orava et al. 29 showed that the highest crystal growth velocity is achieved at 670 K, which is consistent with the faster process rate in our simulations at 700 K. The diffusion constant of Sb is highest in liquid GST, 4 and Sb is important in processes at higher temperatures. The definite correlation between the number of ABAB squares and the number of crystalline atoms supports early observations that these structural units were crucial to the speed of crystallization. 4,6 However, as shown in Supplementary Fig. 6, these units can break and reform during crystallization. 9 An unconstrained MD simulation of the original seeded structure at 600 K resulted in a reduction of ABAB squares by 20-30%. These squares are now scattered throughout the sample and do not form a crystallite. The diffusive nature of the atomic motion is an important aspect of crystallization in GST, and atoms can move significant distances during the transition. This applies to all elements, including Te ( Supplementary Fig. 5).
The existence of cavities in both amorphous and crystalline phases is characteristic of materials in the (GeTe) 1−x (Sb 2 Te 3 ) x family, 32 and it is widely assumed that they play an essential role in the rapid phase changes that can occur. Recent simulations of crystallization in GST indicated that cavity diffusion to the crystal/glass interface, followed by Ge/Sb diffusion to these sites, aided the formation of cubic, cavityfree crystallites. 8 We have not performed smaller simulations to check whether this is a consequence of using a 180-atom sample, but there is no evidence of this effect in our 460-and 648-atom simulations. The cavity ordering is particularly clear in the radial distribution function (Fig. 11).
"Wrong bonds" are those that are absent in a widely accepted model of the metastable structure of GST: a rock salt (Fm3m) structure with one sublattice containing Te atoms, the other a random mixture of Ge atoms (40%), Sb atoms (40%), and cavities (20%). 19,20 There may be substantial displacements from the ideal crystallographic positions, particularly for Ge (Ref. 20). The presence of wrong bonds, particularly Te-Te, in our crystallized samples is particularly striking. X-ray diffraction of GST films showed that exchanging all Sb atoms (atomic number 51) with Te (atomic number 52) in this model did not diminish the R factor significantly, 33 confirming that XRD measurements alone cannot distinguish between these two elements.
The model of a perfect Te sublattice is accepted without comment by many authors, 2 but the arrangement of Ge, Sb, and cavities on the other sublattice has often been discussed. Transmission electron microscopy measurements of thin deposited GST films indicated a possible ordering of Ge and Sb atoms, 34 and this was supported by DF calculations that showed ordering of cavities perpendicular to the [111] direction in the rock salt structure and no intermixing of Ge and Sb atoms in these planes. 35 Other DF calculations, however, found intermixing of Ge and Sb and ordering of the cavities perpendicular to the [111] planes, and identified a Te atom surrounded by 3 Ge and 3 Sb atoms as energetically favorable. 36 A study of the stoichiometric defects involved changes in the Ge/Sb sublattice, 37 and recent x-ray absorption and scattering measurements and DF simulations suggest that the displacements of many Ge atoms led to considerable bonding disorder. 38 Our graphical software 27 has allowed us to follow all simulations in detail, and there are indeed clear indications that the total energy favors Te occupancy of one sublattice of the rock salt structure. Nevertheless, the very short crystallization times mean that other contributions to the free energy, particularly configurational entropy, play crucial roles. Vibration frequencies are low in GST (typically 100 cm −1 or 3 THz), 39 so atoms vibrate no more than a few thousand times if less than 1 ns is available. The dramatic effect of rapid annealing on the observable structures is provided by studies of small sulfur cluster anions S − n (Ref. 40). Although ring structures are more stable than chains, 41 rapid cooling of the S − n vapor led to only chain structures for n = 6-9. Short annealing times allow only small changes from the initial structures and favor chains, which are not constrained to close and are much more numerous than rings. Nanosecond cooling in GST will similarly lead to structures with "wrong bonds," which are much more common than structures with a perfect Te sublattice.