Viscosity estimates of liquid metals and warm dense matter using the Yukawa reference system

Available online at www.sciencedirect.com

High Energy Density Physics 4 (2008) 49e57 www.elsevier.com/locate/hedp

Viscosity estimates of liquid metals and warm dense matter using the Yukawa reference system Michael S. Murillo* Physics Division, MS D410, Los Alamos National Laboratory, Los Alamos, NM 87545, USA Received 26 September 2007; received in revised form 13 November 2007; accepted 28 November 2007 Available online 14 December 2007

Abstract A practical method of computing the viscosity of liquid metals and warm dense matter over wide ranges in parameters is proposed. The method is based on mapping the system of interest onto the Yukawa model, for which the viscosity is well known and can be written in a quasiuniversal form. Comparisons are made with quantum molecular dynamics results for compressed iron relevant to the Earth’s core, experimental data for many liquid metals, and simulation results for dense deuterium relevant to inertial confinement fusion experiments. Finally, the dispersion and damping of ion-acoustic waves in warm dense matter are considered in this context. Ó 2007 Elsevier B.V. All rights reserved. PACS: 52.27.Gr; 52.25.Fi; 91.35.x Keywords: Viscosity; Dense plasma; Warm dense matter; Yukawa; Yukawa system; Yukawa model; Reynolds number; Ion-acoustic wave; ICF; Earth’s core; Wave damping; Ionization level

1. Introduction The hydrodynamic behavior of materials is determined in part by the material viscosity. There are several physical systems for which theoretical predictions of the viscosity are needed over a wide parameter range due to the inaccessibility of experimental methods in the regime of interest. For example, liquid iron under high pressure is difficult to achieve in the laboratory, but its viscosity is important for modeling planetary bodies [1] such as the Earth [2,3]. Turbulent mixing in inertial confinement fusion (ICF) targets also requires knowledge of the viscosity [6] over a large range of parameter space; viscosity controls the growth of instabilities that arise from imperfections in the radiation source and the target, as well as those at the ablation surface [7]. Knowledge of the viscosity is also important for understanding wave damping [8] in dense plasmas, which may play a role in X-ray Thomson scattering

* Tel.: þ1 505 667 6767; fax: þ1 505 667 7684. E-mail address: [email protected] 1574-1818/$ - see front matter Ó 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.hedp.2007.11.001

experiments in warm dense matter (WDM) [9] that, in the future, may attempt to probe ion-acoustic waves. Although various methods exist to compute the viscosity, results for liquid metals (LMs) can vary over many orders of magnitude [2]. Similarly, for dense ICF plasmas, different models [6] vary over orders of magnitude. In principle, a very accurate and straightforward method to obtain transport coefficients is to use molecular dynamics (MD) simulation to evaluate the relevant KuboeGreen [10] expression. For the cool to warm materials of interest here, this requires a somewhat sophisticated treatment of the electrons, which are degenerate for liquid metals and partially degenerate for WDM. Moreover, for moderately heavy elements, atomic physics plays a role in determining the ionization level. In this context, quantum MD (QMD) is a promising approach, and has been used to compute the viscosity for LM [2,3] and WDM [11]. However, QMD remains too computationally expensive to provide data over a wide range. This is particularly true for WDM, since QMD is only efficient on the cool side below about T w 10 eV, even for simple quantities such as the energy or pair correlation function. Because of this, in

M.S. Murillo / High Energy Density Physics 4 (2008) 49e57

50

some previous works, only the diffusion coefficient was obtainable, and the StokeseEinstein relation was relied upon to get the viscosity [2]. In cases where the viscosity could be obtained directly [3,11], limited results were computed with significant error bars. An alternative approach is to map the system of interest onto a similar system for which the viscosity is well known [12]. For liquid iron, a simple reference system has been used [3] in which a power-law pair potential was found to yield viscosities in rough agreement with the QMD results. Such approaches allow for more data to be obtained, convergence issues are more readily addressed, and scaling studies are possible. For ICF plasmas, the dense plasma has been mapped onto the one-component plasma (OCP) model [6], and the known MD results of that model were used. Here, we follow such approaches by mapping LM and WDM systems onto the Yukawa model, which is the screened version of the OCP, and employ accurate viscosities from recent MD studies. Our focus is on the ionic viscosity, since that dominates the overall fluid viscosity. 2. The Yukawa model LM and WDM are both strongly coupled Coulomb systems, as characterized by the Coulomb coupling coefficient: G¼

Q2 ; aT

ð1Þ

where Q is the charge, a ¼ (3/4pn)1/3 is the ion-sphere radius in terms of the ionic density n, and T is the temperature in energy units. Moreover, because both systems are partially ionized, atomic physics plays a role in determining the charge state of the ionic core and the strength of the effective ioneion interaction. The difference between the two systems is that WDM has thermal core excitation and is only partially degenerate. In other words, in terms of the degeneracy parameter: T q¼ ; EF

ð2Þ

where EF is the usual Fermi energy, LM have q z 0 whereas WDM has q w 1. Relative to the OCP model, LM and WDM are effectively less weakly coupled due to screening by free electrons. Recently, it was shown [11] that a simple viscosity model [12] was able to predict QMD results for the viscosity of dense hydrogen with reasonable accuracy. Since the simple model used the Yukawa system, this suggests that mapping the systems of interest onto the Yukawa system has some merit. The Yukawa model is defined in terms of the pair potential: uðrÞ ¼

G kr e ; r

ð3Þ

where the potential v(r) ¼ Tu(r) is measured in units of temperature and distances are measured in terms of a. In these

units, the static properties of the Yukawa model (YM) are completely characterized by the two parameters G and k. The viscosity of the YM has been computed by MD to reasonable accuracy [14]. The results can be put into a very simple, quasiuniversal form if appropriate dimensionless units are used. We write the viscosity in terms of the characteristic viscosity: pffiffiffi ð4Þ h0 ¼ 3uE Mna2 ; where uE ¼ uE(k) is the Einstein frequency and M is the mass of an ion. The motivation for these units is given in Appendix A. In these units, the viscosity is given by h Gm G ¼ 0:0051 þ 0:374 þ 0:022; h0 Gm G

ð5Þ

where the ratio G/Gm characterizes the state of the system; here, Gm is the melt boundary in the Gek plane. The coefficients were obtained by a fit of the MD data [14] for the cases k ¼ 2 and k ¼ 3, taken together, to place an average emphasis on k values in this range. Although this intentional bias has been used, we can expect that the form (5) still has wide applicability as a result of the quasiuniversal form that is used. Note that this scaling is different from, but similar in spirit to, other scaling arguments [3]. The advantages of using (5) are that the StokeseEinstein relation need not be used, the quasiuniversal form ensures that the result can be used over a very wide range of parameter space, accurate results are already known, generalizations (to mixtures, for example) can be easily obtained, and the model is physically motivated by the physics of screened Coulomb systems. In practice, it is useful to have simple fits to some of the quantities needed to evaluate (5). The melt boundary [13] is located approximately at   1:38 Gm ðkÞ ¼ 171:8 þ 82:8 e0:565k  1 ; ð6Þ as obtained from a fit to the MD data. Similarly, the Einstein frequency can be written approximately in terms of the ionic plasma frequency as pffiffiffi 1:62 3uE ðkÞ ¼ ui e0:2k ; ð7Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi from a fit to MD data [14], where ui ¼ 4pnQ2 =M is the ion plasma frequency. 3. Mapping liquid metals and warm dense matter Here, we propose to use the YM for both LM and WDM. It remains to be determined how to map the systems onto the YM so that the results of the previous section can used; unfortunately, there is no unique way to do this. One could, for example, ensure that the pair correlation functions are similar, which would require a QMD simulation for the LM or WDM, although it would be less intensive than directly obtaining the viscosity. Similarly, one could use energy fluctuations [3] to determine the optimal YM parameters.

M.S. Murillo / High Energy Density Physics 4 (2008) 49e57

Here, we take a physically based approach in which the parameters G and k are treated not as fitting parameters, but as physical quantities. As mentioned above, LM and WDM are characterized by partial ionization. We therefore assume that the ionization state Z can be computed and define Q ¼ Ze; this assumes that the polarization of the atomic core can be neglected. We take the ion-sphere radius to be related to the mass density in an obvious way through r ¼ Mn. The temperature T is the ion temperature Ti, which need not be equal to the electron temperature Te. Specifically, this leads to G¼

½Zðr; Te ; Z; AÞe Ti

2

4pr 3M

1=3 :

ð8Þ

Here, A is the atomic number and Z is the nuclear charge. In practice, we will obtain Z from a simple ThomaseFermi (TF) model, for which a good fit is already available [15]. This TF model has a self-consistent finite-density and finitetemperature treatment of all electronic states. Next, we need to compute k. When there are bound electrons, the interaction between the valence and the core electrons is weak [16] and we can expect that a linear screening model can yield k. Again, the TF model provides a good starting point. The simple interpolation form, vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u 4pZne2 k ¼ atqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð9Þ 2ffi; 2 2 Te þ ð3EF retains the very low and very high temperature limits. The Fermi energy is given by 2=3  Z2 3p2 Zr EF ¼ ; M 2m

ð10Þ

where m is the electron mass. Given the mass density r, the ion mass M, the atomic number A, the nuclear charge Z, and the temperatures Te and Ti, all parameters are completely determined and the viscosity is given by Eq. (5). Taken together, we will refer to Eqs. (4)e (10) as the ‘‘Yukawa viscosity model’’ (YVM). 4. Applications Given the basic YVM result, we now turn to several applications. These applications were chosen for their intrinsic interest and, in most cases, because other (experimental or computational) methods have yielded results. This section starts with liquid metals, which represent the most challenging application of the YVM because the temperatures are so low that the YM is suspect. Compressed iron relevant to the Earth’s core is considered first, followed by a wide range of liquid metals under normal conditions. Then, dense hydrogen is considered, which should also represent a challenge for the YVM because TF theory is not accurate for very low-Z elements; but, QMD results are available to compare with.

51

Finally, the YVM is applied to ion-acoustic waves in WDM, and it is shown that the viscosity plays a central role in both the dispersion and the damping of these waves. 4.1. Liquid iron: the Earth’s core Below the Earth’s thin silicate crust and the mantle layers lies a liquid iron-dominated alloy heated by radioactive decay. Processes such as geomagnetism are governed by hydrodynamic motions of this liquid metal. The inner core is believed to be highly compressed solid iron [17]. It is therefore of great interest to understand the viscosity of compressed iron [2]. Initial computations [2] on the viscosity of liquid iron (l-Fe) relied upon the StokeseEinstein relation: h¼

T ; 2pRD

ð11Þ

where D is the diffusion coefficient and R is the ‘‘diameter’’ of a particle. Since D is obtained from the velocity autocorrelation function, a single particle property, it is more easily obtained than the stress autocorrelation function needed to get h directly. There are many known issues with Eq. (11) (mixing particle and fluid descriptions, boundary conditions, etc.), and such results may only be correct to about a factor of 3 [2]. Subsequent QMD simulations [3] were able to evaluate h directly, albeit only for one representative point, and a reference system was introduced that employed the pair potential f(r) ¼ B/ra to investigate size effects. A priori, it is not clear that the YVM can model l-Fe well. In particular, compressed iron has a complex phase diagram that is determined by complex electronic structure [17,3]. The QMD case considered previously [3] was for r ¼ 10.7 g/cm3 and T ¼ 4300 K, which is representative of the outer core. Their result is hQMD ¼ 8.5 mPa s, whereas the YVM yields hYVM ¼ 3.2 mPa s. This suggests that the YVM is somewhat predictive. Moreover, the YVM can be used to obtain results over the Earth’s entire range of densities r w 4e12 g/cm3 and temperatures T w 4000e5000 K with negligible computational effort. Alternately, a few more QMD points could examine the quasiuniversal scaling behavior given here, or elsewhere [3], which would limit the number of simulations that are necessary to cover the entire mantle and core regions. 4.2. Liquid metals near melt Given the success, relative to QMD, of the YVM for compressed l-Fe, we can quantify how predictive the YVM is for liquid metals for which experimental data exist [18]. In Fig. 1 the viscosity of several liquid metals at the melting point is shown; clearly, the YVM captures the general trend over the range from Al to Pu. Interestingly, the YVM value for the Pu viscosity is closer to the data [4] than modified embedded atom method (MEAM) MD results [5]. Overall, the YVM predicts viscosities that are too low, but qualitative trends are captured.

52

M.S. Murillo / High Energy Density Physics 4 (2008) 49e57

the agreement tends to improve at higher temperatures, where we would expect the Yukawa model to be more appropriate; note, however, that the temperatures in Fig. 2 are still quite low (below vaporization). 4.3. Hydrodynamic instabilities in inertial confinement fusion

Fig. 1. The viscosity of various elements (Al, K, Fe, Cu, Ag, Yb, Au, Pb, U, and Pu) at their melting temperatures is shown. Experimental data are shown in red and predictions of the YVM model are shown in blue. (For interpretation of colors in this article, the reader is referred to the web version of the article.)

Because of the interest in higher temperature systems, it is also of interest to compare experimental data to the YVM at elevated temperatures; this is shown in Fig. 2 for Au and Ag. Again, the general trend is captured with the YVM. Note that

Fig. 2. The viscosities of l-Ag and l-Au are shown at elevated temperatures. Red circles are experimental data, whereas the blue lines reflect the YVM. Note that the YVM greatly improves with increasing temperature. (For interpretation of colors in this article, the reader is referred to the web version of the article.)

Implosion symmetry is important in ICF to minimize mix [19]. Asymmetry is caused by the growth of various instabilities, including RayleigheTaylor, RichtmeyereMeshkov, and KelvineHelmholtz [7]. Such instabilities are characterized by the relative size of the inertial and viscous terms in the NaviereStokes equation, as measured by the Reynolds number. To quantify these hydrodynamic flows, an OCP viscosity model similar to the one presented here has been presented [6] and Reynolds number has been computed for various ICF scenarios. The primary advantage of the YVM is that partial ionization is accounted for via Z and electron screening is included via k. There are two important material regions in ICF implosions: the hydrogen fuel and the pusher, which itself may be hydrogen or some other moderate-Z material (e.g., SiO2). It is mixing at the fuelepusher interface that can spoil the burn conditions of the fuel. Very little is currently known about the viscosity of various pushers, which explore a wide regime of parameter space [6]. We therefore focus on the viscosity of the fuel, for which QMD calculations of the viscosity have been obtained [11]. These QMD data are especially useful for validating the YVM, since they span the molecular, atomic, and dense plasma regimes. Results are shown in Table 1 for the three densities r ¼ 0.665, 1.0, 1.5 g/cm3 over a wide range of temperatures. It is seen that the YVM agrees very well with the QMD data except when the deuterium is cool and expanded. The results are also shown in Fig. 3. Note that the QMD has considerable fluctuations relative to the YVM, and that both methods agree very well above about T w 104 K. The largest fluctuations in the QMD, where the largest disagreement occurs, are at low temperature. The disagreement is expected, of course, since the YVM is too simple to treat neutral atomic interactions and the formation of molecules. It has been suggested [11] that the molecules dissociate to a plasma state when T > 104 K; the results of the present model tend to be accurate above about half that temperature. In general, the present model is within the error bars of the QMD data, with essentially no computational cost. In the parameter regime treated, the ionization state was nearly constant around Zw0:58  0:67, because the density was not varied over a very large range (recall that a w n1/3) and q is small for all temperatures. Note also that k is fairly large, which obviates the use of an OCP model [6] in this regime. At higher temperatures, we expect that the limits Z/1 and k/0 will be obtained and OCP results are applicable. At very high temperatures, where thermonuclear burn is to be achieved, MD simulations

M.S. Murillo / High Energy Density Physics 4 (2008) 49e57

53

Table 1 Viscosity of dense deuterium as compared with detailed MD simulations [11]

4.4. Ion-acoustic wave damping in warm dense matter

r (g/cm3)

T (K)

h (MD) (mPa s)

h (YVM) (mPa s)

G

k

0.665

792 1350 3560 4780 5878 11 020 19 800 30 250 50 100

0.28 0.49 0.24 0.23 0.21 0.16 1.18 1.5 1.33

0.11 0.11 0.18 0.23 0.27 0.48 0.83 1.22 1.9

68 40 15 11.4 9.3 5.0 2.8 1.9 1.2

2 2 2 2 2 2 2 1.9 1.8

1.0

896 2019 3091 5700 7760 11 300 19 800 29 677 49 700

0.61 0.71 0.44 0.2 0.37 0.65 1.4 1.23 1.56

0.19 0.18 0.21 0.32 0.41 0.56 0.94 1.4 2.2

79 35 23 12.5 9.2 6.3 3.6 2.5 1.5

1.9 1.9 1.9 1.9 1.9 1.9 1.88 1.85 1.8

WDM is matter with G w 1 and q w 1, which tends to correspond to conditions just below solid density with a temperature in the few tens of eV range. Theoretical modeling in this regime is particularly difficult because of simultaneous atomic core excitation and ionization, electronic and ionic correlations, and partial electron degeneracy. Because of the energy density needed to form WDM, experiments are equally difficult, and diagnostics are largely unavailable. The advent of X-ray Thomson scattering (XRTS) promises to change this situation [9] by providing an external probe that can reveal most of the interesting properties of WDM. To date, most of the focus on XRTS has been toward features of the electron plasma wave (EPW). The reasons for this can be understood by noting that the cross-section is given by [20]

1.5

7700 11 800 19 900 29 422 49 700

0.35 0.64 1.2 1.33 2.5

0.48 0.68 1.09 1.56 2.55

11.8 7.8 4.6 3.2 1.9

1.8 1.8 1.8 1.8 1.7

are unreliable, and standard Braginskii results become useful [6]. Given the success of the present model for the cases considered so far, it is reasonable to assume that it will yield useful predictions of the viscosity of the pusher as well, which will tend to be denser, cooler, and possibly of higher Z than the fuel. For some materials, such as SiO2, the present model will need to be generalized to mixtures.

Fig. 3. The data from Table 1, for dense deuterium, are shown. The MD (red) and YVM (blue) viscosity predictions are shown versus temperature for the three densities r ¼ 0.665 (solid), r ¼ 1.0 g/cm3 (long dash), and r ¼ 1.0 g/cm3 (dash-dot). (For interpretation of colors in this article, the reader is referred to the web version of the article.)

d2 s ¼ sT Stot ee ðq; uÞ; dU du

ð12Þ

where sT ¼ ðe4 =m2e c4 Þ ðbe1 $be2 Þ2 ðk1 =k2 Þ is in terms of the polarizations be and wavevectors k of the incident and scattered waves. The important quantity here is the electroneelectron dynamic structure factor Stot ee (q, u), taken here for all of the electrons in the system, which can be written as a sum over terms that account for ion motion, pure free-electron motion (i.e., independent of the ions), and core electron motion [20]. The pure free-electron term gives rise to peak in the cross-section nearpthe EPW, ffior, equivalently, near the freffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi quency uEPW ¼ 4pne e2 =me ; clearly, this peak reveals the free-electron density. Temperature can be obtained from the asymmetry of the peaks at uEPW [22]. The simplicity of this situation is the reason that XRTS could be an important WDM diagnostic, especially in the context of next generation light sources. There are details that are of interest, such as the electron dynamic structure factor [21], but stringent tests of WDM physics are likely to occur in the ionic component of the XRTS signal. As mentioned above, WDM is largely of interest because of atomic physics issues, which are the result of electrone ion coupling. For moderate-Z elements, the ions to which the electrons are coupled to have strong, liquid-like correlations. Thus, the physics contained in the ionic terms of XRTS is much richer in physics. Moreover, since the basic collective mode of the ions is the ion-acoustic wave (IAW), issues arise from the fact that a wide range of frequencies is explored over the dispersion of the waves [8]. At long wavelengths the modes are of very low frequency and reflect hydrodynamic behavior, whereas, at short wavelengths, the modes are high frequency elastic-like modes that reflect the underlying quasilattice formed by strong ioneion coupling. In simple models [8], the time scale separation is set by the viscoelastic relaxation time, which reveals that the viscosity of WDM plays a central role. In some systems, this simple picture is not sufficient [23]. Experimentally, WDM IAWs will be difficult to measure because of their low frequency,

54

M.S. Murillo / High Energy Density Physics 4 (2008) 49e57

although developments on next generation light sources will help mitigate this. The IAW is best described by the ioneion dynamic structure factor Sii(q, u), which enters into the XRTS cross-section through the electroneion coupling terms [20]. We can construct a model for this structure factor using a modified NaviereStokes equation, which leads to pffiffiffi ~ 8 3 uE q4 h ui Sii ðq; uÞ ¼  2 ;  2  pffiffiffi 2 9Gui q 4 3u E 2 2 u  þ q u~ h 3GSðqÞ 3ui ð13Þ 1 ~¼ where pffiffiffi u is in2units of ui, wavevectors are in units of a , h h= 3uE Mna , and S(q) is the static structure factor. (Here, we neglect the bulk viscosity.) The dispersion of the IAW is then sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi  pffiffiffi 2 q2 4 3 uE 2 ~ :  ð14Þ qh uðqÞ ¼ 3GSðqÞ 3ui

Note how the static structure and the viscosity enter directly into the dispersion relation. Moreover, is clear from Eq. (14) that damping of the IAW is controlled by the viscosity. Despite this richness of physics, the form (14) is not complete, since it does not include, for example, Landau damping and does not satisfy higher-frequency sum rules. Thus, experimental measurement of the IAW will likely provide many challenges to WDM theory. We now turn to the specific example of warm dense Be at half solid density, r ¼ 0.925 g/cm3, and with a temperature of T ¼ 20 eV. We will treat the system as a Yukawa system and assume that Eq. (14) is accurate. For these conditions, we have Z ¼ 2:0, G ¼ 1.86, k ¼ 1.6, and h ¼ 3.06 mPa s. Dispersion relations are shown in Fig. 4 according to three different approximations. The top (orange) curve is the simplest form of the theory in which viscosity is neglected and the long wavelength isothermal compressibility (S(0) ¼ nTcT) is taken to be a constant. (The structure factor is here computed in the hypernetted chain approximation.) If nonlocal compressibility cT(q) is included, one obtains the magenta curve, which reflects the ionic ordering. Finally, the full result (14) is shown as a cyan curve, which reveals how the wave frequency is downshifted due to viscous damping, and a negative dispersion regime appears above about q w 0.4. Because Eq. (14) is based on the NaviereStokes equation, it becomes less accurate at larger q, so only a small range (q ¼ 0e0.5) is shown. In the broader context of XRTS in WDM, it should be noted that the result of Fig. 5 will be modified further by q-dependent atomic form factors [20].

Fig. 4. Three dispersion relations are shown for the IAW in warm dense Be. The orange curve neglects viscosity and nonlocal compressibility. The magenta curve again neglects viscosity but includes the full S(q). Finally, the cyan curve is the full result. (For interpretation of colors in this article, the reader is referred to the web version of the article.)

more computationally challenging because some of the elements may exist in trace amounts that cannot be represented well in simulations with a few hundred particles. In this section, a simple method for extending the YVM to mixtures is discussed. Consider a mixture of elements with concentrations xi ¼ ni/n, where ni is the density of species i and n is the total nuclear density. For such a mixture, we can define two average charge states. Each species will have an ionization level denoted by Zi , whereas the average ionization of the material is given by P Zg ¼ Z gi xi , with g ¼ 1. The species ionization can be obtained i iteratively by ensuring that the electron density is continuous

5. Mixtures So far, we have been concerned with pure substances. In many situations, however, a mixture of various elements is present. For example, in the case of ICF, the Rayleighe Taylor instability can lead to mixing of heavier pusher elements into the lighter fuel regions. Such cases are even

Fig. 5. The full spectrum of Sii(q, u) is shown for the wavevectors q ¼ 0.1, 0.3, 0.5. Note the significant viscous damping of the IAW at larger q. (For interpretation of colors in this article, the reader is referred to the web version of the article.)

M.S. Murillo / High Energy Density Physics 4 (2008) 49e57

across ion-sphere cells of each species. Using the TF model discussed above, one varies the ionic density of each individual species (equivalently, the ion-sphere radii ai) until the electron density ne(ai) is the same at each ion-sphere boundary; then, one uses the relation ð4=3Þpne ðai Þa3i ¼ Zi to obtain the species ionization states. This simple procedure yields a self-consistent set of ionization states for each species under the partial pressures of the other species. To obtain the viscosity, we now use the species ionization states Z i g to construct an effective single-component fluid. Although there are several methods for reducing the mixture to a single component, we will adopt the ion-sphere mixing rule [25]. Using this rule, we add the energy contributions from each species, treated as independent ion spheres, to obtain the total energy (per particle in temperature units) P utot ¼ 0:9 Gi xi i

e2 P 2 a ¼ 0:9 xi Z i ; ai aT i 2 e 1=3 5=3 ¼ 0:9 Z Z aT

ð15Þ

where Gi ¼ Z2i e2 =ai T are the species coupling parameters and a is the ion-sphere radius for the total nuclear density. In the second line, the sum plays the role of an effective charge squared, which is defined in the third line. Note that there is an inconsistency between the uses of the ion-sphere radii ai in this step of the calculation and that above for the ionization state. In the ion-sphere mixing method, it is assumed that the electrons are uniform and, therefore, that ai in Eq. (15) refers to the radius an ion requires to be charge neutralized for a given Zi . In contrast, for the Zi calculation, the ion-sphere radii are chosen such that the nonuniform electron density is continuous at the cell boundaries. Non-linear mixing rules, including nonuniform screening, beyond this prescription are also available [26], although they do not have the simplicity of Eq. (15). Thus, a simple procedure to treat mixtures is to: (1) find the self-consistent ionization states Z i g of each species, (2) compute the effective coupling Geff ¼ Z1=3 Z 5=3 e2 =aT of the mixture, and P (3) use the effective coupling in the YVM, with M ¼ xi Mi in h0. This procedure has been tested for binary ionici mixtures (pure Coulomb with fixed ionization) [27,28] with the result that the predicted viscosity is somewhat lower than simulation results in regions with small fractions of the higher-Z element. 6. Conclusions In summary, a simple model for the viscosity of LM and WDM has been constructed based on a mapping of the system onto the Yukawa system, and using the known Yukawa viscosity. This model is referred to as the YVM. There are three primary motivations for constructing the YVM: (1) analytic approaches have yielded orders of magnitude differences in viscosity predictions, (2) molecular dynamics results that

55

include ionization (e.g., QMD) can only yield viscosities over a limited range of parameter space, which are very expensive to compute (relying in many cases on the approximate StokeseEinstein relation), and are limited to low temperatures (below 10 eV), and (3) many previous works did not include the effects of ionization and screening [6,27,28]. In contrast to these previous works, the YVM: (1) uses molecular dynamics results for the viscosity directly, (2) uses the TF ionization model (and its extension to mixtures) and (3) accounts for electron screening. As such, the YVM applies over an extremely wide range of parameters where other methods are very difficult or impossible. Although the YVM is intended for warmer regimes, such as WDM, where there are no reliable methods, comparisons with detailed MD studies and experimental data for LM suggest that the YVM works within a factor of a few even for systems just above the melting point. Further comparisons with QMD results for cool, dense hydrogen and LM at elevated temperatures suggest that the YVM becomes more accurate at higher temperatures. Therefore, the YVM appears to be of good accuracy for a very wide range of materials over large regions of temperature and density where there are no reliable methods, with essentially no computational cost. The results of the YVM have been applied to predictions of IAW damping in WDM. It is argued here that the IAW is much richer in physics than the EPW; therefore, despite the experimental difficulties in measuring properties of the IAW, much more stringent tests of theory can be made. Using a simple NaviereStokes model for the IAW, it has been shown that the IAW dispersion depends directly on the ionic structure and the viscosity, with the viscosity determining the damping of the wave. Much more work in this area is warranted; in particular, reconciliation of the NaviereStokes model with models based on use of the fluctuationedissipation theorem and sum rules will provide insights into viscoelastic phenomena [8,29] in WDM. MD can play an important role in calibrating our knowledge of IAW in WDM. The YVM itself can be improved and extended. First, the MD results [14] used here could be improved statistically, and with the use of many more particles. The simple form for the fit (5) could be trivially generalized. Issues associated with the computation of Z are discussed in Appendix B. The choice of the TF screening length could be improved; it has been suggested [11] that the TF model overestimates the strength of the screening. More work in mixtures is warranted; here a simple linear mixing model was suggested. The use of non-linear mixing rules, including screening, [26] and effective medium theory [28] may generally provide more accurate results, especially for extreme mass-charge asymmetries.

Acknowledgments The author is very grateful to Stephanie Hansen (LLNL) for supplying the Purgatorio Z results. This is article LA-UR-076565.

M.S. Murillo / High Energy Density Physics 4 (2008) 49e57

56

Appendix A. Quasiuniversal forms A feature of the YVM is the quasiuniversal form (5). It is interesting to note that the r.h.s. depends only upon structural quantities (G and Gm), whereas dynamical quantities (M, uE) appear only in the l.h.s. of the equation in h0. Thus, it is important to obtain a form for h0 such that the remaining terms are as quasiuniversal as possible; unfortunately, the choice is not unique [24]. Here, arguments are given for obtaining good choices of the reduced transport coefficient. For simplicity, the diffusion coefficient is discussed, although the method is more general. The diffusion coefficient is defined by [10] ZN D ¼ dt hvðtÞvð0Þi;

ðA1Þ

0

where v(t) is the velocity of a tagged particle in some particular direction. We wish to put this in a dimensionless form such that the final result can be written in terms of as few reduced variables as possible. There are three terms that need to be considered. First, there are the explicit variables t and v, and we can write 2

D ¼ ð~t~v Þ

ZN

dt hvðtÞvð0Þi=~v2 : ~t

h/i ¼

d Ge

1 Z

d GebH /;

ðA7Þ

which introduces the ion plasma frequency. This can be improved, however, by recalling that the relevant frequency in strongly coupled systems is the Einstein frequency, which results from the lowest-order sum rule for the velocity autocorrelation function [29]; this frequency accounts for caging and screening of the Coulomb interaction. If we use the inverse Einstein frequency as our time scale, we finally arrive at 1 D ¼ uE a2 : 3

ðA8Þ

By writing the diffusion coefficient in terms of this quantity, we have ensured that the integral of the autocorrelation function is quasiuniversal. Similar arguments follow for other transport coefficients, such as viscosity.

Appendix B. Issues with Z

Now, we must also scale the variables that are implicit in the averaging process, viz. bH

4pnQ2 a2 ; D ¼ ~t M 3

ðA2Þ

0

Z

We now have derived our reduced transport coefficient D ¼ ~t~u=M, and we are free to choose the scales as we wish. Based on Eq. (1), it is reasonable to take the energy scale to be Q2/a; however, this will be modified below to account for screening. Since this energy scale can also be written as 4pnQ2a2/3 we get

ðA3Þ

At the heart of the YVM is the average ionization state Z. Here, the simple TF model [15] has been employed as a rapid way of estimating Z. The strengths of the TF ionization model are that: (1) it is an ‘‘all-electron’’ model in which all electrons (bound and free) are treated on an equal footing, (2) it

where G represents the phase-space variables, b ¼ 1/T, and H is the Hamiltonian, and the evolution operator b vðtÞ ¼ eiLt vð0Þ;

ðA4Þ

which is in terms of the Liouville operator b L [29]. For the case considered here, it is sufficient to treat only the averaging process, since scaling the evolution operator is redundant. Then, we can write H ¼ M~v2

X1 2

i

v2i =~v2 þ ~ u

X

Uij =~ u:

ðA5Þ

i
~ ¼ M~v2 , the Note that, if we choose our energy scale to be u Hamiltonian depends on one dimensionful variable, which in turn leads to ~t~u D¼ M

! ZN dt hvðtÞvð0Þi=~v2 ; ~t

ðA6Þ

0

where now the averaging is in terms of the same dimensionless variables as the explicit variables.

Fig. 6. A comparison of three models for Z ¼ hZi is shown for liquid metals just above the melting point. The dark yellow and green curves result from two definitions using the Purgatorio code, whereas the blue curves arise from the TF model. (For interpretation of colors in this article, the reader is referred to the web version of the article.)

M.S. Murillo / High Energy Density Physics 4 (2008) 49e57

Zedge ¼

4pa3 ne ðaÞ: 3

57

ðB2Þ

Both definitions can be obtained from Purgatorio. A comparison of these definitions is given in Fig. 6; here, the TF results are shown in blue, the Purgatorio results according to Eq. (B2) are shown in green, and the Purgatorio results according to Eq. (B1) are shown in dark yellow. As is well known [30], there are large differences between these definitions. Note that the TF Z values, which were computed according to Eq. (B2), are consistently higher than the corresponding Purgatorio values. The way in which these uncertainties in Z affect the viscosity calculation is shown in Fig. 7. Interestingly, it is the TF result that is closest to the actual data (red dots), suggesting that the good results using that model might be fortuitous. References

Fig. 7. Viscosity predictions are shown using the results from Fig. 6. (For interpretation of colors in this article, the reader is referred to the web version of the article.)

is a true finite-temperature model in that all electrons (bound and free) are thermally excited, (3) it describes finite compression through the boundary condition at the WignereSeitz radius, and (4) all of these features are accounted for selfconsistency. There are, however, many limitations to the TF model of Z. The two most important limitations are (1) TF theory does not properly describe atomic shell structure and (2) the spherical boundary condition cannot be justified on most situations. Because of such limitations, it is useful to understand how sensitive the YVM is to the TF result. In practice, it is difficult to relax all of the assumptions of the TF model, but we can relax the assumption surrounding the shell structure issue. In particular, we can compare two YVMs based on TF and the Purgatorio model [30]. The Purgatorio model is similar to the TF model in nearly all respects, except that it computes the electronic structure through a solution of the Dirac equation, and therefore contains atomic shell structure. The ionization state is an ill-defined quantity in dense matter. It is therefore useful to compare two definitions of Z. The most natural definition is to simply count the number of electrons in positive energy states, according to Z d3 r npos ðB1Þ Zpos ¼ e ðr Þ; WS

where npos e (r) is the density of electrons in positive energy states, and integration is over the WignereSeitz cell. The TF model used here, however, uses the alternate definition:

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

[16] [17] [18]

[19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

M.D. Rutter, et al., Phys. Rev. B 66 (2002) 060102. G.A. de Wijs, et al., Nature 392 (1998) 805. D. Alfe, G. Kresse, M.J. Gillan, Phys. Rev. B 61 (2000) 132. L.J. Wittenberg, D. Ofte, W.G. Rohr, Nucl. Appl. Technol. 3 (1967) 550. F.J. Cherne, M.I. Baskes, B.L. Holian, Phys. Rev. B 67 (2003) 092104. D. Galmiche, S. Gauthier, Jpn. J. Appl. Phys. 35 (1996) 4516. T. Desai, H.C. Pant, Laser Part. Beams 18 (2000) 119. M.S. Murillo, Phys. Plasmas 7 (2000) 33. S.H. Glenzer, et al., Phys. Rev. Lett. 98 (2007) 065002. D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, New York, 1996. J. Cle´rouin, J.-F. Dufreˆche, Phys. Rev. E 64 (2001) 066406. M.S. Murillo, Phys. Rev. E 62 (2000) 4115; G. Faussurier, M.S. Murillo, Phys. Rev. E 67 (2003) 046404. S. Hamaguchi, R.T. Farouki, D.H.E. Dubin, Phys. Rev. E 56 (1997) 4671. T. Saigo, S. Hamaguchi, Phys. Plasmas 9 (2002) 1210. R. More, Atomic physics in inertial confinement fusion. LLNL-Report UCRL-84991, 1981. D. Salzman, Atomic Physics in Hot Plasmas, Oxford University Press, Oxford, 1998. A.A. Louis, N.W. Ashcroft, Phys. Rev. Lett. 81 (1998) 4456. C.S. Yoo, J. Akella, A.J. Campbell, H.K. Mao, R.J. Hemley, Science 270 (1995) 1473. L. Battezzati, A.L. Greer, Acta Metall. 37 (1989) 1791; CRC Handbook of Chemistry and Physics, 86th ed. Taylor and Francis Group, New York, 2005. J. Lindl, Phys. Plasmas 2 (1995) 3933. J. Chihara, J. Phys.: Condens. Matter 12 (2000) 231. A.A. Kugler, J. Stat. Phys. 12 (1975) 35; K. Utsumi, S. Ichimaru, Phys. Rev. B 22 (1980) 1522. S. Ichimaru, Statistical Plasma Physics, , In: Condensed Plasmas, vol. II, Addison-Wesley, New York, 1994, p. 15. T. Scopigno, U. Balucani, G. Ruocco, F. Sette, Phys. Rev. Lett. 85 (2000) 4076. Y. Rosenfeld, J. Phys.: Condens. Matter 13 (2001) L39. J.P. Hansen, G.M. Torrie, P. Vieillefosse, Phys. Rev. A 16 (1977) 2153. Y. Rosenfeld, Phys. Rev. E 47 (1993) 2676. J.G. Cle´rouin, M.H. Cherfi, G. Zerah, Europhys. Lett. 42 (1998) 37. S. Bastea, Phys. Rev. E 71 (2005) 056405. J.P. Boon, S. Yip, Molecular Hydrodynamics, Dover Publications, New York, 1991. P.A. Sterne, S.B. Hansen, B.G. Wilson, W.A. Isaacs, High Energy Density Phys. 3 (2007) 278.