De Gennes Narrowing and Hard-Sphere Approach - The Journal of

Physical Chemistry, Georg-August-University of Göttingen, Tammannstrasse 6, D-37077 Göttingen, Germany. Heinz Maier-Leibnitz Zentrum, Lichtenber...
10 downloads 9 Views 2MB Size
Article pubs.acs.org/JPCB

De Gennes Narrowing and Hard-Sphere Approach Oleg Sobolev* Institute for Physical Chemistry, Georg-August-University of Göttingen, Tammannstrasse 6, D-37077 Göttingen, Germany Heinz Maier-Leibnitz Zentrum, Lichtenbergstrasse 1, D-85748 Garching, Germany ABSTRACT: The energy width Δω of the quasielastic coherent dynamic structure factor S(Q, ω) for a simple liquid exhibits the oscillating dependence on wavenumber Q with the sharp minimum at Qmax corresponding to the maximum of the structure factor S(Q). The only known expression for Δω(Q) was derived for a dense hard-sphere (HS) fluid (Cohen et al., Phys. Rev. Lett. 1987, 59, 2872). Though this expression has been frequently used for the analysis of the experimental data obtained for liquid metals, until now, it has never been tested against a true HS fluid. A test performed by means of HS molecular dynamic simulations reveals a considerable discrepancy between the simulations results and the examined model. The main output of the analysis is the finding that the ΔωHS(Q) behavior is defined in terms of the average cage size, ⟨Lc⟩, rather than of the HS diameter, σHS. The simulated ΔωHS(Q) has been compared with the results for the soft-spherical potential. The microscopic dynamics of the soft-sphere fluid shows significant difference in comparison to the HS system. Nevertheless, the diffusive mobility of soft spheres can be characterized within the HS approximation using an effective diameter, σeff, and this parameter can be found from Δω(Q) at Q ≈ Qmax. A similar result has been obtained for the neutron scattering data measured for liquid Rb.

1. INTRODUCTION The self-diffusion coefficient, D, is a quantity of fundamental importance not only for the liquid state physics but also for a great number of technological applications. It can be obtained from incoherent quasielastic neutron scattering (QENS) data using the simple diffusion model, which has often been applied for liquid metals and alloys consisting of incoherent scattering atoms.1−6 The situation is more complicated in the case of liquids that scatter neutrons mostly coherently. Though the coherent QENS contains rich information on individual and cooperative diffusive motions of molecules in a sample, the set of analytical expressions for the detailed analysis of the corresponding experimental data is very limited, that brings molecular dynamics (MD) simulation methods to the fore. The dynamical structure factor S(Q,ω) measured in an inelastic neutron experiment is the distribution of neutrons that have undergone an energy exchange ℏω = Ei − Ef, and a wave vector transfer, Q = ki − kf, after scattering by the sample. S(Q,ω) is the Fourier transform of the intermediate scattering function I(Q,t) 1 S(Q, ω) = ∫ I(Q, t ) exp(−iωt ) dt (1) 2π In the case of coherent scattering on a monatomic sample I(Q,ω) reflects the pair correlations between the respective positions ri(0) and rj(t) of the scattering particles i and j at successive times 0 and t 1 I (Q , t ) = ∑ ∑ ⟨exp(−iQri(0)) exp(iQrj(t ))⟩ N i j (2) © 2016 American Chemical Society

In a liquid, the diffusive motion of molecules results in the loss of the correlation between ri(0) and rj(t). Due to this, I(Q,t) decays with time, and S(Q,ω) shows correspondingly the quasielastic broadening. De Gennes7 found that the normalized second frequency moment of the coherent dynamical structure factor S(Q,ω) for a classical simple fluid is in inverse proportion to the static structure factor S(Q) ⟨ωQ2 ⟩ =

kBT Q2 mS(Q )

(3)

where kB is Boltzmann’s constant, T is temperature of the fluid, and m is the mass of the fluid particle. The quantity ⟨ω2Q⟩1/2 gives an approximate measure of the spectral width of S(Q,ω). The presence of S(Q) at the denominator means that the line width of S(Q,ω) exhibits the oscillating Q-dependence with the sharp minimum at Qmax corresponding to the maximum of S(Q). This minimum, known as de Gennes narrowing, originates from the strong spatial correlations between the nearest-neighboring fluid particles and is an indication of the sensitivity of the relaxation processes with respect to the fluid structure. The only known analytical expression for the coherent quasielastic line width Δω (half width at half-maximum) of S(Q,ω) was derived approximately for a dense hard-sphere (HS) fluid within the framework of the extended heat mode approximation of Enskog’s kinetic theory8−12 Received: May 9, 2016 Revised: September 1, 2016 Published: September 1, 2016 9969

DOI: 10.1021/acs.jpcb.6b04685 J. Phys. Chem. B 2016, 120, 9969−9977

Article

The Journal of Physical Chemistry B ΔωHS

D Q2 ≈ E d(QσHS) S(Q )

d(QσHS) =

1 [1 − j0 (QσHS) + 2j2 (QσHS)]

collected for noble gas fluids12 and liquid metals,19−30 a numerical simulation is more suitable for such a test because, even for simple liquids, σeff is an ill-defined quantity. On the other hand, one cannot expect that the HS model can be applied to real liquids without any limitations. In order to explore these limitations, as the first step, a comparison of ΔωHS(Q) with the simulation result for the soft-spherical potential shall be done. The paper presents a test of eq 4 performed by means of HS molecular dynamic simulations and the analysis of the simulation results for the soft-sphere system in terms of the HS approach. Also, the neutron scattering data reported for liquid rubidium31 were analyzed using the HS approximation. The Q-range for the hard sphere and the soft-sphere simulations was restricted by the first maximum of S(Q). This was dictated by the practical focus of the work, because Δω(Q) can be measured with good accuracy at Q values that are not far from Qmax.

(4)

(4a)

where j0 and j2 are the first two even-order Bessel spherical functions. Enskog’s diffusion coefficient DE can be expressed13 in terms of the packing fraction φ = πnσHS3/6 DE =

1 σHS 16 φ

πkBT (1 − φ)3 m (1 − φ /2)

(5)

where σHS and n are the hard-sphere diameter and the atomic density, correspondingly. Equation 5 takes into account only binary collisions. The neglect of correlated motions results in significant deviations of calculated DE from the macroscopic diffusion constants. The MD simulations14 provided the ratio between the real hard-sphere diffusion coefficient and the value predicted by eq 5 as a function of the reduced density V0/V = nσHS3/√2

2. SIMULATION DETAILS The detailed description of the event driven algorithm used for the HS simulations is given in ref 32. The term “event driven” means that elastic collisions between HS particles are considered in chronological order using the calendar of forthcoming events updated by the program. This method is different from those used for continuous interaction potentials, when the dynamical equations are integrated numerically with a constant time step. The usual dimensionless units were used for the HS simulations and representation of the results. The hard-sphere diameter σHS = 1 and the mass of fluid particles m = 1 are the units of length and mass, respectively, and the unit of time is σHS(m/kBT)1/2. Also, the temperature kBT was set to unity. The simulated system consisted of 10 976 particles in a cubic cell with periodic boundary conditions. The initial configuration was obtained by melting the FCC structure during the equilibration run of ∼2 × 107 collisions (∼4000 collisions per particle). The size of the cell, L, was varied to give the desired packing fraction in the range φ = 0.222−0.49 (L = 22.72σ − 29.58σ). The duration of the production run was typically ∼104 collisions per particle, that was enough to perform good time averaging of I(Q,t). According to ref 33, no more than a few hundred collisions are needed for complete decay of I(Q,t) near Qmax. The coordinates and velocities of the particles were recorded with the time step Δt = 0.02. The Weeks-Chandler-Andersen34 (WCA) potential was used for the soft-sphere simulations

⎛ V V2 V 3⎞ DHS = DE⎜1 + 0.05403 0 + 6.3556 0 − 10.9425 0 ⎟ V V V ⎠ ⎝ (6)

It was demonstrated that the HS approach can be used for the calculation of the diffusion constants D in liquid metals.15−18 One of the problems in this approach lies in determining an effective HS diameter, σeff, that can approximate the repulsive part of a real pair potential. Since σeff cannot be estimated unambiguously using other methods, it is interesting to found out whether this parameter can be properly determined from the experimentally measured Δω(Q). Recently, eq 4 was employed to extract the diffusion coefficient of liquid aluminum in the temperature range of 943−1193 K.19 The parameter σeff derived from a Percus− Yevick fit to the main peak of S(Q) was used as an input parameter. Though the obtained values agree well with the viscosity data and with ab initio calculations, the practical applicability of the method proposed in ref 19 looks somewhat questionable. Indeed, if the σeff value is known, diffusion constants can be calculated directly using eqs 5 and 6 without the need to perform a neutron scattering experiment. On the other hand, having reliable experimental S(Q) data one can try to obtain the parameter σeff from the fit of eq 4 to the experimental Δω(Q). This has only been done for liquid Ga;20−22 that is, however, a priori a non-hard-sphere liquid. A very good agreement with inelastic X-ray scattering data measured near the melting point was demonstrated using the diameter σeff and turned out to be larger than the one that can produce a reasonable DE value using eq 5. A similar result was obtained in the neutron scattering experiment.23 This caused a discussion21,22 on the physical meaning of the parameter σeff when it was used in eq 4. Recently, the analysis of the Δω(Q) for liquid Ga measured in a wide temperature range was reported.24 The fit of eq 4 to the experimental Δω(Q) showed an increase of σeff with temperature from σeff = 2.8 Å at T = 313 K up to σeff = 3.03 Å at T = 793 K. The latter result also looks meaningless from the point of view of the HS approach. Any speculations based on the fit of eq 4 to experimental data are not very meaningful, if it is not known how precisely eq 4 can describe the ΔωHS(Q) for a true HS system. Although eq 4 has been repeatedly tested against neutron scattering data

⎧ ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ ⎪ ⎢⎜ σ ⎟ σ ⎥ 1/6 ⎪ 4ε⎢⎜ ⎟ − ⎜⎜ ⎟⎟ ⎥ + ε rij < 2 σ ⎝ rij ⎠ ⎦ U (rij) = ⎨ ⎣⎝ rij ⎠ ⎪ ⎪ 0 rij ≥ 21/6σ ⎩

(7)

The WCA simulations were performed using DL_POLY35 code, with a time step of 0.943 × 10−4 in reduced units of time (σ(m/ε)1/2). The parameter m represents the particle mass. The cubic simulation box consisted of 21 952 soft spheres. During an equilibration period of 1.5 × 106 time steps the velocities of the particles were rescaled every 50 steps in order to maintain the desired temperature. After equilibration, the simulation was performed in the NVE ensemble. The 9970

DOI: 10.1021/acs.jpcb.6b04685 J. Phys. Chem. B 2016, 120, 9969−9977

Article

The Journal of Physical Chemistry B coordinates and velocities of the particles were recorded every 100 time steps during the run of 106 steps. The WCA system was simulated in the range of the densities nσ3 between 0.424 (L = 37.27σ) and 0.917 (L = 28.82σ) for the temperature T* = kBT/ε = 2/3 and in the temperature range T* = 2/3 − 11.3 for nσ3 = 0.917. The recorded trajectory files for HS and WCA simulations were analyzed using the nMOLDIN package.36 The selected Qrange for the computation of I(Q,t) covered the whole region corresponding to the first maximum of S(Q).

3. RESULTS AND DISCUSSION 3.1. Hard-Sphere Fluid. Before exploring specifically I(Q,t), some structural and dynamical properties, such as S(Q) and DHS, are worth consideration. This is helpful for subsequent analysis, especially for the analysis of the WCA simulations. The Percus−Yevick (PY) equation for the HS structure factor is given as follows:37,38 SPY(x) = 1/[1 − C PY(x , φ)]

C PY(x , φ) = −

Figure 1. Simulated structure factor (points) compared with calculations according to eq 8 (dashed line) and eq 9 (solid line).

(8)

24φ {αx 3(sin x − x cos x) x6

+ ζx 2(2x sin x − (x 2 − 2)cos x − 2) + γ[(4x 3 − 24x)sin x − (x 4 − 12x 2 + 24) cos x + 24]} x = QσHS

(1 + 2φ)2 (1 − φ)4

α=

2

ζ = −6φ

(1 + φ /2) (1 − φ)4

γ=

Figure 2. DHS obtained in the present work (filled circles) compared with the data of Heyes42 (open circles). The line represents the polynomial fit used in eq 10.

1 φα 2

It is well-known that the PY approximation becomes increasingly inaccurate for φ > 0.4. Several semiempirical corrections to the PY equation were proposed by different authors.39−41 The most convenient for practical use is the “Percus-Yevick plus tail” (PY-T) expression:41 SPYT(x) = 1/[1 − C PYT(x , φ)]

C PYT(x , φ) = aC PY(x , φ) +

(9)

⎞ 24φb ⎛ d ⎜cos x + sin x⎟ 2 2⎝ ⎠ x x +d

Explicit formulas for the parameters a, b, and d are given in ref 41. Figure 1 shows that the PY-T approximation, indeed, describes better the simulated S(Q) than eq 8. The self-diffusion coefficients, DHS, for the HS fluid were obtained by integration of the velocity autocorrelation function (VACF) ψ(t) = 1/3⟨v(⃗ 0)v(⃗ t)⟩. They show a good agreement with very accurate and detailed simulation data of Heyes42 (Figure 2). Therefore, for further analysis, DHS will be calculated as follows: DHS = σHS kBT /m D(φ)

Figure 3. Velocity correlation functions for different packing fractions: φ = 0.37 (dashed line), φ = 0.426 (dotted line), φ = 0.444 (dashdotted line), and φ = 0.49 (solid line).

forward motion and results in a higher diffusion coefficient in comparison to the Enskog’s value DE.43−45 The diffusive dynamics at φ = 0.49 is dominated by the cage effect, which is indicated by a negative part of the ψ(t) and reduced the D/DE ratio. In this case, the diffusing particle is reflected with high probability back after a collision with the neighbors, and the diffusive motion is coupled to the density fluctuations.13,46 The ψ(t) for φ = 0.426 and φ = 0.444 are in a transition zone: they show a pronounced ″cage″ minimum and a small, but clearly visible, positive tail. According to eq 6, DHS = DE at φ ≈ 0.435. The shape of the simulated I(Q,t) was taken into account for the proper determination of the line width Δω (see Figure 4).

(10)

where D(φ) is a polynomial fit of the data tabulated in ref 42. The simulated VACF plotted in Figure 3 reproduced the pioneering results of Alder and Wainwright43−45 and illustrated the different microscopic dynamics regimes realized at different densities. At lower density (φ = 0.37) the ψ(t) shows a longtime positive tail caused by a vortex flow pattern generated by the diffusing particle in its neighborhood that supports the 9971

DOI: 10.1021/acs.jpcb.6b04685 J. Phys. Chem. B 2016, 120, 9969−9977

Article

The Journal of Physical Chemistry B

Figure 5. ΔωHS(Q) as a function of QσHS for different φ. The points depict the simulation results and the curves are calculated according to eq 4.

φ at the right side of Qmax, where eq 4 strongly overestimates the ΔωHS values. If we seek an expression that can describe ΔωHS(Q) better than eq 4, eq 3 gives good reasons to search for a solution in a form similar to eq 4

Figure 4. Simulated I(Q,t) for QσHS = 6 (points) fitted by the KWW function (lines). The shape parameter β increases with decreasing φ.

The most convenient function for this purpose is Kohlrausch− Williams−Watts (KWW) function KWW(t ) = exp( −(t /τw)β )

ΔωHS =

(11)

with the mean relaxation time ⟨τw⟩ defined as ∞ τ ⟨τw ⟩ = dt KWW(t ) = w Γ(β −1) β 0



(14)

One can expect that the correction function h(QσHS,φ) for the HS fluid does not depend on temperature, but depends on the density. The simulation results shown in Figure 6 confirm

the ΔωHS(Q) was determined as follows Δω = 1/⟨τW ⟩

DHSQ 2 h(QσHS , φ) S(Q )

(12)

(13)

I(Q,t) demonstrates the stretched exponential decay (β < 1) at high φ, where cage diffusion dominates. At the same time, the shape parameter β goes up above unity at lower φ. The nonexponential decay of I(Q,t) is not unusual for simple liquids. For example, the stretched exponential behavior was observed in liquid Rb near the melting point.31 In the free particle limit, I(Q,t) is a Gaussian both in the wavevector and in time exp(−kBTQ2t2/2m).33 Therefore, one should expect that I(Q,t) becomes more Gaussian-like (β > 1) at lower density due to increasing free-path lengths between collisions, l. Such transition from a Lorentzian-like (β ≈ 1) to a more Gaussianlike shape of the incoherent S inc (Q,ω) was detected experimentally for liquid Na.47 According to ref 47, this transformation starts at Q values far below those corresponding to the free particle limit, and this effect is more pronounced at higher temperatures (lower densities). Although the coherent scattering represents a more complicated case, a qualitatively similar behavior can be expected, because, at large Ql, the pair correlations tend to vanish, and the self-correlation part of I(Q,t) becomes dominant. The linewidths, ΔωHS, obtained for different φ are shown in Figure 5. Equation 4 is consistent with the simulated data in the φ-range from 0.37 to 0.46 but only for the limited Q-interval at the left side of Qmax. The model curves calculated for φ > 0.47 and φ < 0.33 are above the simulated points in the whole examined Q-region. The discrepancy is especially strong for all

Figure 6. Ratio ΔωHS(Q)S(Q)/DQ2 (connected circles) for different φ. The red line represents eq 4a.

this expectation. Not only the intensity of h(QσHS,φ), but also its peak position, depends on φ and coincides approximately with Qmax. This finding contrasts with eq 4a that supposes for d(QσHS) the maximum fixed at QσHS ≈ 7.5 and not related to φ. In addition, this QσHS value is larger than the values for the maxima of h(QσHS,φ) (Figure 6). Therefore, if one tries to construct h(QσHS,φ), as an empirical modification of eq 4a, one has to introduce a density dependent coefficient b(φ) that shifts the maximum of h(QσHS,φ) to the left in the Q-scale. Also, it would be surprising if DHS (or DE) could stay in the numerator of eq 14 9972

DOI: 10.1021/acs.jpcb.6b04685 J. Phys. Chem. B 2016, 120, 9969−9977

Article

The Journal of Physical Chemistry B without any correction factors. Thus, the simplest upgrade of eq 4a can be designed as follows: h(QσHS , φ) =

a (φ ) [1 − j0 (b(φ)QσHS) + 2j2 (b(φ)QσHS)] (14a)

Equation 14a fits the simulation results well in the φ-region where the cage effect dominates the diffusion process, but the quality of the fit degrades rapidly at φ < 0.43 (Figure 7). This

Figure 8. Values for a(φ) and b(φ) (symbols) found from the fit of eq 14a to the simulated h(QσHS,φ). The solid lines represent eq 14b. The green dashed line shows (σHS+lE)/σHS (eq 15). The vertical dashed line (φ = 0.435) depicts an approximate borderline between different regimes of diffusion, cage diffusion, and vortex diffusion.

(where lE is the Enskog’s mean free path between collisions) due to the nonvanishing probability of forward motion after a collision. Figure 8 shows that the parameter b(φ), indeed, is very close to the ratio σHS + lE (1 − φ)3 +1 = 6 2 (1 − φ /2)φ σHS

(15)

The present analysis is basically phenomenological, and it is focused on the Q-region of the main peak of S(Q). Thus we cannot expect that eq 14 can describe the second minimum of ΔωHS(Q) as well as the first one (Figure 9). Moreover, eq 14

Figure 7. Fit of eq 14a (lines) to the results of the HS simulations (circles). The left panel shows the data for the φ-region where the cage effect dominates the diffusion process. The right panel corresponds to the vortex diffusion region.

fact is understandable, because one cannot expect that this expression can take into account both mechanisms: cage diffusion and vortex diffusion. Note that the initial model (eq 4) was designed for describing cage diffusion in a dense fluid.12 The results of the fit, the values for a(φ) and b(φ) for the high densities (φ > 0.435), are shown in Figure 8. These two parameters can be expressed, as a function of φ, using the following empirical formulas:

Figure 9. ΔωHS(Q) for φ = 0.444 in the Q-range corresponding to the first two maxima of S(Q): simulation results (circles); eq 14 (solid line); eq 4 (dashed line).

(as well as eq 4) has a limited ability to reproduce the ΔωHS(Q) behavior at large Q, where pair correlation effects tend to vanish, and the difference between Δω(Q) and incoherent line width Δωinc(Q) becomes negligible. It is known that Δωinc(Q) deviates strongly from the DQ2 dependence in this Qregion,13,46 whereas the discussed formula predicts the DHSQ2 limit for ΔωHS(Q) at Q →∞. Therefore, in spite of the fact that eq 14 represents a substantial improvement of eq 4, it cannot be considered the final version of the expression for ΔωHS(Q). In order to calculate ΔωHS(Q) with arbitrary parameters σHS, m, and T using the HS simulation data (obtained for σHS = 1, m = 1, and kBT = 1), it is convenient to collect all φ-dependent terms in one function and to express ΔωHS(Q) as follows:

a(φ) = 7.658 − 34.835φ + 44.08φ 2 b(φ) = 1.323 − 0.552φ

(14b)

Note that, in contrast to eq 4, which supposes the dominance of the stochastic binary collisions, eq 14 contains DHS instead of DE in the numerator. It is easy to see that the values found for a(φ)DHS are smaller than the corresponding values for DE by 16−33%. The decrease of b(φ) with increasing φ suggests that the quantity ⟨Lc⟩ = b(φ)σHS should be considered a characteristic length related to the average cage size occupied by the particle. At high densities, ⟨Lc⟩ should be slightly larger than σHS + lE 9973

DOI: 10.1021/acs.jpcb.6b04685 J. Phys. Chem. B 2016, 120, 9969−9977

Article

The Journal of Physical Chemistry B ΔωHS(Q ) = Q 2σHS

kBT HHS(QσHS , φ) m

(16)

The HS simulations results represented in the form of the function HHS(QσHS,φ) (Figure 10) will be used for the analysis of the WCA simulations and experimental data.

Figure 11. DWCA for different nσ3 as a function of σeff = πnσeff3/6 (left) and DWCA for nσ3 = 0.917 at different temperatures (right). The filled circles represent the results of the present work; open circles show the data of Heyes.42 The lines are DHS calculated using eq 10 with σHS = σeff given by eq 18.

Figure 10. Results of the HS simulations represented in the form of HHS(QσHS,φ).

3.2. Soft-Sphere Fluid. An attempt to describe structural and dynamical properties of the WCA fluid in terms of the HS approximation reveals some inconsistency of this approach. A temperature-rescaled effective diameter, σeff, of the WCA potential that accounts very well for the temperature dependence of the compressibility factor was proposed in ref 48 ⎞1/6 ⎛ 2 σeff = σ ⎜ ⎟ ⎝ 1 + T* ⎠

Figure 12. Simulated S(Q) for the WCA fluid at nσ3 = 0.917 (circles) compared with calculations according eq 9 using the parameter σeff given by eq 17 (dash−dot line) or by eq 18 (solid line).

(17)

It was already shown that if σeff is given by eq 17, the selfdiffusion coefficient of the WCA fluid, DWCA, is larger than DHS at the same effective packing fraction φeff = πnσ3eff/6.42 Equation 17 calculates σeff as the distance at which the WCA potential equals the thermal energy kBT. If one assumes that the potential U(rij) should be equal to the average kinetic energy of the particle 3kBT/2 at rij = σeff, the diameter σeff has to be defined as follows: ⎛ ⎞1/6 2 ⎜ ⎟ σeff = σ ⎜ ⎜ 1 + 3 T * ⎟⎟ ⎝ ⎠ 2

The linewidths, ΔωWCA, obtained for the WCA fluid at different densities and temperatures are compared with the results of the HS simulations, ΔωHS, in Figure 13 and Figure 14. The curves for ΔωHS were calculated according to eq 16 with σHS = σeff and φ = σeff; the parameter σeff was determined by eq 18. The values of ΔωWCA are found to be larger than those of ΔωHS in the whole simulated Q-range except the narrow region at Q ≈ Qmax where both curves coincide. The difference between ΔωWCA and ΔωHS becomes smaller with decreasing σeff (Figure 13). In addition, ΔωWCA must tend to its HS limit ΔωHS at T* → 0 (which means either T → 0 or ε → ∞). This dependence is illustrated in Figure 15. According to eq 18, the simulation results for nσ3 = 0.917 at T* = 2.35 and 11.3 are related to σeff = 0.4 and σeff = 0.3, correspondingly; they can be compared with the ΔωWCA obtained for the same values of σeff at T* = 2/3. In order to do a correct comparison, the data should be properly scaled taking into account the different values of σeff. Also, as has just been demonstrated, ΔωWCA(Q) = ΔωHS(Q) at Q ≈ Qmax. With eq 16 kept in mind for ΔωHS, the simulation results

(18)

Figure 11 shows that eq 10 with σHS = σeff determined according to eq 18 describes the values for DWCA obtained in the present work and the results of ref 42 with very good accuracy. On the other hand, the use of the parameter σeff given by eq 17 results in much better agreement between eq 9 and the simulated SWCA(Q) than the use of σeff defined by eq 18 (Figure 12). Thus even for WCA fluid, the σeff is a quantity that is meaningless without specifying a physical property it accounts for. Since we are interested in diffusion mobility of WCA particles, eq 18 shall be used for the further analysis. 9974

DOI: 10.1021/acs.jpcb.6b04685 J. Phys. Chem. B 2016, 120, 9969−9977

Article

The Journal of Physical Chemistry B

Therefore, the results of comparison between ΔωWCA(Q) and ΔωHS(Q) can be summarized by analogy with eq 16 using the following set of expressions: Δω WCA (Q ) = Q 2σeff

kBT HWCA(Qσeff , φeff , T *) m

(19)

In contrast to HHS(QσHS,φ), HWCA(Qσeff,φeff,T*) depends on temperature and the softness of the potential through the variable T*. HWCA(Qσeff , φeff , T *) ≥ HHS(Qσeff , φeff ) HWCA(Qσeff , φeff , T *) → HHS(Qσeff , φeff ) for T * → 0, φeff → 0 HWCA(Qσeff , φeff , T *) = HHS(Qσeff , φeff ) Figure 13. ΔωWCA(Q) for T* = 2/3 and different densities (points) compared with ΔωHS(Q) calculated according to eq 16 using the σeff determined by eq 18 (lines). The arrows indicate Qmax positions.

at Q ≈ Q max

(20)

The microscopic dynamics of the simulated WCA fluid shows significant difference in comparison to the HS system due to the fact that, in contrast to particles with a continuous pair potential, the mutual motions of hard spheres are uncoordinated between collisions. Nevertheless, if diffusion in a liquid is dominated by the repulsive part of the pair potential, the diffusive mobility of atoms depends entirely on their ability to push apart the neighbors that can be characterized using the effective diameter σeff. As has just been demonstrated, this parameter can be found from Δω(Q) using the HS approximation at Q ≈ Qmax, corresponding to the average distance between the nearest neighbors. In other words, the probability for a WCA particle to escape its equilibrium position is the same as for a hard sphere with σHS = σeff, and this results in equality of diffusion constants for both fluids. 3.3. Diffusion Coefficient of Liquid Rb. Further work is needed to find out whether the result obtained for the WCA fluid can be generalized for more complicated cases. In this section, it is demonstrated that, like in the case of a soft-sphere fluid, the HS approach works at least for simple liquid metals. Of course, a comparison of a model with arbitrary selected experimental data cannot be very convincing. In order to explore the problem systematically, the HS analysis of MD simulations with more realistic pair potentials should be performed. It should be noted that eq 4 can sometimes describe Δω(Q) measured for a real liquid with a reasonable accuracy using very realistic σeff and DE parameters (see for example Figure 8 in ref 25). The reason for this is that, at moderate densities, eq 4 provides a good agreement with ΔωHS(Q) at Qmax (Figure 5). At the same time, the strong overestimation of ΔωHS(Q) by eq 4 at the right side of Qmax mimics the behavior of Δω(Q) typical for a soft-core potential. Even if the applicability of the HS model is restricted to the narrow Q-range around the minimum of Δω(Q), this approach is still very attractive due to its simplicity. If particle density n for given T is known, σeff is the only unknown quantity that has to be found from Δωexp(Q). Indeed, supposing the validity of the HS approach, one should not consider σeff, DHS, and SHS(Qmax) as independent parameters. As an example, one can try to find the diffusion constants for liquid Rb using the relaxation times of I(Q,t) measured for this metal at Qmax and published in ref 31. All required parameters, except the atomic mass of Rb, are listed in Table 1. Table 1 also

Figure 14. ΔωWCA(Q) for nσ3 = 0.917 and different temperatures (points) compared with ΔωHS(Q) calculated in the same way as for Figure 13 (lines).

Figure 15. Simulations results for nσ3 = 0.917 at T* = 2.35 and T* = 11.3 compared with the data obtained at T* = 2/3 for nσ3 = 0.765 (left) and for nσ3 = 0.573 (right).

are represented in the form of Δω WCA /(Q 2 T * σeff ) as a function of Qσeff in Figure 15. 9975

DOI: 10.1021/acs.jpcb.6b04685 J. Phys. Chem. B 2016, 120, 9969−9977

Article

The Journal of Physical Chemistry B

diameter, σeff, and this parameter can be found from Δω(Q) at Q ≈ Qmax. The similar result has been obtained for the neutron scattering data measured for liquid Rb.

Table 1. Parameters Used for the Evaluation of σeff for Liquid Rb and the Results of the Analysis. T (K)

n (Å−3)

Qexp (Å−1)

320

0.01034

1.52

330

0.01031

370

0.01018

400

0.01007

450

0.00989

Δωexp (meV) 0.315 ± 0.005 0.347 ± 0.005 0.430 ± 0.006 0.484 ± 0.008 0.598 ± 0.014

σeff (Å)

D (10−5 cm2/s)

4.412 ± 0.008

3.01 ± 0.11

4.392 ± 0.008

3.42 ± 0.12

4.357 ± 0.01

4.51 ± 0.17

4.338 ± 0.015

5.32 ± 0.27

4.29 ± 0.02

7.13 ± 0.4



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Telephone: +49 89 289 14754. Notes

The author declares no competing financial interest.



ACKNOWLEDGMENTS I am grateful to Prof. A.G. Novikov for fruitful discussions. I also thank the organizers and lecturers of MDANSE school 2014 for the excellent training in running and analyzing MD simulations.

shows the resulting σeff and D. The atomic density n is interpolation of the values given in ref 49. Qexp is taken from ref 31. The experimental linewidths Δωexp are calculated according to eqs 12 and 13 using the values for τw and β given in ref 31. The parameters σeff were found by numerically solving the equation Δωexp(Qexp) = ΔωHS(Qexp), where ΔωHS(Q) was calculated according to eq 16. Figure 16 demonstrates that the resulting diffusion constants D found using eq 10 are very close to the values measured by the isotope electrotransport method.50



REFERENCES

(1) Pilgrim, W.-C.; Morkel, C. State dependent particle dynamics in liquid alkali metals. J. Phys.: Condens. Matter 2006, 18, R585−R633. (2) Meyer, A. Self-diffusion in liquid copper as seen by quasielastic neutron scattering. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 012102. (3) Meyer, A. Atomic transport in dense multicomponent metallic liquids. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 66, 134205. (4) Meyer, A.; Horbach, J.; Heinen, O.; Holland-Moritz, D.; Unruh, T. Self Diffusion in Liquid Titanium: Quasielastic Neutron Scattering and Molecular Dynamics Simulation. Defect Diffus. Forum 2009, 289292, 609. (5) Mavila Chathoth, S.; Meyer, A.; Koza, M. M.; Juranyi, F. Atomic diffusion in liquid Ni, NiP, PdNiP, and PdNiCuP alloys. Appl. Phys. Lett. 2004, 85, 4881−4883. (6) Meyer, A.; Stüber, S.; Holland-Moritz, D.; Heinen, O.; Unruh, T. Determination of self-diffusion coefficients by quasielastic neutron scattering measurements of levitated Ni droplets. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 092201. (7) de Gennes, P. G. Liquid Dynamics and Inelastic Scattering of Neutrons. Physica 1959, 25, 825−839. (8) de Schepper, I. M.; Cohen, E. G. D. Collective modes in fluids and neutron scattering. Phys. Rev. A: At., Mol., Opt. Phys. 1980, 22, 287−289. (9) de Schepper, I. M.; Cohen, E. G. D.; Zuilhof, M. J. The width of neutron spectra and the heat mode of fluids. Phys. Lett. A 1984, 101, 399−404. (10) Cohen, E. G. D.; de Schepper, I. M.; Zuilhof, M. J. Kinetic theory of the eigenmodes of classical fluids and neutron scattering. Physica B+C (Amsterdam) 1984, 127, 282−291. (11) Kamgar-Parsi, B.; Cohen, E. G. D.; de Schepper, I. M. Dynamical processes in hard-sphere fluids. Phys. Rev. A: At., Mol., Opt. Phys. 1987, 35, 4781−4795. (12) Cohen, E. G. D.; Westerhuijs, P.; de Schepper, I. M. Half Width of Neutron Spectra. Phys. Rev. Lett. 1987, 59, 2872−2874. (13) Montfrooij, W.; de Schepper, I. Excitations in simple liquids, liquid metals and superfluids; Oxford University Press: Oxford, 2010. (14) Erpenbeck, J. J.; Wood, W. W. Self-diffusion coefficient for the hard-sphere fluid. Phys. Rev. A: At., Mol., Opt. Phys. 1991, 43, 4254− 4261. (15) Shimoji, M. Liquid Metals; An Introduction to the Physics and Chemistry of Metals in the Liquid State; Academic Press: London, 1977. (16) Protopapas, P.; Andersen, H. C.; Parlee, N. A. D. Theory of transport in liquid metals. I. Calculation of selfdiffusion coefficients. J. Chem. Phys. 1973, 59, 15−25. (17) van Loef, J. J. Hard core size at the melting point and Lindemann law. J. Chem. Phys. 1974, 61, 1605. (18) van Loef, J. J. Atomic and electronic transport properties and the molar volume of monoatomic liquids. Physica 1974, 75, 115−132.

Figure 16. Diffusion constants D for liquid Rb obtained from the neutron scattering data31 using the HS approximation (circles) compared with the experimental values measured by the isotope electrotransport method, represented as 0.66 × 10−3 exp(1.98(kcal/ mol)/RT)50 (line).

4. CONCLUSION The test performed by means of HS molecular dynamic simulations reveals a significant discrepancy between the simulation results and eq 4. The main output of the analysis is the finding that the ΔωHS(Q) behavior is defined in terms of the average cage size ⟨Lc⟩, rather than of σHS. The empirical improvement of the model is proposed, which describes ΔωHS(Q) with good accuracy in the Q-range corresponding to the first maximum of S(Q) and for the packing fractions φ = 0.43−0.49. One cannot expect that the HS approach can be applied to real liquids without any limitations. In order to explore these limitations, as the first step, the simulated ΔωHS(Q) has been compared with the results for the soft-spherical potential. The microscopic dynamics of the simulated soft-sphere fluid shows significant difference in comparison to the HS system. Nevertheless, the diffusive mobility of soft spheres can be characterized within the HS approximation using an effective 9976

DOI: 10.1021/acs.jpcb.6b04685 J. Phys. Chem. B 2016, 120, 9969−9977

Article

The Journal of Physical Chemistry B (19) Demmel, F.; Szubrin, D.; Pilgrim, W.-C.; Morkel, C. Diffusion in liquid aluminium probed by quasielastic neutron scattering. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 014307. (20) Scopigno, T.; Di Leonardo, R.; Comez, L.; Baron, A. Q. R.; Fioretto, D.; Ruocco, G. Hard-Sphere-like Dynamics in a Non-HardSphere Liquid. Phys. Rev. Lett. 2005, 94, 155301. (21) Bermejo, F. J.; Bustinduy, I.; Cabrillo, C.; Levett, S. J.; Taylor, J. W. Comment on ‘‘Hard-Sphere-Like Dynamics in a Non-Hard-Sphere Liquid’’. Phys. Rev. Lett. 2005, 95, 269601. (22) Scopigno, T.; Di Leonardo, R.; Comez, L.; Baron, A. Q. R.; Fioretto, D.; Ruocco, G.; Montfrooij, W. Reply. Phys. Rev. Lett. 2005, 95, 269602. (23) Bermejo, F. J.; Bustinduy, I.; Levett, S. J.; Taylor, J. W.; Fernández-Perea, R.; Cabrillo, C. Experimental evidence for the presence of nonacoustic excitations in molten Ga. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 104103. (24) Blagoveshchenskii, N.; Novikov, A.; Puchkov, A.; Savostin, V.; Sobolev, O. Self-diffusion in liquid gallium and hard sphere model. EPJ Web Conf. 2015, 83, 02018. (25) Bodensteiner, T.; Morkel, Chr.; Gläser, W.; Dorner, B. Collective dynamics in liquid cesium near the melting point. Phys. Rev. A: At., Mol., Opt. Phys. 1992, 45, 5709−5720. (26) Ruiz-Martín, M. D.; Jiménez-Ruiz, M.; Plazanet, M.; Bermejo, F. J.; Fernández-Perea, R.; Cabrillo, C. Microscopic dynamics in molten Ni: Experimental scrutiny of embedded-atom-potential simulations. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 224202. (27) Demmel, F.; Diepold, A.; Aschauer, H.; Morkel, C. Temperature dependence of the de Gennes narrowing in liquid rubidium. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 104207. (28) Badyal, Y. S.; Bafile, U.; Miyazaki, K.; de Schepper, I. M.; Montfrooij, W. Cage diffusion in liquid mercury. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 68, 061208. (29) Ruiz-Martín, M. D.; Jiménez-Ruiz, M.; Stunault, A.; Bermejo, F. J.; Fernández-Perea, R.; Cabrillo, C. Stochastic dynamics in molten potassium explored by polarized quasielastic neutron scattering. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 174201. (30) Blagoveshchenskii, N. M.; Novikov, A. G.; Savostin, V. V. Selfdiffusion in liquid lithium and lead from coherent quasi-elastic neutron scattering data. Phys. Solid State 2014, 56, 120−124. (31) Demmel, F.; Morkel, C. Nonexponential relaxation in a simple liquid metal. Phys. Rev. E 2012, 85, 051204. (32) Rapaport, D. C. The art of molecular dynamics simulation; Cambridge University Press: Cambridge, 1995. (33) Alley, W. E.; Alder, B. J.; Yip, S. The neutron scattering function for hard spheres. Phys. Rev. A: At., Mol., Opt. Phys. 1983, 27, 3174− 3186. (34) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of Repulsive Forces in Determining the Equilibrium Structure of Simple Liquids. J. Chem. Phys. 1971, 54, 5237−5247. (35) Smith, W.; Forester, T. R. DL_POLY_2.0: a general-purpose parallel molecular dynamics simulation package. J. Mol. Graphics 1996, 14, 136−141. (36) Róg, T.; Murzyn, K.; Hinsen, K.; Kneller, G. R. nMoldyn: a program package for a neutron scattering oriented analysis of molecular dynamics simulations. J. Comput. Chem. 2003, 24, 657−667. (37) Wertheim, M. S. Exact Solution of the Perus- Yevick Integral Equation for Hard Spheres. Phys. Rev. Lett. 1963, 10, 321−323. (38) Thiele, E. Equation of State for Hard Spheres. J. Chem. Phys. 1963, 39, 474−479. (39) Verlet, L.; Weis, J. J. Equilibrium Theory of Simple Liquids. Phys. Rev. A: At., Mol., Opt. Phys. 1972, 5, 939−952. (40) Henderson, D.; Grundke, E. W. Direct correlation function: Hard sphere fluid. J. Chem. Phys. 1975, 63, 601−607. (41) Colot, J. L.; Baus, M.; Xu, H. The freezing of hard spheres. Mol. Phys. 1986, 57, 809−823. (42) Heyes, D. M. System size dependence of the transport coefficients and Stokes−Einstein relationship of hard sphere and Weeks−Chandler−Andersen fluids. J. Phys.: Condens. Matter 2007, 19, 376106.

(43) Alder, B. J.; Wainwright, T. E. Velocity Autocorrelations for Hard Spheres. Phys. Rev. Lett. 1967, 18, 988−990. (44) Alder, B. J.; Gass, D. M.; Wainwright, T. E. Studies in Molecular Dynamics. VIII. The Transport Coefficients for a Hard-Sphere Fluid. J. Chem. Phys. 1970, 53, 3813−3826. (45) Alder, B. J.; Wainwright, T. E. Decay of the Velocity Autocorrelation Function. Phys. Rev. A: At., Mol., Opt. Phys. 1970, 1, 18−21. (46) Balucani, U.; Zoppi, M. Dynamics of the Liquid State; Clarendon Press: Oxford, 1994. (47) Morkel, C.; Gläser, W. Single-particle motion in liquid sodium. Phys. Rev. A: At., Mol., Opt. Phys. 1986, 33, 3383−3390. (48) Heyes, D. M.; Okumura, H. Equation of state and structural properties of the Weeks-Chandler-Andersen fluid. J. Chem. Phys. 2006, 124, 164507. (49) Waseda, Y. The Structure of Non-Crystalline Materials; McGrawHill: New York, 1980. (50) Larsson, S. J.; Roxbergh, C.; Lodding, A. Self-Diffusion in liquid alkali metals. Phys. Chem. Liq. 1972, 3, 137.

9977

DOI: 10.1021/acs.jpcb.6b04685 J. Phys. Chem. B 2016, 120, 9969−9977