Dynamic structure factor of liquid and amorphous Ge from ab initio simulations Jeng-Da Chai

Dynamic structure factor of liquid and amorphous Ge from ab initio simulations

Jeng-Da Chai*and D. Stroud^†

Department of Physics, The Ohio State University, Columbus, Ohio 43210 J. Hafner and G. Kresse

Institut fu¨r Materialphysik and Center for Computational Materials Science, Technische Universita¨t Wien, Sensengasse 8/12, A-1090 Wien, Austria

!Received 23 June 2002; revised manuscript received 15 January 2003; published 28 March 2003"

We calculate the dynamic structure factor S(k,#) of liquid Ge (l-Ge" at temperature T!1250 K, and of amorphous Ge (a-Ge" at T!300 K, using ab initio molecular dynamics. 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 ultrasoft 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. $Phys. Rev. B 63, 134205 !2001"% 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. $Phys. Rev. Lett. 56, 1720 !1986"%. Our results suggest that the ab initio approach is sufficient to allow approximate calculations of S(k,#) in both liquid and amorphous materials.

DOI: 10.1103/PhysRevB.67.104205 PACS number!s": 61.20.Ja, 61.20.Gy, 61.20.Lc, 61.43.Dq

I. INTRODUCTION

Ge is a well-known semiconductor in its solid phase, but becomes metallic in its liquid phase. Liquid Ge (l-Ge" has, near its melting point, an electrical conductivity characteris- tic of a reasonably good metal (&1.6"10^#4'^#1cm^#11), 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 residual tetrahedral short- range order. This shoulder is absent in more conventional liquid metals such as Na or Al, which have more of a close- packed structure in the liquid state and a shoulderless first peak in the structure factor. Similarly, the bond-angle distri- bution function just above melting is believed to have peaks at two angles, one near 60° and characteristic of close pack- ing, and one near 108°, indicative of tetrahedral short-range order. This latter peak rapidly disappears with increasing temperature in the liquid state.

These striking properties of l-Ge have been studied theo- retically by several groups. Their methods fall into two broad classes: empirical and first principles. A typical empirical calculation is that of Yu et al.,² who calculate the structural properties of l-Ge assuming that the interatomic potentials in l-Ge are a sum of two-body and three-body potentials of the form proposed by Stillinger and Weber.³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 cor- respondingly smaller than their actual values.

In the second approach, the electronic degrees of freedom are taken explicitly into account. If the electron-ion interac- tion is sufficiently weak, it can be treated by linear-response theory.⁴ 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,⁴or unless electronic fluctuation forces are included.⁵ Such interactions are, however, in- cluded 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⁶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 coeffi- cients.

This combined approach, usually known as ab initio mo- lecular dynamics, was pioneered by Car and Parrinello,⁷and, in somewhat different form, has been applied to a wide range of liquid metals and alloys, including l-Ge,^{8 –10} l-GaxGe1#x,¹¹ stoichiometric III-V materials such as l-GaAs, l-GaP, and l-InP,^12,13 and nonstoichiometric l-GaxAs1#x,¹⁴ l-CdTe,¹⁵ and l-ZnTe,¹⁶ among other materi- als which are semiconducting in their solid phases. It has been employed to calculate a wide range of properties of these materials, including the static structure factor, bond- angle distribution function, single-particle electronic density of states, dc and ac electrical conductivity, and the ionic self- diffusion coefficient. The calculations generally agree quite well with data from available experiments.

A similar ab initio approach has also been applied exten- sively to a variety of amorphous semiconductors, usually ob- tained 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 fic-

titious classical variables" to obtain many structural and elec- tronic properties of amorphous Si.^17,18 A similar approach has been used by Lee and Chang.¹⁹Kresse and Hafner⁸ ob- tained both S(k) and g(r), as well as many electronic prop- erties, of a-Ge, using an ab initio approach similar to the one 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.,²⁰ also making use of a plane-wave basis and treating the electron-density functional in the generalized gradient approximation and other amorphous semiconductors have been carried out by !GGA".²¹More recently, a number of calculations for a-Si Sankey and Niklewsky,²² and by Drabold and collaborators.²³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.²⁴

Finally, we mention a third approach, intermediate be- tween empirical and ab initio molecular dynamics, generally known as tight-binding molecular dynamics. In this ap- proach, the electronic part of the total energy is described using a general tight-binding Hamiltonian for the band elec- trons. The hopping matrix elements depend on separation between the ions, and additional terms are included to ac- count for the various Coulomb energies in a consistent way.

The parameters can be fitted to ab initio calculations, and forces on the ions can be derived from the separation depen- dence of the hopping matrix elements. This approach has been used, e.g., to treat l-Si,²⁵ a-Si,²⁶ and liquid compound semiconductors such as l-GaAs and l-GaSb.²⁷Results are in quite good agreement with experiment.

In this paper, we extend the method of ab initio molecular dynamics !MD" 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. The present work appears to be the first to calculate S(k,#) using ab initio molecular dynamics. Here, we will calculate S(k,#) for l-Ge, where some recent experiments²⁸ 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 vibrational density of states of as-quenched a-Ge near tem- perature T!300 K in reasonable agreement with experi- ment. The calculated S(k,#) for the liquid also agrees quite well with experiment, especially considering the computa- tional uncertainties inherent in an ab initio simulation with its necessarily small number of atoms and limited time inter- vals.

The remainder of this paper is organized as follows. A brief review of the calculational method is given in Sec. II.

The results are presented in Sec. III, followed by a discussion and a summary of our conclusions in Sec. IV.

II. METHOD

Our method is similar to that described in several previ- ous papers,^10,11,14 but uses the Vienna ab initio simulation

package !VASP", whose workings have been extensively de- scribed in the literature.²⁹ 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 inte- grated 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. 10.

The VASP code uses ultrasoft Vanderbilt pseudopotentials,³⁰ a plane-wave basis for the wave func- tions, with the original Monkhorst-Pack (3"3"3) k-space meshes³¹and a total of 21 952 plane waves, corresponding to an energy cutoff of 104.4 eV. We use a finite-temperature version of the Kohn-Sham theory,³²in which the electron-gas Helmholtz free energy is calculated on each time step. This version also broadens the one-electron energy levels to case the k-space sums converge more rapidly. Most of our calcu- lations are done using the GGA !Ref. 21" for the exchange- correlation energy !we use the particular form of the GGA developed by Perdew and Wang²¹", but some are also carried out using the local-density approximation !LDA".

In our iteration of Newton’s laws in liquid Ge (l-Ge", we typically start from the diamond structure !at the experimen- tal liquid state density for the temperature of interest", then iterate for 901 time steps, each 10 fs, using the LDA. To obtain S(k) within the GGA, we start from the LDA configu- ration after 601 time steps, then iterate 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 at a time t2 after the start of the GGA simulation. We average over all t₂’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 t₁ after the start of this simulation;

the calculated LDA S(k) is also averaged over all t₁’s from 1.0 ps to 3.0 ps.

To calculate quantities for amorphous Ge (a-Ge", we start with Ge in the diamond structure at T!1600 K but at the calculated liquid density for that temperature, as given in Ref. 10. Next, we quench this sample to 300 K, cooling at a uniform rate, 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 at a time t1 after the system has reached 300 K; we also average this S(k) over all t₁’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 at a time t₂ 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 unrepresentative 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 de- fects 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 experiments on a-Ge themselves show some varia- tion, depending on the exact method of sample preparation.

We have used the procedure outlined above to calculate various properties of l-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.³³ However, our results for the dy- namic 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

! ⁽

"

0

#N+_k,0, !1"

where ri is the position of the ith ion 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 calcu- lations, we have used a cubic cell with N!64 and periodic boundary conditions in all three directions. The choices of particle number and cell shape are compatible with any pos- sible 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 "/0!k,t "0!#k,0"1dt, !2"

where the Fourier component0(k,t) of the number density is defined by

0!k,t "!

(

N

exp$#ik•ri!t "%. !3"

In our calculations, the average /. . .1 is computed as /0!k,t "0!Àk,0"1! 1

2t₁

#

0!k,t₁$t "0!Àk,t1"dt1 !4"

over a suitable range of initial times t₁. Because of the ex- pected 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

#

!5"

For large enough t₂, S(k,#,t₁,t₂) should become indepen- dent of t2 but will still retain some dependence on t1. There- fore, in the liquid, we obtain our calculated dynamic struc- ture factor, S_calc(k,#), by averaging over a suitable range of t1 from 0 to 2t1:

Scalc!k,#"! 1 2t₁

#

2t₁

dt1S!k,#,t1,t2". !6"

We choose our initial time in the t₁ integral to be 1 ps after the start of the GGA calculation, and !in the liquid" 2t1

!6 ps. We choose the final MD time t₂!16.41 ps. For a-Ge, we use the same procedure but t2!18.11 ps in our simulations. For our finite simulational sample, the calcu- lated 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 re- duces statistical error in both S(k,#) and S(k).

Finally, we have also incorporated the experimental reso- lution functions into our plotted values of S(k,#). Specifi- cally, we generally plot S_av(k,#)/S(k), where Sav(k,#) is obtained from the !already orientationally averaged" S(k,#) using the formula

S_av!k,#"!

#

.

R!###!"S!k,#!"d#!, !7"

where the resolution function R(#) $normalized so that 3_#.^. R(#)d#!1] is

R!#"! 1

!

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 nega- tive # for clarity, but they do not provide any additional information.

III. RESULTS A. S„k… for l-Ge and a-Ge

In Fig. 1, we show the calculated S(k) for l-Ge at T

!1250 K, as obtained using the procedure described in Sec.

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 pre-

vious simulations.^8,10Most 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 experi- mental results of Waseda et al.;³⁴ agreement between simu- lation 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 300 K. We prepare our sample of a-Ge as described in the previous section. For l-Ge, we average the calculated S(k) over different k vectors of the same length.

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 pre-

pared and the averages obtained as described in Sec. II. In contrast to l-Ge, but consistent with previous simulations,^8,24 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;³⁵in particular, the split prin- cipal peak seen in experiment is accurately reproduced by the simulations.

We have also calculated a number of other quantities for both l-Ge and a-Ge, including the pair distribution function g(r), and the electronic density of states n(E). For l-Ge, we calculated n(E) using the Monkhorst-Pack mesh with gamma point shifting !one of the meshes recommended in theVASPpackage". The resulting n(E) is generally similar to that found in previous calculations,^8,10 provided that an av- FIG. 1. !Color online" Static structure factor S(k) for l-Ge at T!1250 K, just above the ex- perimental melting temperature. Full curve: pre- sent 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. 34.

FIG. 2. !Color online" Full curve: Calculated S(k) for a-Ge at T!300 K, as obtained using the GGA. Structure is prepared as described in the text. Dashed curve: same as full curve, but calcu- lated in the LDA. Open circles: measured S(k) for a-Ge at T!300 K, as given in Ref. 35.

erage is taken over at least five to ten liquid state configura- tions. Our n(E) for a-Ge $calculated using a shorter averag- ing time than that used below for S(k,#)] is also similar to that found previously.⁸ Our calculated g(r)’s for both l-Ge and a-Ge, as given by theVASPprogram, are similar to those found in Refs. 8 and 10. 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 l-Ge, if we count as ‘‘nearest neighbors’’ all those atoms within 3.4 Å of the central atom !the larger of the cutoffs used in Ref. 8"

we find approximately 7.2 nearest neighbors, quite close to the value of 6.9 obtained in Ref. 8 for that cutoff. Finally, we have recalculated the self-diffusion coefficient D(T) for l-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. 10.

B. S„k,!… for l-Ge and a-Ge 1. l-Ge

Figure 3 shows the calculated ratio S(k,#)/S(k) for l-Ge at T!1250 K, as obtained using the averaging procedure described in Sec. II. We include a resolution function $Eqs.

!7" and !8"% of width 4#!2.5 meV, the same as the quoted experimental width.²⁸In Fig. 4, we show the same ratio, but without the resolution function !i.e., with #0!0). Obvi- ously, 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 l-Ge with hydrodynamic predictions, which should be appropriate at small k and#. This prediction takes the form !see, for example, Ref. 36"

2-^S!k,#"

S!k " !5#1

5

$

%

$1

5

$

$ 6k²

!##csk "²$!6k²"²

%

Here 5!c_P/cV is the ratio of specific heats, constant pres- sure, and constant volume; D_T is the thermal diffusivity; c_s is the adiabatic sound velocity; and 6 is the sound attenua- tion constant. D_Tand 6 can in turn be expressed in terms of other quantities. For example, DT!7_T/(0cP), where 7_T is the thermal conductivity and0 is the atomic number density.

Similarly, 6!¹₂$a(5#1)/5$b%, where a!7_T/(0c_V) and b is the kinematic longitudinal viscosity !see, for example, Ref.

36, pp. 264 –266".

FIG. 3. Calculated ratio of dynamic structure factor S(k,#) to static structure factor S(k) for l-Ge at T!1250 K 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 cal- culated 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 4#₀!2.5 meV, as in Eqs. !7" and !8". This value of #0

is chosen to equal the quoted experimental resolution !Ref. 28".

FIG. 4. Same as in Fig. 3, but without the resolution func- tion.

Equation !9" indicates that S(k,#) in the hydrodynamic regime should have two propagating peaks centered at #

!%c_sk, and a diffusive peak centered at #!0 of width determined by D_T. 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 4#&10 meV for k

!5.60 nm^#1, 4#&11 meV for k!7.92 nm^#1, and !some- what less clearly" 4#&13 meV for k!9.70 nm^#1. The value of c_s estimated from the lowest-k value is c_s&2.7

"10⁵ cm/sec. !The largest of these three k values may al- ready be outside the hydrodynamic, linear-dispersion regime."

These predictions agree reasonably well with the mea- sured S(k,#) obtained by Hosokawa et al.,²⁸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^#1occur at 4#!17.2 meV. Furthermore, the in- tegrated relative strength of our calculated sound-wave peaks, compared to that of the central diffusion peak, de- creases between k!7.92 nm^#1 and 12.5 nm^#1, consistent with both Eq. !9" and the change in experimental behavior²⁸ 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 D_T, 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 24#₀

!5.5 meV. Thus, a rough estimate of the intrinsic full width is 8

!

Using this estimate, one obtains D_T81.3"10^#3cm²/sec.

The hydrodynamic expression for S(k,#)/S(k) was origi- nally obtained without consideration of the electronic de- grees of freedom. Since l-Ge is a reasonably good metal, one might ask if the various coefficients appearing in Eq. !9"

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, and7_T, or from only the ionic contributions to these quantities? For l-Ge, the question is most relevant for7_T, since the dominant contribution to c_P and c_V should be the ionic parts, even in a liquid metal.⁴ However, the principal contribution to7_T is expected to be the electronic contribution.

We have made an order-of-magnitude estimate of DTus- ing the experimental liquid number density and the value C_P!(5/3)k_Bper ion, obtaining the electronic contribution to 7_T from the Wiedemann-Franz law³⁷ together with previ- ously calculated estimates of the electronic contribution.¹⁰ This procedure yields D_T&0.1 cm²/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 D_Twhich should be used in Eq. !9" for l-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 ef- fective pairwise screened ion-ion interactions.⁴ 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 l-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 sup- ported by our numerical results.

For k beyond around 12 nm^#1, the hydrodynamic model should start to break down, since the dimensionless param- eter#9!where9is the Maxwell viscoelastic relaxation time"

becomes comparable to unity. At these larger k’s, both our calculated and the measured²⁸ curves of S(k,#)/S(k) con- tinue 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.³⁸ 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.²⁸ The likely physics behind the de Gennes narrowing is straightforward. The half width of S(k,#) is inversely pro- portional 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 attempt to describe these fluctuations by a very oversimplified Langevin model.

We suppose that the Fourier component 0(k,t) $Eq. !3"% is governed by a Langevin equation

0˙ !k,t"!#:0!k,t "$;!t ". !10"

Here the dot is a time derivative,:is a constant, and;(t) is a random time-dependent ‘‘force’’ which has ensemble aver- age /;1!0 and correlation function /;(t);*(t!)1!A+(t

#t!). Equation !10" can be solved by standard methods !see, e.g., Ref. 36 for a related example", with the result !for suf- ficiently large t)

/0!k,t "0*!k,t$9"1!-A

: ^exp!#&9&:". !11"

According to Eq. !2", S(k,#) is, to within a constant factor, the frequency Fourier transform of this expression, i.e.,

S!k,#"<

#

or, on carrying out the integral,

S!k,#"< -A

#²$:^{2. !13"}

This is a Gaussian function centered at#!0 of half width:. On the other hand, the static structure factor

S!k "< lim

9→0

/0!k,t"0*!k,t$9"1!-A

: ^{. !14"}

Thus, if the constant A is independent of k, the half width: of the function S(k,#) at wave number k is inversely pro- portional to the static structure S(k). This prediction is con- sistent with the ‘‘de Gennes narrowing’’ seen in our simula- tions and in experiment.²⁸

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!300 K. 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 incorpo- rated a resolution function of width 4#0!2 meV into S(k,#). This width is a rough estimate for the experimental

resolution function in the measurements of Maley et al.;³⁹ 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³⁹

G!k,#"!

'

( ^$

^%

Here C is a k- and #-independent constant, and n(#,T)

!1/$e^4#/kBT#1% is the phonon occupation number for phonons of energy E!4# at temperature T. The quantity plotted by Maley et al.³⁹ is an average of G(k,#) over a range of k values from 40 to 70 nm^#1. These workers as- sumed that this average is proportional to the vibrational density of states n_vib(#). The measured n_vib(#) as obtained in this way³⁹is shown in Fig. 8 for two different amorphous structures, corresponding to two different methods of prepa- ration, and having differing degrees of disorder.

In order to compare our calculated S(k,#) to experiment, we use Eq. !15" to infer G(k,#), then average over a suit- able range of k. However, in using Eq. !15", we use the classical form of the occupation factor, n(4#)$1 8kBT/4#. This choice is justified because we have calcu- FIG. 5. Calculated S(k,#)/S(k) for a-Ge at T!300 K at k

=35 nm^#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 4#0

!2 meV. This value is chosen to give the best results for nvib(#) as measured by Ref. 39. The time step here is 10 fs.

FIG. 6. Same as Fig. 5, but at k>35 nm^#1.

lated S(k,#) using classical equations of motion for the ions.

We thus obtain for the calculated vibrational density of states

n_vib!#"8

'

( '

(

In Fig. 8 we show two such calculated plots of n_vib(#), obtained by averaging Eq. !16" over two separate groups of k’s, as indicated in the caption⁴⁰. For comparison, we also show n_vib(#) for a-Ge as calculated in Ref. 8 directly from the Fourier transform of the velocity-velocity autocorrelation function.

The calculated plots for n_vib(#) in Fig. 8 have some dis- tinct 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 simi- lar.

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 dif- ferences probably arise from several factors. First, there are several approximations involved in going from the calculated and measured S(k,#)’s to the corresponding n_vib(#)’s, and FIG. 7. !Color online" Calculated S(k,#)/S(k) as in Figs. 5 and 6 but including higher frequencies #. At 4#!60 meV, the curves are arranged vertically in order of increas- ing frequency. Each curve is vertically displaced by 0.0005 units from the one below it.

FIG. 8. !Color online" Full curve: calculated vibrational density of states n_vib(#), in units of 10^#3 states/meV. n_vib(#) is obtained from the resolution broadened S(k,#) and S(k) of Fig. 7 using the formula n_vib(#)!/^Gcalc(k,#)1^{, where} G_calc(k,#) is given by the right-hand side of Eq.

!16", and the averaging is carried out over the three magnitudes of k near 40 nm^#1for 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 n_vib(#), as re- ported in Ref. 39 for two forms of a-Ge. Finally, the open diamonds denote n_vib(#) as calculated in Ref. 8 from the ionic velocity-velocity autocor- relation function !dot-dashed curves". In all plots except that of Ref. 8, n_vib(#) is normalized so that 3₀^#maxn(#)d(4#)!1. #_maxis the frequency at which n_vib(#)→0, and is estimated from this figure by extrapolating the right-hand parts of the solid and dashed curves linearly to zero.

these may be responsible for some of the discrepancies. Sec- ondly, there may actually be differences between the particu- lar 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. Consid- ering all these factors, agreement between calculation and experiment is quite reasonable.

Previous ab initio calculations for a-Ge !Ref. 8" 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. 8 do not require computing S(k,#). In the present work, by contrast, we start from S(k,#) !which is calculated here in an ab initio calculation for a-Ge", and we work backwards to get n_vib(#). 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 po- sitions. In Fig. 8, we also show the results of Ref. 8 for n_vib(#) 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 assum- ing that the quenched configuration is a local energy mini- mum and calculating the potential energy for small positional deviations from that minimum using ab initio molecular dy- namics. This procedure has been followed for a-GeSe2, for example, by Cappelletti et al.⁴¹These workers have then ob- tained 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 does n_vib(#).

IV. DISCUSSION AND CONCLUSIONS

The present results show that ab initio molecular dynam- ics can be used to calculate the dynamic structure factor S(k,#) for both liquid and amorphous semiconductors. Al- though the accuracy of the calculated S(k,#) is lower than that attained for static quantities, 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 l-Ge agrees qualitatively with that measured by inelastic x-ray scattering,²⁸and that the one calculated for a-Ge leads to a vibrational density of states qualitatively similar to the quoted experimental one.³⁹ Since such calcu- lations are thus shown to be feasible, our work 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 underlying processes at the atomic scale which give rise to specific fea- tures in the measured and calculated S(k,#).

ACKNOWLEDGMENTS

This work has been supported by NASA, Division of Mi- crogravity Sciences, through Grant No. NCC8-152, and by NSF Grant Nos. DMR01-04987 !D.S." and CHE 01-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. We are very grateful to Professor David H.

Matthiesen for his continual support and encouragement through the course of this work, and to Professor Neil Ash- croft and Sergey Barabash for valuable conversations.

Jeng-Da Chai wishes to thank Lan Bi for her endless support.

*Present address: Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742. Email address: [email protected]

†Corresponding author. FAX: !614"292-7557. Email address:

[email protected]

1V. M. Glazov, S. N. Chizhevskaya, and N. N. Glagoleva, Liquid Semiconductors !Plenum, New York, 1969".

2W. B. Yu, Z. Q. Wang, and D. Stroud, Phys. Rev. B 54, 13 946

!1996".

3F. H. Stillinger and Weber, Phys. Rev. B 31, 5262 !1985".

4For a review, see, e. g., N. W. Ashcroft and D. Stroud, in Solid State Physics !Academic, New York, 1978", Vol. 33, p. 1.

5N. W. Ashcroft, Nuovo Cimento D 12, 597 !1990".

6W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 !1965".

7R. Car and M. Parrinello, Phys. Rev. Lett. 35, 2471 !1985".

8G. Kresse and J. Hafner, Phys. Rev. B 49, 14 251 !1994".

9N. Takeuchi and I. L. Garzo´n, Phys. Rev. B 50, 8342 !1995".

10R. V. Kulkarni, W. G. Aulbur, and D. Stroud, Phys. Rev. B 55, 6896 !1997".

11R. V. Kulkarni and D. Stroud, Phys. Rev. B 57, 10 476 !1998".

12Q. Zhang, G. Chiarotti, A. Selloni, R. Car, and M. Parrinello,

Phys. Rev. B 42, 5071 !1990".

13L. Lewis, A. De Vita, and R. Car, Phys. Rev. B 57, 1594 !1998".

14R. V. Kulkarni and D. Stroud, Phys. Rev. B 62, 4991 !2000".

15V. Godlevsky, J. Derby, and J. Chelikowsky, Phys. Rev. Lett. 81, 4959 !1998"; V. Godlevsky, M. Jain, J. Derby, and J. Che- likowsky, Phys. Rev. B 60, 8640 !1999".

16M. Jain, V. V. Godlevsky, J. J. Derby, and J. R. Chelikowsky, Phys. Rev. B 65, 035212 !2001".

17R. Car and M. Parrinello, Phys. Rev. Lett. 63, 204 !1988".

18I. Stich, R. Car, and M. Parrinello, Phys. Rev. Lett. 63, 2240

!1989"; Phys. Rev. B 44, 4261 !1991"; ibid. 44, 11 092 !1991".

19I. Lee and K. J. Chang, Phys. Rev. B 50, 18 083 !1994".

20N. C. Cooper, C. M. Goringe, and D. R. McKenzie, Comput.

Mater. Sci. 17, 1 !2000".

21See, e. g., J. P. Perdew, in Electronic Structure of Solids, edited by P. Ziesche and H. Eschrig !Akademic Verlag, Berlin, 1991".

22O. Sankey and D. J. Niklewsky, Phys. Rev. B 40, 3979 !1989".

23D. A. Drabold, P. A. Fedders, O. F. Sankey,and J. D. Dow, Phys.

Rev. B 42, 5135 !1990".

24F. Alvarez, C. C. Dı´az, A. A. Valladares, and R. M. Valladares, Phys. Rev. B 65, 113108 !2002".

25C. Z. Wang, C. T. Chan, and K. M. Ho, Phys. Rev. B 45, 12 227

!1992"; I. Kwon, R. Biswas, C. Z. Wang, K. M. Ho, and C. M.

Soukoulis, ibid. 49, 7242 !1994".

26G. Servalli and L. Colombo, Europhys. Lett. 22, 107 !1993".

27C. Molteni, L. Colombo, and Miglio, J. Phys.: Condens. Matter 6, 5255 !1994".

28S. Hosokawa, Y. Kawakita, W.-C. Pilgrim, and H. Sinn, Phys.

Rev. B 63, 134205 !2001".

29G. Kresse and J. Hafner, Phys. Rev. B 47, 558 !1993"; G. Kresse, Ph.D. Technische Universita¨t Wien, 1993; G. Kresse and J.

Furthmu¨ller, Comput. Mater. Sci. 6, 15 !1996"; G. Kresse and J.

Furthmu¨ller, Phys. Rev. B 54, 11 169 !1996".

30D. Vanderbilt, Phys. Rev. B 32, 8412 !1985".

31H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 !1976".

32N. D. Mermin, Phys. Rev. 137, A1141 !1965".

33Our method is very similar to that used in Ref. 8 to treat l-Ge and a-Ge, but with the following differences: !i" we use both the GGA and the LDA, while Ref. 8 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. 8. Our a-Ge structure may also differ slightly from that of Ref. 8, since we use a somewhat faster quench rate (&3"10¹⁴K/sec rather than 1.7

"10¹⁴K/sec). Both codes use the conjugate gradient method to

find energy minima and the same Vanderbilt ultrasoft pseudopo- tentials.

34Y. Waseda, The Structure of Noncrystalline Materials: Liquids and Amorphous Solids !McGraw-Hill, New York, 1980".

35G. Etherington, A. C. Wright, J. T. Wenzel, J. T. Dore, J. H.

Clarke, and R. N. Sinclair, J. Non-Cryst. Solids 48, 265 !1982".

36See, for example, J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, 2nd ed. !Academic Press, London, 1986", pp.

222–228.

37See, for example, N. W. Ashcroft and N. D. Mermin, in Solid State Physics !Harcourt-Brace, Fort Worth, 1976", Chap. 2.

38P. G. de Gennes, Physics !Long Island City, N.Y." 25, 825 !1959";

ibid. 3, 37 !1967".

39N. Maley, J. S. Lannin,and D. L. Price, Phys. Rev. Lett. 56, 1720

!1986".

40Two recent groups $V. Martı´n-Mayor, M. Me´zard, G. Parisi, and P. Verrocchio, J. Chem. Phys. 114, 8068 !2001"; T. S. Grigera, V.

Martı´n-Mayor, G. Parisi, and P. Verrocchio, Phys. Rev. Lett. 87, 085502 !2001"% have derived approximate expressions connect- ing S(k,#) to n_vib(#) in an amorphous solid, which closely resemble Eq. !16". In fact, their expressions can be shown to become equivalent to Eq. !16" at sufficiently large k.

41R. L. Cappelletti, M. Cobb, D. A. Drabold, and W. A. Kamitaka- hara, Phys. Rev. B 52, 9133 !1995".