Atomistic Origin of the Enhanced Crystallization Speed and n-Type Conductivity in Bi-doped Ge-Sb-Te Phase-Change Materials

Phase-change alloys are the functional materials at the heart of an emerging digital-storage technology. The GeTe-Sb 2 Te 3 pseudo-binary systems, in particular the composition Ge 2 Sb 2 Te 5 (GST), are one of a handful of materials which meet the unique requirements of a stable amorphous phase, rapid amorphous-to-crystalline phase transition, and signiﬁ cant contrasts in optical and electrical properties between material states. The properties of GST can be optimized by doping with p-block elements, of which Bi has interesting effects on the crystallization kinetics and electrical properties. A comprehen-sive simulational study of Bi-doped GST is carried out, looking at trends in behavior and properties as a function of dopant concentration. The results reveal how Bi integrates into the host matrix, and provide insight into its enhancement of the crystallization speed. A straightforward explanation is proposed for the reversal of the charge-carrier sign beyond a critical doping threshold. The effect of Bi on the optical properties of GST is also investigated. The microscopic insight from this study may assist in the future selection of dopants to optimize the phase-change properties of GST, and also of other PCMs, and the general methods employed in this work should be applicable to the study of related materials, for example, doped chalcogenide glasses.


Introduction
Phase-change materials (PCMs) are the functional materials at the heart of phase-change random-access memory (PCRAM), one of a handful of digital-storage technologies which are contenders to replace current Si-based fl ash memory in the future. [ 1 ] PC memory operation is based on the fast, reversible switching of the PCM between amorphous and (meta)stable crystalline states, with measurable differences in optical refl ectivity and/or electrical resistivity, thereby allowing binary information to be encoded into the material structure. PCRAM doping concentration of ≈10 at%, the dominant charge-carrier type of the material changes from the native p-type to n-type.
Aside from its effects on GST, Bi is also well known as a dopant for tuning the properties of other chalcogenide glasses, for example, the GeSe system [ 16,17 ] and its alloys with Sb, [ 18 ] and GeTe. [ 17,19 ] Similar p-to n-type transitions have been observed in these materials at critical dopant concentrations, [ 16,20 ] and the phenomenon has been connected to a percolation threshold being reached. [ 20 ] The effect of doping on electrical and optical properties has been linked to the large polarizability of Bi, [ 21 ] and also to the nature of the Bi-chalcogen bonding. [ 16 ] Bi-doped chalcogenide alloys are frequently used as luminescent materials, [ 22 ] for example, for solar spectral converters, [ 23 ] and thus their range of applications extends beyond PCMs, and atomiclevel insight into the property enhancements that Bi doping imparts is of general interest.
In recent years, ab initio molecular-dynamics (AIMD) simulation has become an important tool for studying PCMs. [24][25][26][27][28] AIMD melt-quench simulations are more or less routinely used to model the structure of amorphous PCMs, [ 25,29,45 ] and have also been successfully applied to understanding the unusually fast amorphous-to-crystalline phase transition in GST, [ 24,27,28 ] as well as how the crystallization kinetics can be manipulated to achieve sub-nanosecond switching. [ 3 ] Such studies have also been used to investigate the behavior of dopants in GST, and their effect on its properties. [ 9,10,30,31 ] In this work, we have performed a detailed, fi rst-principles simulational study of the effect of Bi doping on the structure, dynamics and electrical and optical properties of GST, and how these effects vary with dopant loading. We focus in particular on how the dopant integrates into the host structure, and on the possible origin of its enhancement of the crystallization speed and change in the dominant charge-carrier type.

First-Principles Calculations
MD simulations of the complete melt-quench-anneal PC cycle were carried out on fi ve 180-atom models with compositions of the form Ge 2 Sb 2x Te 5 Bi x , for x = 0, 0.15, 0.3, 0.45 and 0.6 (up to 6.67 at% Bi). The system with x = 0.15 is comparable to one of the three materials studied elsewhere, [ 13 ] while that with x = 0.3 can be compared to the composition investigated elsewhere. [ 11 ] Initial confi gurations of atoms in cubic simulation cells with periodic boundary conditions were fi rst randomized by heating at 3000 K for 20 ps. The models were then maintained in a liquid state at 1200 K for 100 ps, before being quenched to 300 K at a rate d T /d t = -15 K ps −1 , which invariably led to amorphization. These amorphous models were subsequently crystallized by thermal annealing at 600 K for 1 ns (doped models) or 1.5 ns (undoped model).
Our simulations were carried out, as in many other studies, [ 24,25 ] under constant-volume conditions (i.e., an NVT ensemble), since this matches typical device confi gurations in which the PCM is confi ned within a fi xed volume by dielectric and electrode materials. At a fi xed density of 6.11 g cm −3 , intermediate between the experimental densities of undoped amorphous and crystalline Ge 2 Sb 2 Te 5 , we recorded a maximum variation in the pressure in the simulation cell of 1.5 GPa during the simulations ( Figure S13, Supporting Information). We also found that the average pressure was reduced with increasing Bi content, suggesting that the equilibrium density of the doped systems would be lower than the undoped one. Nonetheless, these differences in pressure are small, and, based on previous work, are not expected to infl uence the behavior of the material (e.g., crystallization dynamics) to a large extent. [ 28 ] To investigate the structure and physical properties in more detail, amorphous and crystalline confi gurations of each system, taken from the end of the quench and annealing periods of the simulation, respectively, were energy-relaxed, with the cell volume and shape being allowed to vary. Further structural characterization and property calculations were then carried out on these optimized models.

Local Environments of Bi in Glassy Ge 2 Sb 2 Te 5
AIMD modelling has previously been shown to be a powerful tool for studying the local atomic structure of dopant atoms in amorphous host materials. [ 10,30,32 ] To investigate the structural environments of Bi in glassy GST, we analyzed the 30 Bi-centered complexes present in the four Bi-doped simulated amorphous models. On classifying the dopant atoms by coordination number, four distinct types of Bi complex were identifi ed, with between three and six neighbors ( Figure 1 a; see also Figure S1, Supporting Information). In order to obtain a full octet, early group V elements often adopt three-coordinate pyramidal geometries. Three-coordinate complexes were observed in these simulations, albeit infrequently, but the three-atom bond angles were invariably closer to 90°, rather than the pyramidal angle of ≈107°. This is indicative of a large p-character to the bonding, which is likely due to the fact that the Bi 6s orbitals are heavily contracted, and thus do not participate signifi cantly in forming covalent bonds. In support of this hypothesis, the average bond angle across all 40 Bi complexes was found to be 85.7° (Figure 1 b).
Four-coordinate confi gurations were similarly rare. The vast majority of the Bi atoms formed fi ve-and six-coordinate complexes, with the former being the most common. Interestingly, in complexes containing two Te atoms geometrically opposite one another, we observed a "Peierls-like" alternation in the bond lengths, with one of the two Bi-Te bonds consistently being shorter than the other by around 0.2 Å (Figure 1 c). Our defi nition of bonding in the present case is based on interatomic distance, whereas it has been shown in other work [ 34 ] that, particularly in amorphous materials, proximity does not necessarily imply chemical bonding. Indeed, from chargedensity isosurfaces through the equatorial plane of typical fi ve-and six-coordinate complexes (Figure 1 a), it can be clearly seen that the shorter bond has signifi cantly more charge density along the internuclear axis, whereas the longer one represents a much more loosely associated atom. This suggests that complexes with coordination numbers higher than three are formed by a mixture of covalent and dative-covalent bonding, such as through the interaction of a chalcogen lone pair with a vacant Bi d orbital. This is similar to the recently proposed account of the bonding structure of pyramidal Ge complexes in melt-quenched GST. [ 34 ] From the complexes shown in Figure 1 a, it is evident that Bi bonds preferentially to the chalcogen (Te). While this was overwhelmingly the case, a handful of complexes including Ge and Sb as nearest neighbors (i.e., containing "wrong bonds") were observed, with varying coordination numbers ( Figure 2 ). Plotting the average fraction, x , of the Bi coordination shells made up by Ge/Sb and Te atoms as a function of Bi concentration, it can be seen that Ge/Sb-Bi bonding becomes more prevalent as more dopant is added. This suggests that Bi-Te bonding is energetically preferred, but that other types of bonds are also possible, and will be most prevalent when the concentration of the dopant is such that its preference for higher coordination gives rise to competition with other host atoms for bonding to Te.
Comparing the calculated average Bader atomic charges [ 35 ] on the Bi atoms in the four kinds of BiTe n complexes, although there is considerable variation, an approximate linear increase in the charge with n can be seen ( Figure S2a, Supporting Information). This trend is less pronounced, although nonetheless still evident, if complexes containing Ge/Sb are included in the average ( Figure S2b, Supporting Information), a fi nding which suggests that the formation of wrong bonds can have a significant impact on the electron density around the dopant centers. This latter observation can be rationalized from the relative differences in electronegativity between Bi and Ge (+0.1), Bi and Sb (-0.3), and Bi and Te (-0.8), which suggest that bonds to Ge/ Sb should withdraw less electron density from Bi than bonds to Te.
Finally, it is an established issue with generalized-gradient approximation (GGA) functionals, such as PBE, that strongly correlated electrons in occupied d orbitals are often treated incorrectly. Since, as speculated here, the Bi d orbitals may be involved in the formation of dative bonds, we tested the effect of applying a Hubbard U correction to these orbitals on the structure of the Bi complexes. After optimizing the geometry of the four doped amorphous models with a 5 eV correction, we observed very little change to the average bond angles and lengths, with the short-long bond-length alternation pattern similarly being preserved ( Figure S14, Supporting Information). This suggests that these effects do not play a signifi cant role in determining the local geometry of the Bi atoms, and should correspondingly have little infl uence on the dynamics of the models.

Crystallization Dynamics
One of the key effects of Bi doping on chalcogenide PCMs is an enhancement of the crystallization speed. [ 11,12 ] While studies of crystallization kinetics with the small models (10-100 s of atoms) typically used in AIMD simulations on PCMs may exhibit artefacts from fi nite-size effects, previous work [ 28 ] has suggested that AIMD simulations of the crystallization of undoped GST, using models of the size employed in the present work, can show the microscopic physics of both The Bi and Te atoms are colored magenta and green, respectively. The charge-density isosurfaces to the right of the fi ve-and six-coordinate complexes illustrate the short/long bond alternation for pairs of geometrically opposite Te atoms. b) Distribution of Bi-centered three-atom bond angles within the 30 Bi complexes in the four doped amorphous models examined in the present work. The fi tted Gaussian curve has a mean of 85.7°, and a standard deviation of 9.4°. c) Distribution of the short (blue bars, solid black Gaussian) and long (red bars, dashed black Gaussian) bond lengths in pairs of geometrically opposite Te atoms in fourfold-and higher-coordinated Bi complexes. The means and standard deviations of the two fi tted Gaussian curves are 3.07 and 3.26 Å, and 0.07 and 0.1 Å, respectively. The snapshots and isosurfaces in (a) were prepared using the VMD software. [ 33 ] the nucleation and growth processes. In particular, these simulations estimated a critical nucleus size of ≈5-10 connected (Ge, Sb) 4 Te 4 cubes, that is, signifi cantly smaller than the total number of atoms in the simulation, and yielded an extrapolated crystal-growth rate in very good agreement with experimental work. [ 28,36 ] Based on this, and making the assumption that the Bi doping would not drastically alter the criticalnucleus size, it seems reasonable to investigate whether these simulations reproduce the enhanced crystallization speeds reported experimentally. To do so, we adapted the Fourieranalysis method, used by Hegedus and Elliott [ 24 ] to follow the evolution of structural order during crystallization, in order to provide a relative estimate of the crystal-growth rates in our models ( Figure 3 ).
In all fi ve simulations, crystallization is marked by a characteristic sharp increase in the maximum-intensity Fourier component, indicative of the rapid structural ordering occurring during crystal growth (Figure 3 a). The magnitude of the slope of these "growth edges" provides, in principle, an "ordering rate", which functions as an indirect measure of the crystal-growth velocity. To estimate this slope quantitatively, the smoothed Fourier curves were differentiated (Figure 3 b), and the most intense peaks in this fi rst-derivative space, corresponding to the growth edges, were identifi ed and fi tted to Gaussian functions (Figure 3 c).
For concentrations up to 5 at% Bi, the fi tted functions have higher maxima and become narrower with increasing dopant concentration, indicative of a faster peak ordering rate, and of crystallization completing in a shorter period of time. A comparison of the peak heights with Bi concentration (Figure 3 d) shows that the maximum growth rate roughly doubles over this concentration range. The fact that the 6.67 at% system does not follow this trend appears to be due to the model crystallizing via a "two step" process in this particular simulation; this phenomenon is evident as a "pre-edge" feature in several of the models, but is most prominent in the 6.67 at% doped model. While this could be an effect of the high dopant content, it is equally likely to be due simply to a manifestation of the stochastic variability inherent in crystallization, and the single simulation performed in this work prevents further conclusions being drawn as to this effect. In a similar vein, while the curves in Figure 3 a show that the onset of crystallization in the four Bi-doped systems occurred much faster than in the undoped model, this is almost certainly due to the stochastic nature of nucleation, since, for example, one of the GST models studied elsewhere [ 28 ] crystallized in less than 200 ps, despite being of the same size, and being quenched and annealed under the same conditions.
In spite of these caveats, the present simulations do seem to reproduce the trend in experimental growth-rate enhancement imparted by Bi doping. To investigate the mechanism proposed by Simpson et al., [ 14,15 ] we studied the effect of the dopant on the compressibility of the amorphous and crystalline models. To do this, the volume of each model was changed by ±500 Å 3 (corresponding to volume reductions and expansions of up to ≈10%) in eleven equal steps, with the atomic positions being scaled accordingly, and the total energy was obtained after ionic relaxation. The resulting energy/volume data were then fi tted to the Murnaghan equation of state [ 37 ] to obtain estimates of the bulk moduli ( Figure 4 ). In both phases, Bi-doping generally reduces the bulk modulus, that is, makes the material less resistant to compression, which implies, as suggested elsewhere, [ 15 ] that the addition of this dopant may lead to a decrease in the viscosity, thereby facilitating enhanced atomic diffusion.
To investigate this further, we also calculated the average atomic mean-square displacements (MSDs) in the fi ve models during the annealing stage of the simulation ( Figure S3, Supporting Information). While crystallization leads to a fl attening of the MSDs, at the beginning of annealing, before significant structural ordering has taken place, there is a time period during which the displacements are "liquid-like", and the gradients of the curves in this region, which are proportional to the average atomic-diffusion coeffi cients, can be used for a rough comparison of the diffusivities in the materials. There is a general increase in the gradient with Bi concentration, in support of the hypothesis that the Bi dopant may enhance atomic diffusion. Interestingly, the slopes of the MSDs of the 5 and 6.67 at% models appear to be larger than the trend for Adv. Funct. Mater. 2014, 24, 7291-7300 www.afm-journal.de www.MaterialsViews.com For this analysis, the radii of the coordination spheres of the Bi atoms were taken to be 3.5 Å. The color coding of the atoms is as follows: Ge -blue, Sb -red, Te -green, Bi -magenta. The images in (b) were prepared using the VMD software. [ 33 ] the undoped and the 1.67 and 3.33 at% models would predict, suggesting that the enhancement in diffusivity grows disproportionately with dopant concentration. However, as for the anomalously slow growth rate of the 6.67 at% model in Figure 3 d, it is quite possible that this is a peculiarity of this particular set of simulations, and repeat trajectories would need to be performed to verify the signifi cance, or otherwise, of this observation. Finally, it is of interest to investigate whether the presence of Bi may also lead to changes to the crystallization mechanism. To do this, we carried out an analysis similar to that reported in another study, [ 28 ] in which the time evolution of the number of atoms involved in structural units characteristic of the metastable cubic GST crystal structure, viz. fourfold rings, planes and cubes of atoms, was tracked during annealing ( Figure S3, Supporting Information). No signifi cant differences in behavior were observed, suggesting that the crystallization mechanism of the doped material is similar to that of undoped GST. However, it is perhaps noteworthy that, in the doped systems, increases in the number of rings appear to be more closely mirrored by corresponding changes in the number of atoms in planes, suggesting that Bi may help to stabilize planar structures. Also, in the present simulations, prior to the crystallization of the undoped system, there is a slow growth in the number of ordered structural units, which is absent in the doped models; however, this is likely simply to be an artefact of the undoped model crystallizing particularly slowly in this case.

Electronic Structure
The fact that the charge-carrier-type reversal highlighted elsewhere [ 12 ] occurs once a critical dopant concentration is Adv. Funct. Mater. 2014, 24, 7291-7300 www.afm-journal.de www.MaterialsViews.com Figure 3. Estimation of relative crystal-growth rates during the annealing of amorphous Bi-doped Ge 2 Sb 2 Te 5 . The time-evolution of the maximumintensity Fourier component (a) exhibits a sharp rise, indicative of rapid structural ordering, during crystal growth. The prominent "crystallization edges" appear as peaks after differentiation (b). Fitting the peaks to Gaussian functions (c) allows the gradient at the steepest point of the slope to be quantifi ed, which provides a relative measure of the maximum growth rate (d). In plots (a-c), the lines are color-coded as follows: undoped -blue, 1.67 at% -red, 3.33 at% -green, 5 at% -purple, 6.67 at% -orange. In plot (d), the dashed line is a least-squares linear fi t through the fi rst four data points, as a guide to the eye. surpassed suggests that it may be instructive to look at how the electronic structure of GST changes as Sb is systematically substituted by Bi. We therefore took the undoped amorphous and crystalline models and substituted 0-40 of the Sb atoms, at random, by Bi, generating a sequence of structures with compositions ranging from undoped Ge 2 Sb 2 Te 5 to "fully Bi-doped" Ge 2 Bi 2 Te 5 . These models were energy-relaxed, and the form of their electronic density of states (DOS) was examined.
A comparison of the DOS of the amorphous models ( Figure 5 ) shows that the effect of doping on the valence and conduction bands is somewhat subtle. The most prominent difference is a change in the shape of a feature ≈10 eV below the valence-band edge, composed mainly of s states from the four atomic species. However, an expansion of the DOS around the Fermi energy (Figure 5 b) reveals a more-or-less systematic shift in both the valence-and conduction-band edges with increasing Bi concentration.
We note in passing that the present simulations consistently yielded a fi nite density of states at the Fermi level in both phases of GST, likely indicating the absence of the experimentally-observed bandgaps, which is a common failure of DFT functionals. [ 38 ] Previous work has shown that, in the amorphous phase, this may be related to problems in the description of the structure (e.g., bond lengths and angles) arising from shortcomings of the DFT functionals. [ 29 ] This can be mitigated either by using more accurate functionals (e.g., a meta-GGA such as TPSS, [ 48 ] which makes use of the second derivative of the electron density in order to improve the description of nonlocal exchange and correlation effects), or by refi ning the simulated structure against experimental data. [ 29 ] However, while the use of improved functionals (and/or refi nement) yields more quantitative descriptions of atomic and electronic structure, different functionals lead to qualitatively similar behavior (e.g., as found in other work, [ 25,27,29 ] where GST is studied with three different functionals). Therefore, while we expect that the use of an improved functional would allow for a more quantitative description of the atomic and electronic structure of our models, it is likely that the qualitative trends discussed in the present work would be largely unaffected.
Analysis of the Sb and Bi partial DOS (PDOS; Figure 6 ) reveals the origin of the variation observed in the overall DOS. The Sb 5s and Bi 6s states both form twin peaks, but the lowerenergy part of the two is relatively much more intense for Bi than for Sb. This, together with slight changes to the Te 5s Adv. Funct. Mater. 2014, 24, 7291-7300 www.afm-journal.de www.MaterialsViews.com   Figure 5 . Subtle differences between the form of the Sb and Bi states account for the changes in the DOS of Ge 2 Sb 2 Te 5 as the Bi concentration is increased. As in Figure 5 , the color-coding of the bands is shown in the intensity scale to the right of the plots, with the undoped composition (0%) colored blue, and the fully doped one (100%) colored orange.
states ( Figure S5, Supporting Information), gives rise to the observed changes to the low-energy part of the valence band. For both Sb and Bi, p states make the major contributions to the band edges; comparing Figures 6 a and b, the Bi states have a noticeably narrower bandwidth than the corresponding Sb states, and they also peak closer to the Fermi energy. The differences are small, but nonetheless are suffi cient to cause the observed shift in the band edges. Aside from the changes to the Te s states, the form of the Ge and Te PDOS curves are practically unaffected by the Sb/Bi substitution ( Figure S5, Supporting Information), particularly within the vicinity of the Fermi energy. A lowering in energy of the valence-band edge might, in principle, make hole generation less facile, while an increase in the number of states near the bottom of the conduction band could likewise lead to improved electron conduction. If this effect plays a role in the carrier-type reversal, then the systematic changes in the electronic structure of the material with Bi concentration observed in these simulations provide a simple explanation as to why the charge-carrier reversal occurs beyond a threshold doping level.
The same shift of the band edges with respect to the Fermi energy, and the same differences in the Sb and Bi states, are also evident in an identical analysis carried out on the crystalline models ( Figure S6,S7, Supporting Information). This suggests that similar phenomena may occur in both material states, although due to the differences in electronic structure between them, it is quite possible that the threshold concentration for carrier-type reversal may also differ.
In addition to the shift in energy of the band edges observed in the DOS curves, another effect which could play a role in the reversal is changes in the Fermi energy due to defect formation, the energetics of which could be altered by doping. Vacancies on the cation sublattice in the crystalline phase, which are intrinsic to Ge 2 Sb 2 Te 5 , and are easily formed in related materials such as GeTe, [ 49 ] are known to promote p-type conduction. From experimental studies, Bi appears to substitute for Sb in the material, [ 11 ] as modelled in the present simulations, and so one would expect the anion-to-cation ratio largely to be preserved on doping. One possibility is that differences in the strength of Bi-Te and Sb-Te bonding could make the formation energies for cation vacancies higher in the doped material, reducing the concentration of vacancies at a given temperature and hence making the crystal "less p-type" than the undoped equivalent. However, additional calculations would be needed to investigate this further, and it is unclear how these effects, if signifi cant, would translate to the amorphous phase.
A third possibility is changes to carrier mobility through (de)localization of states at the band edges. To investigate this, we computed the inverse participation ratios (IPRs) of states around the Fermi energy for a subset of the amorphous and crystalline Bi-substitution models ( Figure S8, Supporting Information). We observed a slight increase in the IPR of states at the band edges with Bi concentration in both phases, with the increase being most pronounced in the crystalline phase. In this phase, the increase occurs more or less evenly in both the valence and conduction bands, whereas in the amorphous phase, the changes are much more pronounced in the conduction band. From this analysis, it might be tentatively concluded that an additional factor in the carrier-type reversal could be an increased delocalization of states at the band edges, particularly electronic states near the bottom of the conduction band. However, we note that the IPR calculation is likely to be infl uenced to a signifi cant extent by the same problems with the PBE functional which lead to the absence of a gap in the DOS at the Fermi level, and thus should be treated with caution.
The preceding analyses were performed by substituting Bi into the undoped material after quenching/annealing, and thus neglect any possible infl uence of the dopant atom on the structure through alteration of the kinetics of the melt-quenching/ crystallization processes. We therefore computed DOS curves for the models obtained from the MD simulations, and compared them to the substitution models, to see whether the same trends were observed ( Figure S9,S10, Supporting Information). Although we noted more variation among the DOS curves of the MD structures, the main features, viz. the changes in the s -like feature in the valence band, and the shift in the valenceand conduction-band edges, were clearly evident, in support of the fi ndings from the more systematic substitutional compositional study performed here.
Finally, as in Section 3, we investigated what, if any, effect a Hubbard U correction, applied to the Bi d orbitals, would have on the trends established from this analysis. We therefore re-optimized a subset of the substitution models with a 5 eV correction applied, and compared the total DOS and Bi PDOS curves to the corresponding ones from the uncorrected models ( Figure S15,S16, Supporting Information). We observed no signifi cant changes to the electronic structure, and the trends in the DOS with Bi concentration remained clearly evident.

Optical Properties
Finally, in order to investigate the effects of Bi doping on optical properties, we simulated the dielectric functions of the relaxed amorphous and crystalline models, from which we calculated the refractive indices, n ( λ ), extinction coeffi cients, k ( λ ), and refl ectance, R ( λ ), as a function of wavelength ( Figure S11,S12, Supporting Information). Using the defi nition of optical contrast given elsewhere, [ 11 ] we hence computed the optical contrast between the amorphous and crystalline models at the three technologically important optical-disc wavelengths of 405 nm (Blu-Ray), 650 nm (DVD) and 780 nm (CD) ( Figure 7 ).
The overall shape of the simulated n ( λ ) and k ( λ ) curves ( Figure S9, Supporting Information) compares well to the experimental measurements found elsewhere. [ 13 ] The refractive indices span a range of ≈1.75 to ≈5.5 over the 300-1050 nm range, which is comparable to the experimental measurements, although the largest calculated extinction coeffi cients are ≈50% higher than the experimental values. Overall, the agreement with these experiments is reasonable, considering that the simulations were performed only at the generalized-gradient approximation (GGA) level of theory, rather than with more advanced electronic-structure methods, for example, hybrid functionals such as PBE0 [ 39 ] or one of the screened-exchange HSE functionals. [ 40 ] The calculated refl ectance-contrast curves all show a small overall decreasing trend with Bi concentration, although the reduction in contrast is marginal. According to the calculations, the refl ectance contrast in GST is similar at all three wavelengths, being ≈5% higher at 650/780 nm than at 405 nm. Compared to the undoped material, the Bi-doped models show a maximum reduction of around 5-10%, suggesting that, as remarked elsewhere, [ 11 ] the doped materials should still be suitable for optical recording. The large variation between dopant concentrations observed in the literature [ 13,47 ] is not observed in our results, although many of the concentrations investigated in these studies, including those at which anomalous doping effects are observed, [ 46,47 ] are lower than the smallest concentration studied in this work.

Conclusions
By performing AIMD simulations on Bi-doped GST at various dopant concentrations, we have obtained valuable insight into the behavior of Bi within the host material, and into the microscopic mechanisms underpinning the property enhancements which it imparts.
Bi predominantly forms defective-octahedral BiTe n coordination complexes in amorphous GST, which appear to be formed through a combination of p-bonding and dative bonding; it is therefore natural that Bi is able to substitute for Sb in the cubic metastable rocksalt crystal structure. In an analysis of the crystallization dynamics in our models, we observed an enhanced crystallization speed, as is typically seen in Bi-doped PCMs experimentally. Our results suggest that this occurs due to an increase in free volume in the amorphous phase with dopant concentration, which lowers the viscosity and hence allows for easier atomic diffusion during the early stages of crystallization. There is also some evidence that Bi may stabilize planar structures, characteristic of the rocksalt structure, during crystallization, but additional simulations are needed to verify this.
A detailed analysis of the changes in the electronic density of states with Bi concentration suggests that the switch of the dominant charge-carrier type from the native p-type to n-type may result from the narrower bandwidth and lower relative energy of the Bi p states, compared to those of Sb, which leads to a shift of the valence-and conduction-band edges to lower energies with increasing Bi content. There is also some evidence that the dopant may cause a slight delocalization of states at the band edges, leading possibly to changes in carrier mobility. Finally, our models suggest that the optical contrast in the doped materials remains suitable for optical PC recording at technologically relevant optical-disc laser wavelengths.
We close with some general comments on the possible consequences that the small number of simulations conducted in this work might have on the generality of our fi ndings. The limited numbers of simulations (typically only one) and small model sizes employed in most current AIMD studies on PCMs makes them prone to artefacts from insuffi cient statistical sampling of rare events, particularly stochastic processes such as crystal nucleation. In the present work, the most prominent artefacts would likely manifest themselves in the amorphousphase structures-which are inherently random in natureand associated properties, and in the crystallization kinetics. In the former case, the approach adopted here, viz. running several simulations in parallel with slightly different dopant concentrations, would allow signifi cant anomalies to be identifi ed; the fact that most of the computed structural and electronic properties were observed to follow fairly clear trends suggests that qualitative conclusions here can be drawn with reasonable confi dence. On the other hand, crystallization is a rare event, and several anomalies in the crystallization kinetics are evident in the present set of results-in particular, the apparently anomalous relative growth rate of the 6.67 at% doped model, and the large differences in crystallization onset time between all fi ve systems. Again, however, conducting parallel simulations in this case allowed these anomalies to be highlighted, reinforcing the utility of this approach.
In summary, this study has provided a detailed atomic-level insight into how Bi dopant atoms interact with and infl uence the properties of a GST host, and we expect that these results may also be relevant to other Bi-doped telluride glasses, for example, GeTe. The methodology adopted in this work builds on the existing body of AIMD simulation work, including that on doped materials, to study systematically the effect of dopant concentration on material properties, which, as illustrated by the present results, can help more concretely to establish the effects of the dopant. We believe that a similar approach may be useful to study a variety of other materials, in particular other doped chalcogenide glasses. The insight obtained from this study should also be useful in facilitating the rational selection of dopants to optimize the properties of GST for specifi c memory-device applications in the future.

Computational Methods
DFT calculations were carried out using the VASP code. [ 41 ] The PBE exchange-correlation functional [ 42 ] was employed with projector augmented-wave (PAW) pseudopotentials [ 43 ] and an energy cut-off of Adv. Funct. Mater. 2014, 24, 7291-7300 www.afm-journal.de www.MaterialsViews.com Figure 7. Refl ectance contrast between the amorphous and crystalline phases of Bi-doped Ge 2 Sb 2 Te 5 , as a function of Bi concentration, at three technologically relevant laser wavelengths, viz. Blu-Ray (405 nm), DVD (650 nm) and CD (780 nm). The optical contrast was calculated according to the defi nition given in the literature [ 11 ] viz. O. C. = ( R cry -R am )/ R cry . For comparison, the experimental refl ectance contrasts measured in this reference are overlaid as fi lled data points.
175 eV for the plane-wave basis. This cut-off was increased by 1.3× for structural relaxations and property calculations. For the MD simulations and relaxations, the Brillouin zone was sampled at the gamma point, while for electronic-structure calculations and for dielectric-function simulations, a denser 2 × 2 × 2 gamma-centered Monkhorst-Pack [ 44 ] k -point mesh was used. For the MD simulations, a timestep of 5 fs was employed, with a velocity-rescaling algorithm being used to control the temperature. To improve the accuracy of the optical-property calculations, the number of electronic bands in the system was increased to 1728, roughly triple the default, and the number of grid points used to evaluate the DOS was increased to 10 5 . These parameters ensure that a large number of empty conduction-band states and a dense energy grid are used in the linear optics calculations, leading to a good sampling of wavelengths in the visible and near-infrared region.

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