arXiv:cond-mat/0408005v1 [cond-mat.mtrl-sci] 30 Jul 2004
Ab Initio Studies of Liquid and Amorphous Ge
Jeng-Da Chai∗
Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742
D. Stroud†
Department of Physics, The Ohio State University, Columbus, Ohio 43210 (Dated: February 2, 2008)
Abstract
We review our previous work on the dynamic structure factor S(k, ω) of liquid Ge ("-Ge) at temperature T = 1250 K, and of amorphous Ge (a-Ge) at T = 300 K, using ab initio molecular dynamics [Phys. Rev. B67, 104205 (2003)]. The electronic energy is computed using density- functional theory, primarily in the generalized gradient approximation, together with a plane wave representation of the wave functions and ultra-soft pseudopotentials. We use a 64-atom cell with periodic boundary conditions, and calculate averages over runs of up to about 16 ps. The calculated liquid S(k, ω) agrees qualitatively with that obtained by Hosokawa et al, using inelastic X-ray scattering. In a-Ge, we find that the calculated S(k, ω) is in qualitative agreement with that obtained experimentally by Maley et al. Our results suggest that the ab initio approach is sufficient to allow approximate calculations of S(k, ω) in both liquid and amorphous materials.
∗ Electronic address: jdchai@wam.umd.edu
† Electronic address: stroud@mps.ohio-state.edu
I. INTRODUCTION
Ge is a well-known semiconductor in its solid phase, but becomes metallic in its liquid phase. Liquid Ge (!-Ge) has, near its melting point, an electrical conductivity characteristic of a reasonably good metal (∼ 1.6×10−4Ω−1cm−1[1]), but it retains some residual structural features of the solid semiconductor. For example, the static structure factor S(k), has a shoulder on the high-k side of its first (principal) peak, which is believed to be due to a residual tetrahedral short-range order. This shoulder is absent in more conventional liquid metals such as Na or Al, which have a more close-packed structure in the liquid state and a shoulderless first peak in the structure factor. Similarly, the bond-angle distribution function just above melting is believed to have peaks at two angles, one near 60o and characteristic of close packing, and one near 108o, indicative of tetrahedral short range order. This latter peak rapidly disappears with increasing temperature in the liquid state.
These striking properties of !-Ge have been studied theoretically by several groups. Their methods fall into two broad classes: empirical and first-principles. A typical empirical calculation is that of Yu et al[2], who calculate the structural properties of !-Ge assuming that the interatomic potentials in !-Ge are a sum of two-body and three-body potentials of the form proposed by Stillinger and Weber[3]. These authors find, in agreement with experiment, that there is a high-k shoulder on the first peak of S(k) just above melting, which fades away with increasing temperature. However, since in this model all the potential energy is described by a sum of two-body and three-body interactions, the interatomic forces are probably stronger, and the ionic diffusion coefficient is correspondingly smaller, than their actual values.
In the second approach, the electronic degrees of freedom are taken explicitly into ac- count. If the electron-ion interaction is sufficiently weak, it can be treated by linear response theory[4]. In linear response, the total energy in a given ionic configuration is a term which is independent of the ionic arrangement, plus a sum of two-body ion-ion effective inter- actions. These interactions typically do not give the bond-angle-dependent forces which are present in the experiments, unless the calculations are carried to third order in the electron-ion pseudopotential[4], or unless electronic fluctuation forces are included[5] Such interactions are, however, included in the so-called ab initio approach, in which the forces on the ions are calculated from first principles, using the Hellman-Feynman theorem together
with density-functional theory[6, 7, 8] to treat the energy of the inhomogeneous electron gas.
This approach not only correctly gives the bond-angle-dependent ion-ion interactions, but also, when combined with standard molecular dynamics techniques, provides a good account of the electronic properties and such dynamical ionic properties as the ionic self-diffusion coefficients.
This combined approach, usually known as ab initio molecular dynamics, was pioneered by Car and Parrinello[9], and, in somewhat different form, has been applied to a wide range of liquid metals and alloys, including !-Ge[10, 11, 12], !-GaxGe1−x[13], stoichiometric III- V materials such as !-GaAs, !-GaP, and !-InP[14, 15], nonstoichiometric !-GaxAs1−x[16],
!-CdTe[17], and !-ZnTe[18], among other materials which are semiconducting in their solid phases. It has been employed to calculate a wide range of properties of these materials, in- cluding the static structure factor, bond-angle distribution function, single-particle electronic density of states, d. c. and a. c. electrical conductivity, and ionic self-diffusion coefficient.
The calculations generally agree quite well with available experiment.
A similar ab initio approach has also been applied extensively to a variety of amorphous semiconductors, usually obtained by quenching an equilibrated liquid state from the melt.
For example, Car, Parrinello and their collaborators have used their own ab initio approach (based on treating the Fourier components of the electronic wave functions as fictitious classical variables) to obtain many structural and electronic properties of amorphous Si[19, 20]. A similar approach has been used by Lee and Chang[21]. Kresse and Hafner[10]
obtained both S(k) and g(r), as well as many electronic properties, of a-Ge, using an ab initioapproach similar to that used here, in which the forces are obtained directly from the Hellmann-Feynman theorem and no use is made of fictitious dynamical variables for the electrons, as in the Car-Parrinello approach. A similar calculation for a-Si has been carried out by Cooper et al[22], also making use of a plane wave basis and treating the electron density functional in the generalized gradient approximation (GGA)[23]. More recently, a number of calculations for a-Si and other amorphous semiconductors have been carried out by Sankey et al[24], and by Drabold and collaborators [25]. These calculations use ab initio molecular dynamics and electronic density functional theory, but in a localized basis. A recent study, in which S(k) and g(r) were computed for several ab initio structural models of a-Si, has been carried out by Alvarez et al[26].
Finally, we mention a third approach, intermediate between empirical and ab initio molec-
ular dynamics, and generally known as tight-binding molecular dynamics. In this approach, the electronic part of the total energy is described using a general tight-binding Hamiltonian for the band electrons. The hopping matrix elements depend on separation between the ions, and additional terms are included to account for the various Coulomb energies in a consis- tent way. The parameters can be fitted to ab initio calculations, and forces on the ions can be derived from the separation-dependence of the hopping matrix elements. This approach has been used, e. g., to treat !-Si[27], a-Si[28], and liquid compound semiconductors such as !-GaAs and !-GaSb[29]. Results are in quite good agreement with experiment.
In this paper, we review the application of ab initio molecular dynamics to another dynamical property of the ions: the dynamical structure factor, denoted S(k, ω). While no fundamentally new theory is required to calculate S(k, ω), this quantity provides additional information about the time-dependent ionic response beyond what can be extracted from other quantities [30]. The work we review is the first to calculate S(k, ω) using ab initio molecular dynamics. Here, we will describe calculations of S(k, ω) for !-Ge, where some recent experiments[31] provide data for comparison, and also for amorphous Ge (a-Ge). In the latter case, using a series of approximations described below, we are able to infer the vibrationaldensity of states of as-quenched a-Ge near temperature T = 300K in reasonable agreement with experiment. The calculated S(k, ω) for the liquid also agrees quite well with experiment, especially considering the computational uncertainties inherent in an ab initio simulation with its necessarily small number of atoms and limited time intervals.
The remainder of this paper is organized as follows. A brief review of the theory and calculational method is given in Section II and III respectively. The results are reviewed in Section IV, followed by a discussion in Section V.
II. FINITE-TEMPERATURE DENSITY-FUNCTIONAL THEORY
The original Hohenberg-Kohn theorem of density-functional theory (DFT) [6] was gen- eralized to finite temperatures Tel by Mermin [8]. In this section, we give a brief review of the finite-temperature density-functional theory, in a form where it can be used in ab initio molecular dynamics calculations. We use atomic units here and throughout this article.
In the Born-Oppenheimer approximation, the total potential energy U of a system of N classical ions, enclosed in a volume V, and interacting with Nel=NZ valence electrons
through an electron-ion pseudopotential vei(r), is a functional of the valence electron density ρ(r), and also a function of the ion positions Rl (l=1,2,...,N). This potential energy can be expressed as the sum of a direct Coulombic interaction between the ions and the free energy Fel of the electronic system subject to the external potential created by the ions, which we denote Vext(r, {Rl}) =!Nl=1vei(|r − Rl|). This sum may be written
U[ρ(r), {Rl}] = Fel[ρ(r)] +"
l<l!
W (Rl− Rl!), (1)
where W (R − R#) = Z2/|R − R#| is the Coulomb interaction between ions of charge Z located at R and R#.
The free energy functional of the electronic system Fel[ρ(r)] can be partitioned into various terms as follows:[32, 33]
Fel[ρ(r)] = As[ρ(r)] +
#
ρ(r)Vext(r)dr + 1 2
# ρ(r)ρ(r#)
|r − r#| drdr#+ Fxc[ρ(r)]. (2) Here the different terms are As[ρ(r)], the Helmholz free energy functional for a reference system of noninteracting electrons; the interaction energy between the valence electrons and the external potential provided by the instantaneous ionic configuration; the classical electron-electron potential energy; and the exchange-correlation correction Fxc[ρ(r)] to the free energy. As[ρ(r)] has contributions from both the kinetic energy Ts[ρ(r)] and the entropy Ss[ρ(r)] of the reference system, and may be written
As[ρ(r)] = Ts[ρ(r)] − TelSs[ρ(r)], (3) where Tel is the temperature of the electronic system.
For a given set of ionic positions, the equilibrium electronic number density ρ(r) is ob- tained from the following variational principle:
δ
δρ(r)(U[ρ(r), {Rl}] − µ
#
ρ(r)dr) = 0. (4)
Here µ is the electron chemical potential, which is chosen to give the desired number of valence electrons Nel. The Euler equation associated with eq. (4) is
µ = VAs(r; [ρ]) + Vef f(r; [ρ]). (5) Here
VAs(r; [ρ]) = δAs[ρ]
δρ(r) (6)
is a functional derivative of As[ρ], and is the ”Helmholz free energy potential” (or the kinetic potential for Tel = 0) [34]. Vef f(r; [ρ]) is an effective one-body potential, given by
Vef f(r; [ρ]) = Vext(r) +
# ρ(r#)
|r − r#|dr#+ δFxc[ρ]
δρ(r) . (7)
It is usual to make several approximations for the functionals As[ρ] and Fxc[ρ], and the function Vext(r). First, Fxc[ρ] is only a minor contribution to Fel[ρ]. As a result, this func- tional is often calculated using the local density approximation (LDA) and or the generalized gradient approximation (GGA); both are widely used in many applications [32].
Furthermore, Fxc differs negligibly from its zero-temperature value, the exchange- correlation energy Exc, provided Tel < 0.15Tf, where Tf is the Fermi temperature [35].
Most applications involve temperatures satisfying this inequality. The external field Vext(r) contains the electron-ion pseudopotential vei(r), on which other approximations are made.
For example, it is common to approximate vei as local and energy-independent.
Once an accurately approximated As[ρ] is available, one can, in principle, obtain VAs(r; [ρ]) by functional differentiation, and from this functional, solve eq. (5) to obtain the equilibrium electronic density. However, this approach has proved to be very challenging in practice.
Typically, the kinetic energy Ts[ρ] is much larger than the entropy contribution TelSs[ρ]
in the range of temperature interested. Therefore, As[ρ] is generally well approximated by Ts[ρ]. However, simple local density approximations for Ts[ρ] based on the Thomas-Fermi (TF) theory are not very accurate, and make several qualitatively incorrect predictions [36, 37]. More recent developments based on so-called “orbital-free”(OF) approximations for the kinetic energy functional Ts[ρ] [38, 39, 40, 41], and the kinetic potential VTs(r; [ρ]) [42, 43] are still insufficient to make quantitative predictions. On the other hand, although these OF methods are currently unable to make accurate predictions, they are computation- ally appealing, because the computational time required to use them is independent of the number of valence electrons. Therefore, further refinement of such OF methods in DFT may lead to valuable new approaches for investigating the properties of large systems, especially metallic systems.
Instead of using approximate forms for As[ρ] or VAs(r; [ρ]), Kohn and Sham suggested a different way to compute the numerical value of As exactly, by introducing a set of Nel orbitals satisfying the coupled KS equations that describe the model system [7]. Using these
orbitals one can determine both As and the valence electron density ρ(r). The orbitals ψi satisfy the Schr¨odinger-like equation
$
−1
2∇2+ Vef f(r)
%
ψi = &iψi, (8)
and are chosen to be orthonormal. The set of equations is made self-consistent by requiring that
ρ(r) ="
i
fi|ψi(r).|2 (9)
Here fi is the Fermi occupation number, given by fi = f (&i− µ) = (exp(&i− µ
kBTel
) + 1)−1. (10)
The kinetic energy Ts is given by Ts[{ψi(r)}, {fi}] ="
i
fi
#
ψi(−1
2∇2)ψidr ="
i
&ifi−
#
ρ(r)Vef f(r)dr, (11)
and the electronic entropy Ss by
Ss[{fi}] = −2kB"
i
[filnfi+ (1 − fi)ln(1 − fi)] (12)
Eqs. (7,8,9,10) have to be solved self-consistently. The resulting valence electron density can be used to computed the electronic free energy Fel[ρ]; the electronic kinetic energy Ts and entropy Ss are given by eq(11) and eq(12) [32].
The use of finite-temperature DFT is important, in order to eliminate the instability of the KS equations at Tel = 0. In practice, Tel is usually regarded as a fictitious tempera- ture (chosen for maintaining the stability), and does not necessarily have to equal the ionic temperature T. Furthermore, the broadening function in eq. (10) need not be the Fermi func- tion, but may be chosen to be some other convenient function. Several possible broadening functions have been proposed [10, 44].
The forces on the individual ions can now be obtained using the Hellmann-Feynman theorem. According to this theorem, the force FI acting on the Ith ion is found by differ- entiating the potential energy of the system with respect to the ionic coordinate RI[45, 46].
Specifically,
FI = − ∂U
∂RI = −"
i
fi < ψi|∂Fel[ρ(r)]
∂RI |ψi > − ∂
∂RI
"
l<l!
W (Rl− Rl!) (13)
Given the force, one simply applies Newton’s second law to get the equations of motion of the N classical ions:
MIR¨I = FI (14)
In practice, the ions are treated classically, even though the origin of the force is very much quantum-mechanical.
In solving these equations of motion, it is important to maintain the ionic temperature T at a fixed value. This is typically accomplished by rescaling the total kinetic energy of the ions KI=!Nl=1 P
2 l
2Ml to satisfy
T = 2KI
3kB(N − 1) (15)
every few time steps. The factor of N-1, not N, appears here so as to conserve the total momentum of the system of ions.
In summary, for each instantaneous configuration of ions, one computes Fel[ρ(r)] and hence U[ρ(r), {Rl}] by solving the finite-temperature Kohn-Sham equations (7),(8),(9), and (10). The Hellmann-Feynman force FI is then computed from eq. (13). By numerically solving Newton’s equations of motion for the ions and rescaling the kinetic energy of the ions at every time step, one moves the ions to a new configuration. The entire process is repeated for a sufficiently large number of time steps to compute the desired property with good statistical accuracy.
III. METHOD
The method we have used to calculate the dynamic structure factor[30] is similar to that described in several previous papers for other properties of liquid semiconductors[12, 13, 16], but uses the Vienna Ab Initio Simulation Package (VASP), whose workings have been extensively described in the literature[47].” Briefly, the calculation involves two parts.
First, for a given ionic configuration, the total electronic energy is calculated, using an approximate form of the Kohn-Sham free energy density functional, and the force on each ion is also calculated, using the Hellmann-Feynman theorem. Second, Newton’s equations of motion are integrated numerically for the ions, using a suitable time step. The process is repeated for as many time steps as are needed to calculate the desired quantity. To hold the temperature constant, we use the canonical ensemble with the velocity rescaled at each time step. Further details of this approach are given in Ref. [12].
The VASP code uses the ultrasoft Vanderbilt pseudopotentials[48], a plane wave basis for the wave functions, with the original Monkhorst-Pack (3 × 3 × 3) k-space meshes[49] and a total of 21,952 plane waves, corresponding to an energy cut-off of 104.4eV. We have used a finite-temperature version of the Kohn-Sham theory[8] in which the electron gas Helmholtz free energy is calculated on each time step. This version also broadens the one-electron energy levels to make the k-space sums converge more rapidly. The method of Methfessel- Paxton of order 1 is used for the broadening function [44]., with the fictitious temperature for the electronic subsystem kBTel=0.2 eV. Most of our calculations were done using the generalized gradient approximation (GGA)[23] for the exchange-correlation energy (we use the particular form of the GGA developed by Perdew and Wang[23]), but some are also carried out using the local-density approximation (LDA).
In our iteration of Newton’s Laws in liquid Ge (!-Ge), we typically started from the diamond structure (at the experimental liquid state density for the temperature of interest), then iterated for 901 time steps, each of 10 fs, using the LDA. To obtain S(k) within the GGA, we started from the LDA configuration after 601 time steps, then iterated using the GGA for an additional 1641 10-fs time steps, or 16.41 ps. We calculate the GGA S(k) by averaging over an interval of 13.41 ps within this 16.41 time interval, starting a time t2 after the start of the GGA simulation. We average over all t2’s from 1.0 to 3.0 ps.
For comparison, we have also calculated S(k) within the LDA. This S(k) is obtained by averaging over 601 time steps of the 901 time-step LDA simulation. This 601-step interval is chosen to start a time t1 after the start of this simulation; the calculated LDA S(k) is also averaged over all t1’s from 1.0 ps to 3.0 ps.
To calculate quantities for amorphous Ge (a-Ge), we start with Ge in the diamond struc- ture at T = 1600K but at the calculated liquid density for that temperature, as given in Ref. [12]. Next, we quench this sample to 300 K, cooling at a uniform rate over, so as to reach 300 K in about 3.25 ps (in 10-fs time steps). Finally, starting from T = 300 K, we iterate for a further 897 time steps, each of 10 fs, or 8.97 ps, using the LDA. The LDA S(k) is then obtained by averaging over 5.97 of those 8.97 ps, starting a time t1 after the system is reached 300K; we also average this S(k) over all t1’s from 1.0 to 3.0 ps. To obtain S(k) within the GGA, we start the GGA after 5.7 ps of the LDA simulation, and then iterate using the GGA for an additional 18.11 ps in 10-fs time steps. The GGA S(k) is obtained by averaging over a 15.11 ps time interval of this 18.11 ps run, starting a time t2 after the
start of the GGA simulation; we also average over all values of t2 from 1.0 to 3.0 ps.
The reader may be concerned that the 3.25 ps quench time is very short, very unrepre- sentative of a realistic quench, and very likely to produce numerous defects in the quenched structure, which are not typical of the a-Ge studied in experiments. In defense, we note that some of these defects may be annealed out in the subsequent relaxation at low temperatures, which is carried out before the averages are taken. In addition, the static structure factor we obtain agrees well with experiment, and the dynamic structure factor is also consistent with experiment, as discussed below. Thus, the quench procedure appears to produce a material rather similar to some of those studied experimentally. Finally, we note that exper- iments on a-Ge themselves show some variation, depending on the exact method of sample preparation.
We have used the procedure outlined above to calculate various properties of !-Ge and a- Ge. Most of these calculated properties have been described in previous papers, using slightly different methods, and therefore will be discussed here only very briefly[50]. However, our results for the dynamic structure factor S(k, ω) as a function of wave vector k and frequency ω are new, and will be described in detail. We also present our calculated static structure factor S(k), which is needed in order to understand the dynamical results.
S(k) is defined by the relation S(k) = 1
N%"
i,j
exp[ik · (ri(t) − rj(t))]&t0 − Nδk,0, (16)
where riis the position of the ithion at time t, N is the number of ions in the sample, and the triangular brackets denote an average over the sampling time. In all our calculations, we have used a cubic cell with N = 64 and periodic boundary conditions in all three directions. This choice of particle number and cell shape is compatible with any possible diamond-structure Ge within the computational cell.
S(k, ω) is defined by the relation (for k '= 0, ω '= 0) S(k, ω) = 1
2πN
# ∞
−∞exp(iωt)%ρ(k, t)ρ(−k, 0)&dt, (17) where the Fourier component ρ(k, t) of the number density is defined by
ρ(k, t) =
N
"
i=1
exp(−ik · ri(t)). (18)
In our calculations, the average %...& is computed as
%ρ(k, t)ρ(−k, 0)& = 1
∆t1
# ∆t1
0 ρ(k, t1 + t)ρ(−k, t1)dt1 (19) over a suitable range of initial times t1. Because of the expected isotropy of the liquid or amorphous phase, which should hold in the limit of large N, S(k, ω) should be a function only of the magnitude k rather than the vector k, as should the structure factor S(k).
Our calculations are carried out over relatively short times. To reduce statistical errors, we therefore first calculate
S(k, ω, t1, t2) = 1 πN
# t2
0 dtρ(k, t1+ t)ρ(−k, t1) exp(iωt). (20) For large enough t2, S(k, ω, t1, t2) should become independent of t2 but will still retain some dependence on t1. Therefore, in the liquid, we obtain our calculated dynamic structure factor, Scalc(k, ω), by averaging over a suitable range of t1 from 0 to ∆t1:
Scalc(k, ω) = 1
∆t1
# ∆t1
0 dt1S(k, ω, t1, t2). (21) We choose our initial time in the t1 integral to be 1 ps after the start of the GGA calculation, and (in the liquid) ∆t1 = 6 ps. We choose the final MD time t2 = 16.41ps. For a-Ge, we use the same procedure but t2 = 18.11ps in our simulations. For our finite simulational sample, the calculated S(k) and S(k, ω) will, in fact, depend on the direction as well as the magnitude of k. To suppress this finite-size effect, we average the calculated S(k) and S(k, ω) over all k values of the same length. This averaging considerably reduces statistical error in both S(k, ω) and S(k).
Finally, we have also incorporated the experimental resolution functions into our plotted values of S(k, ω). Specifically, we generally plot Sav(k, ω)/S(k), where Sav(k, ω) is obtained from the (already orientationally averaged) S(k, ω) using the formula
Sav(k, ω) =
# ∞
−∞R(ω − ω#)S(k, ω#)dω#, (22) where the resolution function R(ω) (normalized so that &−∞∞ R(ω)dω = 1) is
R(ω) = 1
√πω0 exp(−ω2/ω02). (23)
In an isotropic liquid, we must have S(k, −ω) = S(k, ω), since our ions are assumed to move classically under the calculated first-principles force. Our orientational averaging
procedure guarantees that this will be satisfied identically in our calculations, since for every k, we always include −k in the same average. We will nonetheless show results for negative ω for clarity, but they do not provide any additional information.
IV. RESULTS
A. S(k) for "-Ge and a-Ge
In Fig. 1, we show the calculated S(k) for !-Ge at T = 1250 K, as obtained using the procedure described in Section II. The two calculated curves are obtained using the GGA and the LDA for the electronic energy-density functional; they lead to nearly identical results. The calculated S(k) shows the well-known characteristics already found in previous simulations[10, 12]. Most notably, there is a shoulder on the high-k side of the principal peak, which is believed to arise from residual short-range tetrahedral order persisting into the liquid phase just above melting. We also show the experimental results of Waseda et al[51]; agreement between simulation and experiment is good, and in particular the shoulder seen in experiment is also present in both calculated curves (as observed also in previous simulations).
We have also calculated S(k) for a model of amorphous Ge (a-Ge) at 300K. We prepared our sample of a-Ge as described in the previous section. As for !-Ge, we average the calculated S(k) over different k vectors of the same length, as for !-Ge. In Fig. 2, we show the calculated S(k) for a-Ge at T = 300, again using both the GGA and the LDA. The sample is prepared and the averages obtained as described in Section II. In contrast to !-Ge, but consistent with previous simulations[10, 26], the principal peak in S(k) is strikingly split.
The calculations are in excellent agreement with experiments carried out on as-quenched a- Ge at T = 300 K[52]; in particular, the split principal peak seen in experiment is accurately reproduced by the simulations.
We have also calculated a number of other quantities for both !-Ge and a-Ge, including pair distribution function g(r), and the electronic density of states n(E). For !-Ge, we calculated n(E) using the Monkhorst-Pack mesh with gamma point shifting (one of the meshes recommended in the VASP package). The resulting n(E) is generally similar to that found in previous calculations[10, 12], provided that an average is taken over at least 5-10
liquid state configurations. Our n(E) for a-Ge [calculated using a shorter averaging time than that used below for S(k, ω)] is also similar to that found previously[10]. Our calculated g(r)’s for both !-Ge and a-Ge, as given by the VASP program, are similar to those found in Refs. [10] and [12]. The calculated number of nearest neighbors in the first shell is 4.18 for a-Ge measured to the first minimum after the principal peak in g(r). For !-Ge, if we count as “nearest neighbors” all those atoms within 3.4˚Aof the central atom (the larger of the cutoffs used in Ref. [10]) we find approximately 7.2 nearest neighbors, quite close to the value of 6.9 obtained in Ref. [10] for that cutoff. Finally, we have recalculated the self- diffusion coefficient D(T ) for !-Ge at T = 1250 K, from the time derivative of the calculated mean-square ionic displacement; we obtain a result very close to that of Ref. [12].
B. S(k, ω) for "-Ge and a-Ge
1. "-Ge
Fig. 3 shows the calculated ratio S(k, ω)/S(k) for !-Ge at T = 1250 K, as obtained using the averaging procedure described in Section II. We include a resolution function [eqs. (22) and (23)] of width ¯hω = 2.5 meV, the same as the quoted experimental width[31]. In Fig. 4, we show the same ratio, but without the resolution function (i. e., with ω0 = 0). Obviously, there is much more statistical noise in this latter case, though the overall features can still be distinguished.
To interpret these results, we first compare the calculated S(k, ω) in !-Ge with hydrody- namic predictions, which should be appropriate at small k and ω. The This prediction takes the form (see, for example, Ref. [53]):
2πS(k, ω)
S(k) = γ − 1 γ
' 2DTk2 ω2+ (DTk2)2
(
+ + 1
γ
' Γk2
(ω + csk)2+ (Γk2)2 + Γk2
(ω − csk)2+ (Γk2)2
(
. (24)
Here γ = cP/cV is the ratio of specific heats and constant pressure and constant volume, DT is the thermal diffusivity, cs is the adiabatic sound velocity, and Γ is the sound attenuation constant. DT and Γ can in turn be expressed in terms of other quantities. For example, DT = κT/(ρcP), where κT is the thermal conductivity and ρ is the atomic number density.
Similarly, Γ = 12[a(γ − 1)/γ + b], where a = κT/(ρcV) and b is the kinematic longitudinal
viscosity (see, for example, Ref. [53], pp. 264-66).
Eq. (24) indicates that S(k, ω) in the hydrodynamic regime should have two propagating peaks centered at ω = ±csk, and a diffusive peak centered at ω = 0 and of width determined by DT. The calculated S(k, ω)/S(k) for the three smallest values of k in Fig. 3, does show the propagating peaks. We estimate peak values of ¯hω ∼ 10 meV for k = 5.60 nm−1, ¯hω ∼ 11 meV for k = 7.92 nm−1, and (somewhat less clearly) ¯hω ∼ 13 meV for k = 9.70 nm−1. The value of cs estimated from the lowest k value is cs ∼ 2.7 × 105 cm/sec. (The largest of these three k values may already be outside the hydrodynamic, linear-dispersion regime.)
These predictions agree reasonably well with the measured S(k, ω) obtained by Hosokawa et al[31], using inelastic X-ray scattering. For example, the measured sound-wave peaks for k = ±6 nm−1 occur near 10 meV, while those k = ±12 nm−1 occur at ¯hω = 17.2 meV, Furthermore, the integrated relative strength of our calculated sound-wave peaks, compared to that of the central diffusion peak, decreases between k = 7.92 nm−1 and 12.5 nm−1, consistent with both eq. (24) and the change in experimental behavior[31] between k = 6 nm−1 and 12 nm−1.
Because S(k, ω) in Fig. 3 already includes a significant Gaussian smoothing function, a quantitatively accurate half-width for the central peak, and hence a reliable predicted value for DT, cannot be extracted. A rough estimate can be made as follows. For the smallest k value of 5.6 nm−1, the full width of the central peak at half-maximum is around 7.5 meV. If the only broadening were due to this Gaussian smoothing, the full width would be around 2¯hω0 = 5.5 meV. Thus, a rough estimate of the intrinsic full width is ≈ √
7.52− 5.52 = 5 meV ≈ 2¯hDTk2. This estimate seems reasonable from the raw data for S(k, ω) shown in Fig. 4. Using this estimate, one obtains DT ≈ 1.3 × 10−3 cm2/sec.
The hydrodynamic expression for S(k, ω)/S(k) was originally obtained without consid- eration of the electronic degrees of freedom. Since !-Ge is a reasonably good metal, one might ask if the various coefficients appearing eq. (24) should be the full coefficients, or just the ionic contribution to those coefficients. For example, should the value of DT which determines the central peak width be obtained from the full cP, cV, and κT, or from only the ionic contributions to these quantities? For !-Ge, the question is most relevant for κT, since the dominant contribution to cP and cV should be the ionic parts, even in a liquid metal[4].
However, the principal contribution to κT is expected to be the electronic contribution.
We have made an order-of-magnitude estimate of DT using the experimental liquid num-
ber density and the value CP = (5/3)kB per ion, and obtaining the electronic contribution to κT from the Wiedemann-Franz law[54] together with previously calculated estimates of the electronic contribution[12]. This procedure yields DT ∼ 0.1 cm2/sec, about two orders of magnitude greater than that extracted from Fig. 3, and well outside the possible errors in that estimate. We conclude that the DT which should be used in eq. (24) for !-Ge (and by inference other liquid metals) is the ionic contribution only.
In support of this interpretation, we consider what one expects for S(k, ω) in a simple metal such as Na. In such a metal, ionic motions are quite accurately determined by effective pairwise screened ion-ion interactions[4]. Since the ionic motion is determined by such an interaction, the S(k, ω) resulting from that motion should not involve the contribution of the electron gas to the thermal conductivity. Although !-Ge is not a simple metal, it seems plausible that its S(k, ω) should be governed by similar effects, at least in the hydrodynamic regime. This plausibility argument is supported by our numerical results.
For k beyond around 12 nm−1, the hydrodynamic model should start to break down, since the dimensionless parameter ωτ (where τ is the Maxwell viscoelastic relaxation time) becomes comparable to unity. At these larger k’s, both our calculated and the measured[31]
curves of S(k, ω)/S(k) continue to exhibit similarities. Most notable is the existence of a single, rather narrow peak for k near the principal peak of S(k), followed by a reduction in height and broadening of this central peak as k is further increased. This narrowing was first predicted by de Gennes[55]. In our calculations, it shows up in the plot for k = 20.9 nm−1, for which the half width of S(k, ω)/S(k) is quite narrow, while at k = 28.5 and 30.7 nm−1, the corresponding plots are somewhat broader and lower. By comparison, the measured central peak in S(k, ω)/S(k) is narrow at k = 20 nm−1 and especially at k = 24 nm−1, while it is broader and lower at k = 28 nm−1[31].
The likely physics behind the de Gennes narrowing is straightforward. The half-width of S(k, ω) is inversely proportional to the lifetime of a density fluctuation of wave number k.
If that k coincides with the principal peak in the structure factor, a density fluctuation will be in phase with the natural wavelength of the liquid structure, and should decay slowly, in comparison to density fluctuations at other wavelengths. This is indeed the behavior observed both in our simulations and in experiment.
In further support of this picture, we may attempt to describe these fluctuations by a very oversimplified Langevin model. We suppose that the Fourier component ρ(k, t) [eq.
(18)] is governed by a Langevin equation
˙ρ(k, t) = −ξρ(k, t) + η(t). (25)
Here the dot is a time derivative, ξ is a constant, and η(t) is a random time-dependent “force”
which has ensemble average %η& = 0 and correlation function %η(t)η∗(t#)& = Aδ(t − t#). Eq.
(25) can be solved by standard methods (see, e. g., Ref. [53] for a related example), with the result (for sufficiently large t)
%ρ(k, t)ρ∗(k, t + τ )& = πA
ξ exp(−|τ|ξ). (26)
According to eq. (17), S(k, ω) is, to within a constant factor, the frequency Fourier transform of this expression, i. e.
S(k, ω) ∝
# ∞
−∞(πA/ξ) exp(iωτ ) exp(−|τ|ξ)dτ, (27) or, on carrying out the integral,
S(k, ω) ∝ πA
ω2+ ξ2. (28)
This is a Gaussian function centered at ω = 0 and of half-width ξ. On the other hand, the static structure factor
S(k) ∝ Limτ →0%ρ(k, t)ρ∗(k, t + τ )& = πA
ξ . (29)
Thus, if the constant A is independent of k, the half-width ξ of the function S(k, ω) at wave number k is inversely proportional to the static structure S(k). This prediction is consistent with the “de Gennes narrowing” seen in our simulations and in experiment[31].
To summarize, there is overall a striking similarity in the shapes of the experimental and calculated curves for S(k, ω)/S(k) both in the hydrodynamic regime and at larger values of k.
2. a-Ge
We have also calculated the dynamic structure for our sample of a-Ge at T = 300K.
The results for the ratio S(k, ω)/S(k) are shown in Figs. 5 and 6 for a range of k values, and, over a broader range of ω, in Fig. 7. Once again, both S(k, ω) and S(k) are averaged
over different values of k of the same length, as described above. We have incorporated a resolution function of width ¯hω0 = 2 meV into S(k, ω). This width is a rough estimate for the experimental resolution function in the measurements of Maley et al[56]; we assume it to be smaller than the liquid case because the measured width of the central peak in S(k, ω)/S(k) for a-Ge is quite small.
Ideally, our calculated S(k, ω) should be compared to the measured one. However, the published measured quantity is not S(k, ω) but is, instead, based on a modified dynamical structure factor, denoted G(k, ω), and related to S(k, ω) by[56]
G(k, ω) =
)C k2
*'
¯hω n(ω, T ) + 1
(
S(k, ω). (30)
Here C is a k- and ω-independent constant, and n(ω, T ) = 1/[e¯hω/kBT − 1] is the phonon occupation number for phonons of energy E = ¯hω at temperature T . The quantity plotted by Maley et al[56] is an average of G(k, ω) over a range of k values from 40 to 70 nm−1. These workers assume that this average is proportional to the vibrational density of states nvib(ω). The measured nvib(ω) as obtained in this way[56] is shown in Fig. 8 for two different amorphous structures, corresponding to two different methods of preparation and having differing degrees of disorder.
In order to compare our calculated S(k, ω) to experiment, we use eq. (30) to infer G(k, ω), then average over a suitable range of k. However, in using eq. (30), we use the classical form of the occupation factor, n(¯hω) + 1 ≈ kBT /¯hω, This choice is justified because we have calculated S(k, ω) using classical equations of motion for the ions. We thus obtain for the calculated vibrational density of states
nvib(ω) ≈
'C¯h2 kBT
( 'ω2 k2
(
S(k, ω). (31)
In Fig. 8 we show two such calculated plots of nvib(ω), as obtained by averaging eq. (31) over two separate groups of k’s, as indicated in the caption[57]. For comparison, we also show nvib(ω) for a-Ge as calculated in Ref. [10] directly from the Fourier transform of the velocity-velocity autocorrelation function.
The calculated plots for nvib(ω) in Fig. 8 have some distinct structure, which arises from some corresponding high frequency structure in S(k, ω). The plot of nvib(ω) for the group of smaller k’s has two distinct peaks, near 8 meV and 29 meV, separated by a broad dip with a minimum near 18 meV. The plot corresponding to the group of larger k’s has similar
structure and width, but the dip is less pronounced. The two experimental plots also have two peaks separated by a clear dip. The two maxima are found around 10 and 35 meV, while the principal dip occurs near 16 meV. In addition, the overall width of the two densities of states is quite similar.
The reasonable agreement between the calculated and measured nvib(ω) suggests that our ab initio calculation of S(k, ω) for a-Ge is reasonably accurate. The noticeable differences probably arise from several factors. First, there are several approximations involved in going from the calculated and measured S(k, ω)’s to the corresponding nvib(ω)’s, and these may be responsible for some of the discrepancies. Secondly, there may actually be differences between the particular amorphous structures studied in the experiments, and the quenched, then relaxed structure considered in the present calculations. (However, the similarities in the static structure factors suggest that these differences are not vast.) Finally, our calculations are carried out over relatively short times, using relatively few atoms; thus, finite-size and finite-time effects are likely to produce some additional errors. Considering all these factors, agreement between calculation and experiment is quite reasonable.
Previous ab initio calculations for a-Ge[10] have also obtained a vibrational density of states, but this is computed directly from the ionic velocity-velocity autocorrelation function rather than from the procedure described here. The calculations in Ref. [10] do not require computing S(k, ω). In the present work, by contrast, we start from our calculated S(k, ω), and we work backwards to get nvib(ω). In principle, our S(k, ω) includes all anharmonic effects on the vibrational spectrum of a-Ge, though in extracting nvib(ω) we assume that the lattice vibrates harmonically about the metastable atomic positions. In Fig. 8, we also show the results of Ref. [10] for nvib(ω) as obtained from this correlation function. They are quite similar to those obtained in the present work, but have a somewhat deeper minimum between the two principal peaks.
The quantity nvib(ω) could, of course, also be calculated directly from the force constant matrix, obtained by assuming that the quenched configuration is a local energy minimum and calculating the potential energy for small positional deviations from that minimum using ab initiomolecular dynamics. This procedure has been followed for a-GeSe2, for example, by Cappelletti et al[58]. These workers have then obtained S(q, ω) versus q from their nvib(ω) at selected values of ω, within a one-phonon approximation. However, as noted above, the present work produces the full S(k, ω) and thus has, in principle, more information than
nvib(ω).
V. DISCUSSION AND CONCLUSIONS
The results reviewed here show that ab initio molecular dynamics can be used to calcu- late the dynamic structure factor S(k, ω) for both liquid and amorphous semiconductors.
Although the accuracy of the calculated S(k, ω) is lower than that attained for static quan- tities, such as S(k), nonetheless it is sufficient for comparison to most experimental features.
This is true even though our calculations are limited to 64-atom samples and fewer than 20 ps of elapsed real time.
We have presented evidence that the calculated S(k, ω)/S(k) in !-Ge agrees qualitatively with measured by inelastic X-ray scattering[31], and that the one calculated for a-Ge leads to a vibrational density of states qualitatively similar to the quoted experimental one[56].
Since such calculations are thus shown to be feasible, the work reviewed here should spur further numerical studies, with longer runs on larger samples, to obtain even more detailed information. Furthermore, we can use these dynamical simulations to probe the underly- ing processes at the atomic scale which give rise to specific features in the measured and calculated S(k, ω).
VI. ACKNOWLEDGEMENTS
This work has been supported by NSF grants DMR01-04987 (D.S.), DMR04-12295 (D.S.) and CHE01-11104 (J.D.C.). Calculations were carried out using the Beowulf Cluster at the Ohio Supercomputer Center, with the help of a grant of time. Jeng-Da Chai acknowledges full support from the UMCP Graduate School Fellowship, the IPST Alexander Family Fel- lowship, and the CHPH Block Grant Supplemental Fellowship.
[1] V. M. Glazov, S. N. Chizhevskaya, and N. N. Glagoleva, Liquid Semiconductors (Plenum, New York, 1969).
[2] W. B. Yu, Z. Q. Wang, and D. Stroud, Phys. Rev. B54, 13946 (1996).
[3] F. H. Stillinger and Weber, Phys. Rev. B31, 5262 (1985).
[4] For a review, see, e. g., N. W. Ashcroft and D. Stroud, Solid State Physics 33, pp. 1ff (1978).
[5] N. W. Ashcroft, Il Nuovo Cimento 12D, 597 (1990).
[6] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
[7] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[8] N. D. Mermin, Phys. Rev. 137, A1441 (1965).
[9] R. Car and M. Parrinello, Phys. Rev. Lett. 35, 2471 (1985).
[10] G. Kresse and J. Hafner, Phys. Rev. B49, 14251 (1994).
[11] N. Takeuchi and I. L. Garz´on, Phys. Rev. B50, 8342 (1995).
[12] R. V. Kulkarni, W. G. Aulbur, and D. Stroud, Phys. Rev. B55, 6896 (1997).
[13] R. V. Kulkarni and D. Stroud, Phys. Rev. B57, 10476 (1998).
[14] Q. Zhang, G. Chiarotti, A. Selloni, R. Car, and M. Parrinello, Phys. Rev. B42, 5071.
[15] L. Lewis, A. De Vita, and R. Car, Phys. Rev. B57, 1594 (1998).
[16] R. V. Kulkarni and D. Stroud, Phys. Rev. B62, 4991 (2000).
[17] V. Godlevsky, J. Derby, and J. Chelikowsky, Phys. Rev. Lett. 81, 4959 (1998); V. Godlevsky, M. Jain, J. Derby, and J. Chelikowsky, Phys. Rev. B60, 8640 (1999);
[18] M. Jain, V. V. Godlevsky, J. J. Derby, and J. R. Chelikowsky, Phys. Rev. B65, 035212 (2001).
[19] R. Car and M. Parrinello, Phys. Rev. Lett. 63, 204 (1988)
[20] I Stich, R. Car, and M. Parrinello, Phys. Rev. Lett. 63, 2240 (1989); Phys. Rev. B44, 4261 (1991); Phys. Rev. B44, 11092 (1991).
[21] I. Lee and K. J. Chang, Phys. Rev. B50, 18083 (1994).
[22] N. C. Cooper, C. M. Goringe, and D. R. McKenzie, Comput. Mater. Sci. 17, 1 (2000) [23] See, e. g., J. P. Perdew, in Electronic Structure of Solids, edited by P. Ziesche and H. Eschrig
(Akademic Verlag, Berlin, 1991).
[24] O. Sankey and D. J. Niklewsky, Phys. Rev. B40, 3979 (1989)
[25] D. A. Drabold, P. A. Fedders, O. F. Sankey, and J. D. Dow, Phys. Rev. B42, 5135 (1990)
[26] F. Alvarez, C. C. D´iaz, A. A. Valladares, and R. M. Valladares, Phys. Rev. B65, 113108 (2002).
[27] C. Z. Wang, C. T. Chan, and K. M. Ho, Phys. Rev. B45, 12227 (1992); I. Kwon, R. Biswas, C. Z. Wang, K. M. Ho, and C. M. Soukoulis, Phys. Rev. B49, 7242 (1994).
[28] G. Servalli and L. Colombo, Europhys. Lett. 22, 107 (1993).
[29] C. Molteni, L. Colombo, and Miglio, J. Phys.: Cond. Matt. 6, 5255 (1994).
[30] J.-D. Chai, D. Stroud, J. Hafner, and G. Kresse, Phys. Rev. B67, 104205 (2003).
[31] S. Hosokawa, Y. Kawakita, W.-C. Pilgrim, and H. Sinn, Phys. Rev. B63, 134205 (2001).
[32] R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, (Oxford Uni- versity Press, 1989).
[33] R. M. Dreizler and E. K. U. Gross, Density Functional Theory: An Approach to the Quantum Many Body Problem, (Springer-Verlag, Berlin, 1990).
[34] R. A. King and N. C. Handy, Phys. Chem. Chem. Phys. 2, 5049 (2000).
[35] D. G. Kanhere, P. V. Panat, A. K. Rajagopal, and J. Callaway, Phys. Rev. A33, 490 (1986).
[36] L. H. Thomas, Proc. Cambridge Phil. Soc. 23, 542 (1927).
[37] E. Fermi, Z. Phys. 48, 73 (1928).
[38] See e.g., Y. A. Wang and E. A. Carter, in Theoretical Methods in Condensed Phase Chemistry, edited by S.D. Schwartz, Progress in Theoretical Chemistry and Physics, (Kluwer, Boston, 2000) p. 117, and references therein.
[39] L.-W. Wang and M. P. Teter, Phys. Rev. B45, 13196 (1992).
[40] Y. A. Wang, N. Govind, and E. A. Carter, Phys. Rev. B58, 13465 (1998).
[41] Y. A. Wang, N. Govind, and E. A. Carter, Phys. Rev. B60, 16350 (1999).
[42] J.-D. Chai and J. D. Weeks, J. Phys. Chem. B108, 6870 (2004).
[43] J.-D. Chai and J. D. Weeks, (unpublished).
[44] M. Methfessel and A. T. Paxton, Phys. Rev. B40, 3616 (1989).
[45] M. Weinert and J. W. Davenport, Phys Rev. B45, 13709 (1992).
[46] R. M. Wentzcovitch, J. L. Martins, and P. B. Allen, Phys. Rev. B45, 11372 (1992).
[47] G. Kresse and J. Hafner, Phys. Rev. B47, 558 (1993); G. Kresse, Thesis, Technische Univer- sit¨at Wien (1993); G. Kresse and J. Furthm¨uller, Comput. Mat. Sci. 6, pp. 15-50 (1996); G.
Kresse and J. Furthm¨uller, Phys. Rev. B54, 11169 (1996).
[48] D. Vanderbilt, Phys. Rev. B32, 8412 (1985).
[49] H. J. Monkhorst and J. D. Pack, Phys. Rev. B13, 5188 (1976).
[50] Our method is very similar to that used in Ref. [10] to treat "-Ge and a-Ge, but with the following differences: (i) we use both the GGA and the LDA, while Ref. [10] uses only the LDA;
(ii) we include slightly more bands (161 as compared to 138); and (iii) most of our calculations are carried out using a 10 fs time step, rather than the 3 fs step used in Ref. [10]. Our a-Ge structure may also differ slightly from that of Ref. [10], since we use a somewhat faster quench rate (∼ 3 × 1014 K/sec rather than 1.7 × 1014 K/sec). Both codes use the conjugate gradient method to find energy minima and the same Vanderbilt ultrasoft pseudopotentials.
[51] Y. Waseda, The Structure of Noncrystalline Materials: Liquids and Amorphous Solids (McGraw-Hill, 1980).
[52] G. Etherington, A. C. Wright, J. T. Wenzel, J. T. Dore, J. H. Clarke, and R. N. Sinclair, J.
Non-Cryst. Solids 48, 265 (1982).
[53] See, for example, J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, 2nd edition (Academic Press, London, 1986), pp. 222-228.
[54] See, for example, N. W. Ashcroft and N. D. Mermin, Solid State Physics (Harcourt-Brace, Fort Worth, 1976), ch. 2.
[55] P. G. de Gennes, Physics 25, 825 (1959); Physics 3, 37 (1967).
[56] N. Maley, J. S. Lannin, and D. L. Price, Phys. Rev. Lett. 56, 1720 (1986).
[57] Two recent groups [V. Mart´in-Mayor, M. M´ezard, G. Parisi, and P. Verrocchio, J. Chem.
Phys. 114, 8068 (2001), and T. S. Grigera, V. Mart´in-Mayor, G. Parisi, and P. Verrocchio, Phys. Rev. Lett. 87, 085502 (2001)] have derived approximate expressions connecting S(k, ω) to nvib(ω) in an amorphous solid, which closely resemble eq. (31). In fact, their expressions can be shown to become equivalent to eq. (31) at sufficiently large k.
[58] R. L. Cappelletti, M. Cobb, D. A. Drabold, and W. A. Kamitakahara, Phys. Rev. B52, 9133 (1995).
0 1 2 3 4 5 6 7 8 9 10 0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
k (10 nm−1)
S(k)
GGA LDA Expt.
FIG. 1: Static structure factor S(k) for "-Ge at T = 1250K, just above the experimental melting temperature. Full curve: present work, as calculated using the generalized gradient approximation (GGA; see text). Dashed curve: present work, but using the local density approximation (LDA;
see text). Open circles: measured S(k) near T = 1250 K, as given in Ref. [51].
0 1 2 3 4 5 6 7 8 9 10 0
0.5 1 1.5 2 2.5
k (10 nm−1)
S(k)
GGA LDA Expt.
FIG. 2: Full curve: Calculated S(k) for a-Ge at T = 300K, as obtained using the GGA. Structure is prepared as described in the text. Open circles: measured S(k) for a-Ge at T = 300K, as given in Ref. [52].
−30 −20 −10 0 10 20 30 0
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
ω (meV)
S(k,ω)/S(k)
5.60 7.92 9.70 11.2 12.5 20.9 25.7 28.5 30.7 nm−1
FIG. 3: Calculated ratio of dynamic structure factor S(k, ω) to static structure factor S(k) for "-Ge at T = 1250K for several values of k, plotted as a function of ω, calculated using ab initio molecular dynamics with a MD time step of 10 fs. For clarity, each curve has been vertically displaced by 0.05 units from the curve below. In each case, the plotted curve is obtained by averaging both the calculated S(k, ω) and the calculated S(k) over all values of k of the same length. We also incorporate a Gaussian resolution function of half-width ¯hω0 = 2.5 meV, as in eqs. (22) and (23).
This value of ω0 is chosen to equal the quoted experimental resolution[31].
−30 −20 −10 0 10 20 30 0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
ω (meV)
S(k,ω)/S(k)
5.60 7.92 9.70 11.2 12.5 20.9 25.7 28.5 30.7 nm−1
FIG. 4: Same as in Fig. 3, but without the resolution function.
−300 −20 −10 0 10 20 30 0.5
1 1.5 2
ω (meV)
S(k,ω)/S(k)
5.53 7.82 9.58 11.1 12.4 20.7 25.4 28.2 30.3 34.1 nm−1
FIG. 5: Calculated S(k, ω)/S(k) for a-Ge at T = 300 K at k ≤ 35 nn−1, plotted as a function of ω Again, each curve has been vertically displaced by appropriate amounts from the one below it, as evident from the Figure, and both S(k, ω) and S(k) have been plotted after an average over all k’s of the same length. We also incorporate a Gaussian resolution function of half-width ¯hω0 = 2 meV. This value is chosen to give the best results for nvib(ω) as measured by Ref. [56]. The time step here is 10 fs.
−300 −20 −10 0 10 20 30 0.2
0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
ω (meV)
S(k,ω)/S(k)
39.1 43.2 48.6 84.6 86.6 87.7 90.2 90.9 93.4 nm−1
FIG. 6: Same as Fig. 5, but at k ≥ 35 nm−1.
5 10 15 20 25 30 35 40 45 50 55 60 0
0.005 0.01 0.015 0.02 0.025 0.03 0.035
ω (meV)
S(k,ω)/S(k)
39.1 nm−1 43.2 48.6 84.6 86.6 87.7 90.2 90.9 93.4
FIG. 7: Calculated S(k, ω)/S(k) as in Figs. 5 and 6 but including higher frequencies ω. At ¯hω = 60 meV, the curves are arranged vertically in order of increasing frequency. Each curve is vertically displaced by 0.0005 units from the one below it.
0 5 10 15 20 25 30 35 40 45 50 0
5 10 15 20 25 30 35 40 45 50 55
ω (meV) n vib(ω) (10−3 meV−1 )
3k near 40 nm−1 6k near 90 nm−1 Expt. (disordered) Expt. (ordered) Ab Initio (Kresse)
FIG. 8: Full curve: calculated vibrational density of states nvib(ω), in units of 10−3 states/meV.
nvib(ω) is obtained from the resolution-broadened S(k, ω) and S(k) of the previous Figure, using the formula nvib(ω) = %Gcalc(k, ω)&, where Gcalc(k, ω) is given by eq. (16), and the averaging is carried out over the three magnitudes of k near 40 nm−1 for which we have computed S(k, ω).
Dashed curve: same as full curve, but calculated by averaging over the six magnitudes of k near 90 nm−1 for which we have computed S(k, ω). The open circles and open stars denote the measured nvib(ω), as reported in Ref. [56] for two forms of a-Ge. Finally, the open diamonds denote nvib(ω) as calculated in Ref. [10] from the ionic velocity-velocity autocorrelation function (dot-dashed curves).
In all plots except that of Ref. [10], nvib(ω) is normalized so that &0ωmaxn(ω)d(¯hω) = 1. ωmax is the frequency at which nvib(ω) → 0, and is estimated from this Figure by extrapolating the right hand parts of the solid and dashed curves linearly to zero.
