First‐Principles Calculation of Transport and Thermoelectric Coefficients in Liquid Ge2Sb2Te5

Thermoelectric effects play an important role in the programming operations of phase‐change memories. The electrothermal modeling of the device thus requires the knowledge of transport coefficients such as the electrical and thermal conductivities and the Seebeck coefficient at different temperatures. Herein, these transport coefficients of liquid Ge2Sb2Te5 at a few temperatures above melting are calculated by means of density functional molecular dynamics simulations.


Introduction
Phase-change materials are exploited in key enabling technologies for nonvolatile electronic memories, [1][2][3][4][5][6] neuromorphic computing, [7,8] and a variety of photonic applications including computing devices and switches. [9,10] In phase-change memories (PCM) the two logical states of the memory are encoded in the amorphous and crystalline phases of a chalcogenide phase-change alloy. The two states of the active material feature a difference in electrical resistivity by about three orders of magnitude, which allows discriminating the two states of the memory by a measurement of the resistance at low bias. [1,3] The application of current pulses at higher bias induces a fast and reversible transformation between the amorphous and crystalline phases by Joule heating. In the set process, the amorphous phase recrystallizes, while the reset consists of heating the crystal above the melting temperature with subsequent amorphization due to fast cooling.
Thermoelectric effects have been shown to play an important role in the programming operations of the device due to the presence of large electric fields and thermal gradients. [11][12][13][14] Numerical simulations based on finite-element methods have been performed to reproduce the electrical characteristics of PCMs based on the prototypical Ge 2 Sb 2 Te 5 (GST) phase-change compound in the set and reset regimes. [13,14] This modeling requires as input parameters the transport and thermoelectric coefficients of the material as a function of temperature in the different phases such as the electrical and thermal conductivity and the Seebeck coefficient. These data are available from experiments for the crystalline and amorphous phases, while only partial information is available for the liquid. The electrical conductivity σ of GST has been measured in the temperature range 930-990 K, [15] while the thermal conductivity κ is typically estimated from electrical conductivity and the Wiedemann-Franz (WF) law, [13] although deviation of about 15% in the thermal conductivity estimated in this way has been reported in liquid Sb 2 Te 3 . [16] No experimental information is available instead for the Seebeck coefficient S in liquid GST, which is then estimated from a linear extrapolation with temperature of the values measured in the crystalline phase. [17]

Results
In this respect, first-principles simulations can provide information on the transport and thermoelectric coefficients of the liquid phase needed for the thermoelectrical modeling of the device and not available from experiments. To this aim, in this paper, we report on first-principles simulations based on density functional theory (DFT) of liquid GST that yield estimates for σ, κ, and S at different temperatures above the experimental melting temperature T m ¼ 858-877 K. [18,19] The molecular dynamics (MD) simulations have been performed within the same framework that we used very recently to study the metal-semiconductor (M-SC) transition in supercooled liquid GST with 270-atom supercells. [20] There is evidence in literature that the inclusion of van der Waals (vdW) interactions is important to properly describe the structural properties of liquid GeTe and Ge 2 Sb 2 Te 5 . [21,22] Therefore, we used both the Perdew-Burke-Ernzherof (PBE) [23] exchange and correlation functional without vdW interaction and the functional proposed by Vydrov and van Voorhis [24] and revised by Sabatini et al (rVV10) [25] that includes vdW corrections as implemented in the CP2k code. [26,27] Finally, as it will be discussed later on, we also assessed the dependence of the results on the use of a hybrid functional (HSE06), [28] which is known to better reproduce the electronic bandgap in the crystalline phases of GST. [29] We performed simulations at three temperatures 900, 1000, and 1100 K above T m . After equilibration for 5 ps, we computed average properties on the subsequent 6 ps, as detailed later.
As we are interested in providing data for the electrothermal modeling of the device, the simulations have been performed at constant volume, which is the situation experienced by the active material in the device due to the confinement by the surrounding dielectrics and electrodes. Liquid GST was modeled in an orthorhombic supercell containing 270 atoms with edges of 21.971, 21.971, and 18.643 Å corresponding to a density of 5.683 g cm À3 (0.0300 atoms Å À3 ), which is close to the experimental value in the temperature range of interest (from 0.0307 atoms Å À3 at 893 K to 0.0303 atoms Å À3 at 1093 K). [22] This density is actually equal to the theoretical equilibrium density at 1100 K obtained from constant temperature-constant pressure (NPT) simulations with the rVV10 functional. We note that the volume per atom at normal pressure increases by about 1.5% from 893 to 1093 K according to the experimental data in the study by Schumacher et al. [22] This small change in density has a very tiny effect on the electronic structure, leading to a change in the transport coefficient which is expected to be smaller than the uncertainties due to the other limitations of the DFT calculations described later.
The partial-pair correlation functions are reported in Figure S1-S2 in the Supporting Information (SI) for PBE and rVV10 simulations at different temperatures. The partial and total coordination numbers are given in Table S1 in Supporting Information. The structural properties of the liquid are consistent with previous works. [22,30] The inclusion of vdW interaction leads to an overall sharpening of the first coordination shell with a decrease in the coordination number at the lower temperatures and a lower fraction of "wrong" bonds (Ge-Ge, Ge-Sb, Sb-Sb) than the PBE results. This effect is much more pronounced in the previous simulations by Schumacher et al., [22] where the nonlocal vdW-DF2 [31] functional, which also includes vdW interactions, is used. In the vdW-DF2 functional, the local exchange and correlation energy lacks gradient corrections which are instead included in the rVV10 functional. The sharper coordination shell might then be partially due to the use of the LDA exchange and correlation energy in the vdW-DF2 functional. Partial pair correlation functions obtained at 925 and 1024 K from vdW-DF2 simulations in the study by Schumacher et al., [22] are compared with our PBE and rVV10 results in Figure S3-S4 in the Supporting Information. The differences in the structural properties between the liquid models generated with the PBE and rVV10 functionals have been also discussed in depth in our previous work, albeit at lower temperatures, to which we refer to for further details. [20] Fourier transforming the partial pair correlation functions yields the partial structure factors which are neutron weighted to obtain the neutron total structure factor shown in Figure S5 in Supporting Information. The rVV10 results are in better agreement with experiments than the PBE ones, although they are still worse than previous vdW-DF2 data.
The GST alloy is usually considered a semiconductor also in liquid phase because its conductivity increases with temperature. However, this behavior is not due to an increase in the thermally excited carriers across a pseudogap, but to the increase in the electronic density of states (DOS) at the Fermi level (E F ) due to delocalized states which progressively fill the pseudogap by increasing temperature. This view was discussed in our previous work on the M-SC transition in supercooled GST. [20] The electronic DOS shown in Figure 1 confirms the progressive filling of the pseudogap at E F by increasing temperature also above T m in agreement with previous works. [22] To quantify the localization properties of individual Kohn-Sham (KS) states, we computed the inverse participation ratio (IPR) which is defined for the i-th KS state by P j c 4 ij =ð P j c 2 ij Þ 2 , where j runs over the Gaussian-type orbitals (GTO) of the basis set and c ij are the expansion coefficients of the i-th KS state in GTOs. The IPR takes values varying from 1/N for a completely delocalized electron, where N is the number of atomic-like orbitals in the basis set of the whole supercell, to one for an electron completely localized on a single atomic-like orbital. The IPR superimposed to the DOS in Figure 1 shows that the KS states at E F are indeed delocalized in the temperature range considered here. As shown in our previous work, the IPR at E F increases by approaching the M-SC transition below T m . [20] As the IPR might be affected by finite size effects, we repeated the calculation with a larger supercell 31.072 Å Â 31.072 Å Â 37.286 Å containing 1080 atoms. The model was equilibrated for 5 ps at 1000 K. The corresponding DOS and IPR reported in Figure 2 confirm the absence of a mobility gap as the states at the Fermi level are delocalized in the larger cell as well.
As the states at E F are delocalized at the conditions considered here, the Kubo-Greenwood formula can be consistently used to compute the kinetic coefficients ℒ i,j . For an isotropic system these are given by [32,33] where ε c,k and ε v,k refer to the KS energies of bands at the N k kpoints in the supercell Brillouin zone (BZ), V o is the supercell volume, and f ε is the Fermi-Dirac distribution function at the given temperature T with chemical potential μ. The kinetic coefficients have been computed using the Quantum-Espresso [34] (QE) suite of programs to properly integrate the BZ for configurations extracted from the MD trajectories generated with the CP2k code.
The transport coefficients were obtained in turn from the ω ! 0 limit of the kinetic coefficients ℒ i,j ðωÞ as [35] σ ¼ ℒ 1,1 The kinetic coefficients were computed by substituting the δfunction in Equation (1) with a Gaussian function, as implemented in the KGEC code. [36] The transport coefficients from Equation (1) were averaged over several independent configurations in the liquid at each temperature.
We first checked the convergence of the transport coefficients with respect to number of k-points in the BZ by averaging over nine configurations at T ¼ 1000 K and using a 2 Â 2 Â 2, 3 Â 3 Â 3, 4 Â 4 Â 4, or 5 Â 5 Â 5 uniform meshes which www.advancedsciencenews.com www.pss-rapid.com correspond to 4, 14, 32, and 63 k-points in the whole BZ (calculations are actually performed over half of the points because of time-reversal symmetry). In Figure S6 in Supporting Information, we report the values of σ, κ, and S for the different meshes and different values of δ¼ ffiffi ffi 2 p σ g , where σ g is the variance of the Gaussian. Some of the transport coefficients are not fully converged even with the finer 5 Â 5 Â 5 mesh (63 k-points) for δ ¼ 0.04 eV. To reduce the computational load, we therefore used broadenings larger than 0.06 eV with the 4 Â 4 Â 4 mesh.
The resulting values for σ, κ, and S are reported as a function of the broadening δ in Figure 3 for the simulation at 1000 K with rVV10. The linear extrapolation to zero broadening provides our estimate for the transport coefficients. The final data in Figure 3 actually refer to the values obtained by averaging over ten configurations (separated by 0.5 ps each) according to the analysis of the convergence of the transport coefficient on the number of MD configurations summarized in Figure S7-S9 in the Supporting Information. The resulting values for σ, κ, and S are shown in Table 1 for rVV10 and PBE simulations. The average value and the mean square deviation in the energy and temperature of the trajectories used to compute the average transport coefficients are given in Table S2 in Supporting Information.
The electrical conductivity increases with temperature. As already stated earlier, although this behavior has been often associated in literature with a semiconducting character, this is not due to an increase in the density of thermally excited carriers as it  www.advancedsciencenews.com www.pss-rapid.com is the case for a normal semiconductor but to an increase with temperature of the DOS at the Fermi level. [37] The value of σ is sizably lower with the rVV10 than with the PBE functional at the lowest temperature, while it is similar for both functionals at the higher temperatures. This can be tracked back to the fact that the M-SC transition is only slightly below T m with the rVV10 functional, while it is about 150 K lower with the PBE functional, as shown in our previous work. [20] As a result, the pseudogap at E F is still deep at 900 K with the rVV10 functional. Sufficiently far from the M-SC transition, the two functionals provide very similar electrical conductivity, which is, however, sizably larger than the experimental value for Ge 1.6 Sb 2 Te 5 of 2433 (Ωcm) À1 at 923 K and of 2777 (Ωcm) À1 at 983 K. [15] A value closer to the experiments of 2169 (Ωcm) À1 was instead obtained for σ in our previous work, [38] in which we used the HSE06 hybrid functional with PBE trajectories at 1060 K. In this latter work, however, the BZ integration was restricted to the Γ-point of a 270-atom supercell.
To assess whether the discrepancy with our previous result on σ arises from the use of the Γ-point only or from the HSE06 functional, we repeated the calculations here with the Γ-point and either the PBE or the HSE06 functional. The results are still averaged over ten configurations. The value of σ obtained with the Γ-point only in BZ integration (Figure 4, left panel) is about 13% lower than the value obtained at convergence in the BZ integration (see Table 1 for PBE results). We therefore repeated the same calculations at the Γpoint with the HSE06 functional on snapshots from the PBE trajectories (see Figure 4). The differences due to the change of the functional is much larger than the error due to use of the Γ-point only. We expect the HSE06 functional to be more reliable than the PBE or rVV10 functional for the electronic properties and the transport coefficients because it better reproduces the bandgap in the semiconducting phase and thus also the pseudogap at the Fermi level which controls the electronic properties in the liquid phase. Unfortunately, the HSE06 functional is much more demanding in terms of the computational load than the PBE functional, which prevents us to perform a proper BZ integration with the HSE06 functional. Therefore, we report in Table 2 the   values of the transport coefficients computed at the Γ-point only with the HSE06 functional from either the PBE or rVV10 trajectories. We remark that the KS states at the Fermi level at 900 K are delocalized also with the HSE06 functional, as shown from the IPR in Figure 3 cobelli et al. [20] If we assume that the restriction of the BZ integration to the Γ-point only would lead to the same error in the PBE and HSE06 calculations, we might attempt to estimate the HSE06 value of σ at convergence in the BZ by scaling the Γ-point result by 13%, as it occurs for PBE. This leads to a value of about 1950 (Ωcm) À1 at 1000 K, which is closer, albeit still lower, to the experimental value of 2777 (Ωcm) À1 at 983 K. [15] This discrepancy might be also partially ascribed to the Ge-poor composition (Ge 1.6 Sb 2 Te 5 ) used in the study by Endo et al. [15] As a reference, the conductivity of Sb 2 Te 3 is 2293 (Ωcm) À1 at 992 K. [15] Turning now to the electronic contribution to the thermal conductivity, this is reasonably well reproduced by the application of the WF law to the calculated electrical conductivity, as shown in Tables 1,2. An exception is the rVV10 value of κ at 900 K, which shows the largest deviation from the WF result because of the closeness of the M-SC transition. [20] For the same reason, the Seebeck coefficient at 900 K with the rVV10 functional is somehow larger than the values at all other temperatures and than the PBE value at the same temperature. But for this point, the Seebeck coefficient shows a weak dependence on temperature and also a weak dependence on the choice of the functional always falling in the range 10-16 μV K À1 . The calculated Seebeck coefficient is very close to the value of about 20 μV K À1 used for liquid GST in electrothermal modeling in the study by Faraclas et al. [13] For the sake of completeness we also report the frequency dependent electrical conductivity σðωÞ for the HSE06 and PBE calculations over the PBE trajectories at 1000 K (see Figure 5). The electronic DOS at the three temperatures with the PBE and rVV10 functionals computed with the same 4 Â 4 Â 4 k-point mesh used for the calculation of the transport coefficients are shown instead in Figure S10 in Supporting Information, while the DOS at the Γ-point only calculated with PBE and HSE06 functionals on PBE trajectories at 1000 K is shown in Figure S11 in Supporting Information.

Conclusions
In summary, we provided DFT estimates of the transport and electrothermal coefficients of liquid GST that fall in the range of values used in the current electrothermal modeling of PCM devices reported in literature. In spite of some uncertainties, due to the dependence of the calculated electrothermal coefficients on the choice of the exchange and correlation functionals, the DFT results further support the reliability of the choice of the parameters for electrothermal modeling used so far in the literature and not directly available from experiments.   www.advancedsciencenews.com www.pss-rapid.com