{'fraction_non_alphanumeric': 0.03416239233523662, 'fraction_numerical': 0.009580954224329818, 'mean_word_length': 4.6271555636231625, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiple Scattering Theory for Dense Plasmas

Dense plasmas occur in stars, giant planets and in inertial fusion experiments. Accurate modeling of the electronic structure of these plasmas allows for prediction of material properties that can in turn be used to simulate these astrophysical objects and terrestrial experiments. But modeling them remains a challenge. Here we explore the Korringa-Kohn-Rostoker Green's function (KKR-GF) method for this purpose. We find that it is able to predict equation of state in good agreement with other state of the art methods, where they are accurate and viable. In addition, it is shown that the computational cost does not significantly change with temperature, in contrast with other approaches. Moreover, the method does not use pseudopotentials - core states are calculated self consistently. We conclude that KKR-GF is a very promising method for dense plasma simulation.


I. INTRODUCTION
In stars like our sun [1], in white dwarfs [2], and in inertial fusion experiments [3,4] a state of matter variously known as warm and hot dense matter, or dense plasma, is reached. To understand and model these physical systems it is necessary to be able to accurately predict the material properties, such as equation of state and transport properties, of these dense plasmas.
There has been a serious and sustained effort to produce such accurate models. For temperatures roughly less than the Fermi temperature, plane-wave based Kohn-Sham density functional theory molecular dynamics (DFT-MD) has become a widely used and reliable tool [5][6][7][8][9]. To get to higher temperatures, orbital-free MD has proved useful [10][11][12][13][14], but is limited in its physical accuracy by the approximate kinetic-entropic functional used. For example, we know of no orbital free functional that can predict discrete bound (or core) states [15,16]. Very recently, stochastic DFT [17] has been shown to perform well at high temperature, but still uses pseudopotentials to represent core states, leading to transferability issues. Moreover, stochastic DFT yields the electron density, not the Kohn-Sham orbitals, so standard linear-response theory for optical and transport coefficients can not be applied, and alternative methods such as time-dependent DFT will have to be developed [18].
Another high fidelity method is Path integral Monte Carlo (PIMC) [19,20], which is an excellent method for equation of state (EOS) at high temperature, but becomes prohibitively expensive at lower temperature due to the fermion sign problem. To our knowledge, the method cannot directly produce other quantities of interest like opacity or transport coefficients.
More approximate methods include average atom models and their extensions to include ionic structure [21][22][23][24][25][26]. These are DFT models of one fictitious averaged atom * Electronic address: starrett@lanl.gov whose electronic structure represents an average over all atoms in the plasma. These methods have been enormously useful, and still are the dominant method for production of large tables of data due to their reasonable fidelity and fast computation times. However, due to the extreme difficulty in performing controlled EOS, opacity and transport coefficient measurements on dense plasmas in the laboratory, benchmarking these codes is largely done with reference to more complete models.
In this paper we explore the Korringa-Kohn-Rostoker Green's function (KKR-GF) method, also known as the multiple scattering method, for calculating the electronic structure of dense plasmas [27][28][29]. KKR-GF is a general method for solving the Schrödinger or Dirac equation, but routinely, as here, it is used to solve the Khon-Sham(-Dirac) DFT equations [30][31][32][33][34]. The method has been developed for many years in the solid state community [35][36][37][38], and was introduced to the dense plasma community quite recently [39], where an overview of the method and its potential benefits were presented, as well as some exploratory calculations. Since then, some work has been done exploring the behavior of high temperature crystals, i.e. hot electrons with ions fixed in their crystal positions [40], and some calculations of EOS, again using hot crystals have been presented [41]. However, the challenging question of the behavior of the method for disordered plasma environments has not, until now, been addressed. Here we address this question, and present calculations for dense plasmas. We find good agreement for EOS with existing state of the art methods, where they are accurate and tractable. We find that the key to applying KKR-GF to these disordered environments lies in the concept of extra expansion centers.
In demonstrating that this method is viable for dense plasmas, we are introducing a method with a very high physical fidelity for calculating the electronic structure of dense plasmas. It has the advantage of providing an, in principle, exact solution of the Khon-Sham DFT equations; even core electrons are evaluated self-consistently, in contrast to pseudopotential-based methods. The main weakness, as discussed later in this work, is the conver-gence of the multiple scattering sum. The practical solution, used for many years when applying the method to low temperature systems, is demonstrated here to also give good results for dense plasmas.

II. GENERAL THEORY
The central quantity in the KKR-GF method is the one-electron Green's function G(r, r , z), defined as the solution to the Kohn-Sham inhomogeneous equation [42] [z − H(r)] G(r, r , z) = δ(r − r ) (1) where z is the (complex) energy and H is the Kohn-Sham Hamiltonian. The spectral representation of the Green's function for real energy is where lim η → +0 should be taken, giving the so-called retarded Green's function, and the ψ ν are the eigenstates with eigenvalues ν . Contact with the more familiar wave function formalism is made by taking the imaginary part of the Green's function The electron density is given by where the 2 is due to spin degeneracy, the integral is along the real energy axis, and f is the Fermi-Dirac function with chemical potential µ. Using equation (3), equation (4) clearly recovers the usual expression for n e (r). An important concept in the KKR-GF [27][28][29] method is the multi-site expansion of the Green's function, where the spatial dependence of G(r, r , z) is expanded about a set of fixed centers. These centers are used to define cells that are non-overlapping and space filling. Starting with the known free electron Green's function where p = √ 2mz, one expands the Green's function about a set of centers {R n } using spherical harmonics Y lm (r). Note that we are using atomic units in which h = a B = 4π 0 = e = 1, but we leave the electron mass m dependence as it allows easy conversion between Rydberg and Hartree atomic units. Letting x = r + R n and x = r + R n , we arrive at where is the spherical Bessel (Hankel) function, and G nn 0,LL are the structure constants (appendix A). The superscript × indicates that the complex conjugate of the spherical harmonic should be taken. Also, r > (r < ) means that one should take the vector with greater (lesser) magnitude of r and r .
The expression (6) can be generalized to the case of non-free electrons, giving where R n L is the regular solution to the Kohn-Sham equation in cell n and H n L is the irregular solution, and G nn LL are the structural Green's function matrix elements (appendix A). The double sum term in equation (7) can be thought of as a multisite correction the the single site Green's function, defined by the single summation term. This so-called 'multiple scattering' term corrects the boundary condition applied to the single site solution (i.e. the free electron boundary condition), replacing it with the correct form so that all incoming and outgoing waves are matched at the cell boundaries. Another way to think of it, is that the multiple scattering term introduces quantum diffraction, where electron scattering from multiple sites can interfere. The structural Green's function matrix G is related to the structure constants matrix G 0 and to the so-called t-matrix, which encapsulates the scattering properties of This expression is known as the fundamental equation of multiple scattering theory. It is solved by matrix inversion. Each term is a matrix indexed by the site indices n and n , and by L and L . The t-matrix is site diagonal.

III. PRACTICAL CONSIDERATIONS
In this section we give details on the implementation of this method. A flow chart summary of the method is given in figure 1. At a high level, the method involves dividing space into non-overlapping, space-filling cells, solving the Kohn-Sham equations within each cell independently, and combining these solutions via Dyson's equation to form the Green's function. In turn this Green's function is used to obtain a new electron density to feed back into the Kohn-Sham equations until self-consistency is achieved.

A. Tessellation of Space
To implement the KKR-GF method in practice the first step is to tessellate space. For regular crystals a Voronoi tessellation has been used with success [30]. For crystals with more than one atomic species, using a power tessellation has been shown to be numerically advantageous [43]. In principle the Green's function should be independent of the tessellation, but to minimize the number of spherical harmonic expansion terms, ideally the expansion regions (i.e. the cells) should be as spherically symmetric as possible. This is why the power tessellation is better for mixtures; it allows one to weight cell sizes, in contrast to the Voronoi procedure which tessellates purely on the positions of the expansion centers. In the implementation applied here we have used the powertessellation code of reference [44].
There is another important consideration when making this tessellation: the derivation of equations (6) and (7) rely the following sets of inequalities being satisfied [45,46]: and |r | < |R n − R n | and |r| < |x − R n | These are satisfied for all combinations of r and r if the cells n and n are the largest possible, non-overlapping spheres for a given set of centers; these spheres are called the muffin tin spheres. For space filling tessellations the left hand set of inequalities can always be made to be satisfied by adding extra centers (see figure 2). For power tessellations, one of the right-hand inequalities can also always be satisfied [45]. However, there will always be certain pairs of points where the other right-hand inequality is violated ( figure 3). This leads to slow or even conditional convergence of the double sum in the expansion (7).
The practical solution to this conditional convergence issue is to truncate the double summation to l, l ≤ l ms max [46,47], where l ms max = 2, 3, or 4. This works because using l ms max = 2, 3, or 4 covers the orbital angular momenta of the chemically relevant electrons and bound states. Higher l states correspond to higher energy electron states that are therefore more like free electrons which are less affected by multisite contributions. It is important to note that adding higher order terms in the expansion (larger l ms max ) may actually lead to worse results due to this conditional convergence [45,46]. This is the key weakness of the method. It can be systematically ameliorated by increasing the number of expansion centers. However, as we shall see in section IV B, relative insensitivity of the equation of state to the choice of l max within the range 2 to 4, affirms the validity of this practical solution.
Once the expansion centers have been decided the structural constants G 0 are calculated. These depend only on the positions of the centers and the energies needed (see later). Hence for a given set of centers, they need only be calculated once. We have assumed a periodic supercell in our calculations. This artificial construct is widely used for simulations of disordered systems. The corresponding structural constants are thus calculated in Fourier-space and transformed back to realspace, see appendix A. An alternative method that we do not explore here, is to neglect long range effects on the Green's function, and evaluate the structure constants using a finite cluster of atoms [30,48,49].

B. Single-Site Solutions
The next step is to solve the single-site problem in each cell. This means solving the Kohn-Sham equation for the regular R n L (r) and irregular H n L (r) solutions, as well as the t-matrix for each site t n . This problem is consid-erably simplified if the effective potential inside the cell V ef f,n (r) is spherically symmetric. In cells around a nucleus, the electron density (and effective potential) will be approximately spherically symmetric due to the dominance of the nuclear potential. In cells not containing a nucleus, the electron density will be reasonably close to that of a free electron gas, which is uniform. Hence for this initial evaluation of the method for disordered dense plasmas, it is reasonable to adopt this spherically symmetric approximation for V ef f,n (r) within the cells. This approximation is called the muffin tin approximation.
In the muffin tin approximation the effective potential inside the muffin tin spheres is spherically averaged: where R n mt is the muffin tin radius for cell n. Outside the muffin tin spheres, in the interstitial (is) region, the potential takes an average value -the muffin tin constant. This is determined bȳ where the sum is over all cells, and V s is is the volume of the interstitial region for cell n. Hence the potential inside each cell is spherically symmetric.
An important point worth clarifying is that even though V ef f,n mt (r) is spherically symmetric, the Green's function and electron density inside a cell are not, in general. The symmetry is broken by the multi-site term in the Green's function. This results from the coupling of different spherical harmonics in the double sum term in equation (7).
The regular and irregular solution must be matched to their free electron forms at the edge of the cell where the t n -matrix is now diagonal in the muffin tin approximation. The regular solution is found by integrating out from the origin, and the irregular solution is found by integrating inward [50]. In matching the numerical solutions to these analytic forms the t n -matrix is found. Going beyond the MT approximation is the goal of so-called "full potential" calculations [35,51,52]. This introduces significant complexity to the single-site problem. For example, the t n -matrix ceases to be diagonal in L. The importance of full-potential versus MT for disordered plasma applications will be considered in the future.

C. Efficient Evaluation of Multi-Site Effects
With the single-site problem solved and the structural Green's function evaluated, one can construct the Green's function, equation (7). The multiple scattering term is expensive to evaluate due in part to the double sum (equation (7)) but also because the solution of Dyson's equation for the structural Green's function matrix involves a matrix inversion for each energy and kspace integration point (appendix A). The size of this matrix is N (l ms max + 1) 2 where N is the number of expansion centers and l ms max is the maximum l for both summations. Since dense matrix inversion cost scales as size of the matrix cubed, it is clearly desirable to keep N and l ms max as small as possible 1 . As noted above, only the first few terms are included in the double sum (l ms max = 2, 3 or 4), which serves to limit computational expense. As explored in reference [40], higher l terms are kept in the single site part of the Green's function. This is expected to be a good approximation because higher l states correspond to higher energy states, and hence are more free-electron like, leading to small multiple scattering contributions.
Another way to reduce computational expense is to evaluate the multi-site term only in the energy range where it is significant. We expect tightly bound states to be unaffected by multi-site interactions due to their localized nature. At the other extreme, electrons in highenergy states should be insensitive to the details of multisite interactions due to their free electron-like character. Hence we introduce a minimum E min and maximum energy E ms max for multi-site contributions, outside of which only the single-site term contributes.
Even if the multiple scattering Green's function were evaluated only in the range E min to E ms max , it would be a rapidly varying function of the energy there, requiring many integration points to accurately resolve the multiple scattering contribution to the density, equation (4). This is avoided by deforming the energy integration path into the upper half of the complex plane. This technique results in much smoother integrands, fewer integration points, and most importantly avoids the poles on the real axis, which correspond to the eigenstates of the system. Hence, one does not actually find these eigenstates, nor are they needed. The (retarded) Green's function has no poles in the upper half complex plane. However, the Fermi-Dirac function does have poles, at the so-called Matsubara frequencies: z = µ + ıπ(2j − 1)k B T , for j = 1, 2, . . .. Their residues can be evaluated directly, or avoided by choosing a contour with imaginary part of energy everywhere less than πk B T . Note also, that the poles are avoided if µ < E min or µ > E ms max . This strategy is summarized in figure 4. For large negative energy only the single-site term is evaluated, and and any core states that exist are accounted for using equation (3). Between E min and E ms max the full Green's function expression is evaluated using contour integra-1 Note, however, that there exist so-called 'screening transform' methods which sparsify the matrix inversion and permit more favorable scaling for large systems [53]. tion. Note that, for this region l ss max for the single-site term is as large as needed for convergence, whereas l ms max for the multi-site term is fixed at a small values (e.g. 2, 3, or 4). Above E ms max only the single site term is evaluated, up to some maximum energy E max , determined by integrand having become negligibly small.
With the Green's function determined, the electron density can be evaluated, and the usual self consistent field procedure is followed. The Poisson equation is solved using the technique presented in [54]. The selfconsistent field problem is solved using Eyert's method [55], and normally takes less than 10 iterations to converge, using the Thomas-Fermi cell model to generate an initial guess [56]. For spatial integrations over the cells, we use the isoparametric integration method of reference [57]. Note also, that in the calculations presented here we use the temperature dependent local density approximation [58].

D. Thermodynamics
To study the convergence properties of the KKR-GF method in disordered plasmas, we focus on thermodynamic properties. These require the electron density of states in addition to the density. The density of states per volume V is also found from the Green's function The pressure is given by the Virial expression where F el is the electrostatic free energy [54], P I is the ideal ion pressure, and U k is the internal kinetic energy . (17) and assuming the local density approximation where F xc is the exchange and correlation free energy, Evaluating the plasma equation of state also requires averaging over an ensemble of realistic ion positions. We generate these using the PAMD (pseudo-atom molecular dynamics) model [59]. This model uses a DFT average atom calculation to produce a force field for use in classical molecular dynamics. It gives an accurate set of ion transport coefficients and pair distribution functions for dense plasmas [60,61]. It is numerically inexpensive and we use it to generate sets of ion positions that are well separated in time and hence, are uncorrelated. This is a reliable and accurate procedure for dense plasmas. For lower temperature liquids and warm dense matter with transient chemical bonds, PAMD will not produce accurate ion positions. In principle, it is possible to calculate forces on ions directly from the KKR-GF solution [62]. However, we have not attempted that yet and we note that results shown are firmly in the regime where PAMD is accurate.

A. The Effect of Extra Expansion Centers
In figure 5 we show the total pressure as a function of the number of extra centers for two different temperatures, for aluminum at solid density. For each temperature, we show results for two different sets of ion positions. We are using 8 nuclei in the simulation volume. The number of extra centers corresponds to the number of cells not containing a nucleus. We see that on adding centers the pressure initially changes rapidly and then settles down to a rough constant. There are still some fluctuations in the pressure even after the value settles down. The pressure initially changes rapidly because the inequities (9) are strongly violated. As extra centers are added, the expansion becomes valid and the value settles down. However, as we are using the muffin-tin approximation, adding extra centers treats the electrons in the interstitial region slightly differently to before, hence the fluctuations. At the same time, the region being modeled with the conditionally convergent double summation also reduces in size. The exact interplay of these errors is not yet known but could be examined in detail with a future full-potential treatment.
It is also noted that the 100 eV case reaches its (roughly) constant pressure more quickly than the 10 eV case. We can understand this behavior by noting that at high temperature the multiple scattering effect is of less relative importance as there are more high energy electrons that are insensitive to multi-site effects.  where we have used l ms max = 3. At 10 eV we compare with DFT-MD results [63] and find agreement within the error bars. We also compare with PAMD [64] which also gives reasonable agreement, and the average atom   model Tartarus [23], which gives too low a pressure. Since PAMD is based on the idea of correcting average atom pressures for ionic disorder, the difference between them implies that the lack of this effect in the Tartarus model is the cause of the difference. At 100 eV the relative change with l ms max is smaller than at 10 eV. This is due to the decreased relative importance of multi-site effects, as noted earlier. Agreement with PAMD is again good.
The rather large fluctuations in the pressure seen in figure 6 are due to the small number of nuclei in the simulation, i.e. eight nuclei. There is no restriction of the method to small numbers of atoms. Indeed, linear scaling variations exist that could be exploited for plasma applications [37,49,62]. However, our results are enough to demonstrate the viability of the method, our main aim. Table I shows total pressure along the 2.7 g/cm 3 isochore from 10 to 1000 eV and the wall time for one complete self consistent field cycle (to convergence) for a given set of ion positions. The key point is that, in practice, the wall-time for KKR-GF method does not significantly change with temperature, in stark contrast with plane wave DFT-MD [15]. Also shown is the wall-time for calculation of the structural constants, which takes roughly 5/6 of the total execution time. This can probably be dramatically improved, either through algorithmic improvements, exploitation of modern parallelism, or switching to real-space approximations [48,49].

C. Diffraction and Disorder Effects on Electronic Structure
In figure 7 the density of states per atom is shown for the 10 and 100 eV aluminum cases, focusing on the valence region where multi-site, quantum diffraction effects are most important. The density of states gives detailed information about the electronic structure, and is a less averaged quantity than the EOS. To produce the density of states from the KKR-GF model, we have evaluated the Green's function slightly above the real energy axis. This leads to a broadening of the the density of states [36], equivalent to convolving with a Lorentzian of fullwidth-half-maximum equal to the imaginary part of the energy (here set to 0.1 E H ). This is similar to the artificial broadening used in plane-wave DFT-MD to produce densities of states.
We compare to the average atom density of states using the Tartarus model, which is a spherically symmetric, single center model of an average atom in a jellium. It has only one atom in it and hence no multiple scattering. At 10 eV the Tartarus model yields two bound states in this energy range (2s and 2p), which are represented graphically by vertical lines corresponding to a Diracdelta density of states. In the KKR-GF model, these states broaden into bands of states, at roughly the same energy as the Tartarus bound states. This broadening is the result of two effects. First, since each atom has a different local environment, the eigenenergy varies from site to site. Second, the inclusion of quantum diffraction through the multi-site term broadens the state, as was seen for high temperature ordered lattices [40]. The broadening can crudely be interpreted as a reduction in lifetime of the state due to multiple scattering, an effect that is absent in the average atom model.
For the average atom model, a distinct continuum in observed (figure 7) which corresponds to states with energies greater than the average atom muffin tin level. These states are the charge carrying electron states in conductors. At 10 eV, this continuum is also observed in the KKR-GF result and is more like the free electron form (∝ √ ). The artificially high symmetry in the average atom model leads to the structure in the continuum not seen in the KKR-GF model.
At 100 eV, the average atom model predicts two bound states in the energy range shown (figure 7), as well as a resonance structure in the continuum. The resonance structure in average atom models is well documented [22,23] and is a long lived quasi-stable state, caused by a potential minimum due to the sum of the effective potential and the centrifugal term l(l + 1)/r 2 in the radial Kohn-Sham equation. These bound and resonance states merge into one 'bump' in the KKR-GF result, and the continuum merges with the bound states. Resonances could occur in the KKR-GF model but are likely to be more spread out due to the fluctuations in the local enviroment and the multiple scattering effect. This merging of the average atom bound and resonance states clearly has an influence on the EOS (see earlier). Perhaps more dramatic would be the influence on optical spectra, as lines corresponding to transitions to and from these states in the average atom model would disappear and be replaced with some sort of merged feature. This is the subject of future work.

V. CONCLUSIONS AND DISCUSSION
The KKR-GF method has been demonstrated to be a viable and promising modeling tool for dense plasmas. Calculations of equation of state demonstrate that the method gives good agreement with other state of the art methods, where they are valid and are practical. Two major advantages of the method are 1) no pseudopotential is needed; core states are calculated explicitly and self-consistently, and 2) example calculations show that the computational time does not significantly change with temperature.
There are limitations to our current implementation that are not limitations of the method in general. First, we have used the muffin tin approximation. This does not appear to cause major inaccuracies for the conditions explored in this work (i.e. temperature > 10 eV), but we expect it to be more problematic at lower temperature where there only lower energy electrons, that will be more strongly affected by inaccuracies in the potential. We note that the muffin-tin approximation is unnecessary, but serves us as a useful test vehicle. Second, we have used a small number of atoms in the simulation cell (8), due to computer memory limitations, but a more sophisticated implementation could overcome this, possibly simulating thousands of atoms [37,49,62]. Third, we have imported the ion positions from a different model (PAMD [50]). While this is accurate for the present results, molecular dynamics forces could be calculated directly from the electronic density, as is done in plane-wave DFT-MD.
A limitation of the method that has yet to be fully resolved is the convergence of the multi-site contribution to the Green's function. Recent work has explored the possibility of transforming the conditionally convergent multi-site summation to an absolutely convergence one [46], but to our knowledge, the practicality of this approach remains an open question. We have used the empirically supported method of keeping only the chemically relevant expansion terms, and found only weak dependence of the pressure on the exact number, consistent with past results.
An exploration of the density of states of hot, dense aluminum plasmas reveals significant changes to the density of states compared to an average atom model. There is broad alignment of features with this more approximate model, and the comparison facilitates understanding how multiple scattering affects the electronic structure. It is pointed out that this will have significant impact on spectroscopic lines, for example, in the absorption coefficients.
In this appendix we give the expression for the structural Green's function G nn LL for a periodic system [65,66]. The structural Green's function is first found in Fourier space using Ewald's technique and the transformed back to real space. For the application presented here, we consider a periodic supercell in a simple cubic lattice. Let µ be the index the lattice vectors a µ , and R n be the vector pointing from the origin to the n th atom in the supercell. The real space structural Green's function where the matrices are indexed in n and L, and the integral is over the Brillouin zone of volume V BZ . Note that, in the main text, the omission of the superscripts µ and µ indicates that µ = µ , referring to atoms in the same supercell. The structural constants are given by with k i being a reciprocal space lattice vector, the sum is over all these lattice vectors, and V is the volume of the supercell. D ×Y * lm ( a s − R nn )|a s − R nn | l × ∞ √ η/2 ξ 2l exp[−ξ 2 (a s − R nn ) 2 + (z/(4ξ 2 ))]dξ (A7) where R nn = R n − R n and the sum is over all lattice vectors, with the prime indicating that the lattice vector a s = 0 is ommitted.
η is the Ewald parameter and the result is independent of its value (though computation time is not [67]). In implementing these equations, we found the work of reference [68] helpful.