Electronic-structure calculations for nonisothermal warm dense matter

John Jasper Bekx,1,2,* Sang-Kil Son ,1,3 Beata Ziaja ,1,4 and Robin Santra 1,2,3 1Center for Free-Electron Laser Science, Deutsches Elektronen-Synchrotron DESY, Notkestrasse 85, 22607 Hamburg, Germany 2Department of Physics, Universität Hamburg, Jungiusstrasse 9, 20355 Hamburg, Germany 3The Hamburg Centre for Ultrafast Imaging, Luruper Chaussee 149, 22761 Hamburg, Germany 4Institute of Nuclear Physics, Polish Academy of Sciences, Radzikowskiego 152, 31-342 Kraków, Poland

WDM lies on the border between condensed-matter physics and plasma physics [16]. It is characterized by temperatures of the order of T ∼ 0.1-100 eV and densities ranging quantum kinetic theory [44,45], and quantum Monte Carlo simulations [46,47]. A recent review highlighting and analyzing the last two topics can be found in Ref. [48].
In the current study, we will focus on the description of WDM created by an intense femtosecond XFEL pulse from solid aluminum (Al). Specifically, we will limit ourselves to the description of the electronic states during the first 10-100 fs after the x-ray exposure, which allows us to assume that the ion lattice remains cold. This assumption is well justified as the heating of electrons by the laser pulse occurs much faster than the transfer of energy to ions through electron-ion coupling [49][50][51][52][53][54][55]. Additionally, the dense electronic environment created forces the electrons to rapidly thermalize (on the order of a few fs to a few tens of fs) through electron-electron collisions and impact-ionization processes [56]. These assumptions allow us to describe the system as a hot, thermalized electron distribution embedded in a periodic crystal lattice, constituted of parent ions. The system is modeled by a muffin-tin potential, whereas the description of electronic states is done using a finite-temperature Hartree-Fock-Slater approach. These states are represented with respect to a hybrid basis consisting of plane waves and localized core orbitals solved on a radial pseudospectral grid. A simultaneous interleaved optimization of inner-shell and outer-shell electronic states accounts for inner-shell modifications at high electronic temperatures in an efficient manner. Additionally, the periodicity of the system allows for the implementation of the Bloch formalism [57], increasing computational efficiency. A detailed description of the developed methodology is presented in Sec. II. The model has been implemented into a new toolkit, XCRYSTAL. With the calculated electronic states, various characteristic properties of the transient, nonisothermal WDM state become available for evaluation and for comparison with experimental data. In order to benchmark XCRYSTAL, we will compare our model results to the data from the experiment performed at the LCLS on x-ray excited solid Al [18][19][20]. Specifically, the K-shell threshold energies were measured for this system. From them, the ionization potential depression (IPD) was determined, i.e., the lowering of the ionization potentials of atoms present in the system caused by the dense and charged environment. It was used to check the predictive capabilities of the widely-used EK and SP models [18,19,39]. The predictions of the SP model were found to be unsatisfactory, and a modified EK model [19,58] was proposed for fitting the data.
In what follows, we will provide XCRYSTAL predictions for the K-shell threshold energy and IPD for Al, and compare them to experiment and the aforementioned theoretical models. In addition, we will compare the results from XCRYSTAL with those obtained from the AA model and the two-step HFS model, introduced in Ref. [29]. Because of the incorporation of the periodic crystal structure within our model, it is possible with XCRYSTAL to calculate the band structure of the system at very high electronic temperatures (∼100 eV). Let us also note that, even though band structure calculations at finite temperatures have been performed before [59,60], it has never been done in the regime of such high temperatures for the type of transient WDM state we are considering. Also, the WDM system that we are interested in is strongly nonisothermal. It is characterized by electronic temperatures far above room temperature (up to 10-100 eV or 10 5 -10 6 K) and ion temperatures close to room temperature. Therefore, the finite-temperature band structure calculations that we obtain at these conditions will be the first of their kind.
The paper is structured as follows: In Sec. II we outline the theory framework. In the following Sec. III, we provide the K-shell threshold energies and the IPD values for soliddensity aluminum (Al) at electron temperatures ranging from 0 to 100 eV. We then compare them with experiment and various theoretical models. Finally, we present the temperaturedependent band structure calculations of this system at various temperatures, and comment on the trends observed. In Sec. IV we provide conclusions and an outlook.

II. METHODOLOGY
In this section we outline our theory framework. In particular, we describe the treatment of electrons in our approach. We distinguish between (1) electrons highly localized around an ion (core electrons) and (2) electrons delocalized within the unit cell (valence electrons). Worthy of note is that in our framework we will update the core electrons to the presence of the plasma electrons as opposed to keeping them frozen as is custom in low-temperature condensed-matter calculations. Our approach consists of describing the electronic states as Bloch states, represented in a hybrid basis consisting of plane waves and these updated localized core orbitals. Throughout this paper, atomic units are employed, i.e., m e = e =h = 4π 0 = 1. We also set k B = 1. The unit of length will be given in units of the Bohr radius, a 0 , and energies and temperatures will be given in electronvolts (eV).

A. Schrödinger equation for a periodic system
We consider a perfectly periodic crystal. The crystal lattice is defined by the primitive lattice vectors a 1 , a 2 , and a 3 , thus defining a primitive unit cell volume = |a 1 · (a 2 × a 3 )|. In the independent-particle approximation, the wave function of an electron in a periodic system can be represented as a Bloch wave: where a total volume V = N cell is introduced to ensure periodic boundary conditions, with N cell the number of primitive cells. Furthermore, k is the wave vector of the electron, n is a band index, and μ n,k (r) is a lattice periodic function, i.e., μ n,k (r + a i ) = μ n,k (r), with i = 1, 2, and 3. Both labels k and n specify an eigenstate of the effective one-electron (mean-field) Hamiltonian,Ĥ = [− 1 2 ∇ 2 + V (r)], where V (r) is the mean-field potential of the periodic system. Since μ n,k (r) is periodic, we can decompose it into a Fourier series. In a Bravais lattice, this corresponds to a Fourier series with respect to a sum over reciprocal lattice vectors K = 3 i=1 β i b i , with β i being integer coefficients and b i being the reciprocal lattice basis vectors, which obey b i · a j = 2πδ i j .
The Fourier decomposition is thus where m = (β 1 , β 2 , β 3 ) runs over all reciprocal lattice vectors K m and the latter integral is restricted to a unit cell. This implies that the wave function is expanded in a plane-wave basis of the form After inserting Eq. (3) into the Schrödinger equation, H ϕ n,k (r) = E n,k ϕ n,k (r), we arrive at its representation in k space for a periodic system: where i and j denote specific reciprocal lattice vectors and W (K) = 1 d 3 re −iK·r V (r). So far, this representation is general as all that was assumed to derive it was a perfect crystal structure. This allowed us to use a lattice periodic potential V (r), decomposable in a Fourier series over the reciprocal space, and to describe the electrons with a oneelectron wave function specified by Eq. (1). This implies that we neglect correlation effects.

B. Hybrid basis
A direct attempt to solve Eq. (4) would be computationally extremely costly as an accurate representation of the highly localized core states in a basis of plane waves would require the contribution of large wave vectors. This would result in a very slow convergence of the solution. Many different methods of overcoming this problem, typically involving the separation of the localized core electrons from the delocalized valence electrons, have been developed over the years. They include the utilization of Wannier functions [61][62][63], pseudopotential methods [64,65], and the augmented plane wave (APW) method [66][67][68].
Similar to the APW method, we intend to extend the basis of the plane waves so as to include atomic orbitals. Naturally, we do not wish to include all atomic orbitals in our hybrid basis. Only those that are highly localized, i.e., for which the description using only plane waves is inefficient, are added to this hybrid basis. We will refer to them as core orbitals. This provides an ansatz for our wave functions in the form where m runs from 1 to some cutoff index N K , which represents the number of valence bands. The index n C denotes the band index of the core wave functions, introduced as follows.
Since ϕ n,k (r) is a Bloch function and the first term of Eq. (5) is structured so as to be a Bloch function [like in Eq. (3)], the core wave functions ψ n C ,k (r) must be Bloch functions. We can construct them using the conventional atomic orbitals that are eigensolutions of the Hamiltonian, with one Bloch function per core orbital, per atom in the unit cell. This implies that the general label n C includes information regarding which atom in a unit cell we are considering, as well as the electronic structure information, i.e., n C = (a, n, l, m, s), where a denotes an atomic index, and n, l, m, and s denote the principal, azimuthal, magnetic, and spin quantum number, respectively. Subsequently, n C runs from 1 to the total number of core orbitals, N C , in the unit cell. Suppose that φ n C (r) denotes an atomic orbital that is strongly localized (for a detailed explanation of "strong localization" see Sec. II C) and can be considered as a core orbital. Let R I,n C denote the position of the nucleus in the Ith unit cell at which the spatial orbital of the n C type is located: R I,n C = R I + r n C , with R I being the lattice vector denoting the Ith unit cell, and r n C = R 0,n C . The associated atomic core orbital can then be written as φ n C (r − R I,n C ). Since the unit cells are identical, the φ n C (r − R I,n C ) for all the unit cells I are degenerate eigenfunctions of the mean-field Hamiltonian, as they all have the exact same eigenvalue E n C as a consequence of the assumed strong localization. Therefore, any linear combination of them is also an eigenfunction of the same Hamiltonian. In particular, we can choose a linear combination such that the resulting wave function is a Bloch wave ψ n C ,k (r). We can achieve this by considering the following linear combination of atomic spatial orbitals: where N = N a N b is a normalization constant, with N a = 1/ √ V and N b = √ (see Appendix A for the derivation). In order to prove that ψ n C ,k (r) constructed above is indeed a Bloch wave, we must show that μ n C ,k (r) exhibits the same periodicity as the lattice. Let R denote some lattice displacement: where we used the fact that all cells are identical, resulting in a periodic behavior of the wave functions φ n C (r − R I,n C ). Furthermore, we define the index J such that R J,n C = R I,n C − R. Since I and J run over all unit cells, they contain the same terms.
Having constructed our hybrid basis, we may insert our ansatz, Eq. (5), into the Schrödinger equation. This For the core wave functions, we usedĤ ψ n C ,k (r) = E n C ψ n C ,k (r), where the energy eigenvalue of the core electrons, E n C , is independent of k, because these core electrons are highly localized. We subsequently project this onto our hybrid basis and solve the Schrödinger equation through a matrix diagonalization. Performing a projection onto the plane-wave part of our basis, we calculate k + K i |Ĥ|ϕ n,k . Keeping in mind the proper normalization of our plane waves, r|k = e ik·r / √ V , and using V d 3 re i(K i −K j )·r = V δ i j , this yields Performing a projection onto the core-wave-function part of our basis, we calculate ψ n C ,k |Ĥ |ϕ n,k . By using the fact that ψ n C ,k is an eigenstate of the HamiltonianĤ , with the eigenvalue E n C , we find that ψ n C ,k |Ĥ|ϕ n,k = E n C ψ n C ,k |ϕ n,k , or We can combine these equations into a generalized eigenvalue equation, similar to the Roothaan-Hall equations [69,70] of the form FC = SC , where = diag(E n,k ). The nth column of C is . . .
and F and S are with i and j denoting rows and columns, respectively. Note that we have one such equation FC = SC per k point. It should come as no surprise to encounter equations similar to the Roothaan-Hall equations as these are a representation of the Hartree-Fock equations in a nonorthonormal basis, precisely as we have done with our hybrid basis consisting of plane waves and atomic orbitals. Note that the theory developed in Sec. II B assumed nothing more than a periodic potential V (r), and a one-electron description for the electron wave function under considera-tion, with the ansatz given in Eq. (5). At this point we would like to emphasize that, even though we will not consider it in this work, it is possible to include molecular dynamics in this framework, so as to extend the applicability of XCRYSTAL to beyond the timescale of ∼100 fs after x-ray irradiation. Conceptually, we need not take the assumption of periodicity to imply the use of a crystallographic unit cell as is done in this work. A larger supercell containing many more atoms could be employed. However, this would cause the Brillouin zone to shrink, rendering the matrices that must be diagonalized much larger. Under these conditions, employing the Bloch theorem in combination with a plane-wave basis may no longer be the most efficient approach. We leave this development and investigation for a future endeavor.

C. Calculation of atomic orbitals
In order to calculate the atomic orbitals needed for the construction of ψ n C ,k (r), we used the XATOM toolkit [71][72][73]. This toolkit describes an isolated many-electron atom in a Hartree-Fock-Slater (HFS) framework. This is an independent-particle approximation with a mean-field Hamiltonian of the form where Z is the nuclear charge. The exchange potential V exc (r) is assumed to be of the form [74] V exc (r) = − 3 2 with being the electron density, and φ i (r) denoting a single-particle spin-orbital wave function. In addition, XATOM imposes spherical symmetry, assuming the solution of the Schrödinger equation to be of the form with n, l, m, and s being the principal, azimuthal, magnetic, and spin quantum numbers of the electron with the associated wave function φ nlms (r), respectively. A spherical averaging is done on the electronic density. This implies that both the electronic density and potential are spherically symmetric, such that ρ(r) = ρ(r) and V (r) = V (r). Subsequently, XATOM solves the Schrödinger equation in a self-consistent way. There are numerous computational input parameters required in XATOM calculations [71], the most relevant of which are the following. The radial coordinate r in Eq. (16) is defined with the generalized pseudospectral method on a nonuniform grid [75]. This grid is characterized by the number of radial grid points N r , the maximum radius, r max , and the mapping parameter L, which determines the density distribution of radial points [75]. A larger L invokes that more radial grid points are being pushed towards higher values of r. For the computations of the atomic orbitals presented further, we used N r = 200, L = 10, and r max = 5.0 a 0 , unless specified otherwise. Since XATOM calculates the atomic orbitals in an isolated-atom description, these orbitals will not be accurate representations of orbitals in a charged and dense environment of atoms. We show how we adapt the core orbitals to the environment in Sec. II F. Let us continue by explaining exactly how we define the core atomic orbitals φ n C (r) in the framework of XCRYSTAL.
As we intend to work with the atomic orbitals provided by XATOM, which are calculated in a sphere of size r max , we also construct spheres around the constituent atoms located within our unit cell. The radius of the sphere for atom a will be denoted by r s (a). In order to be able to maximally exploit spherical symmetry, we assume the spheres to be touching each other. The values of r s (a) are then dependent on the crystal geometry. From all orbitals obtained with XATOM, we include only those orbitals into our hybrid basis which are localized within the touching spheres. Thus, we define our core orbitals as those that satisfy where δ C is a cutoff parameter set to 10 −3 for all calculations shown in this work. This implies that at least 99.9% of the norm of the core orbitals of atom a is confined to the sphere of radius r C (a). The radius 0 < r C (a) r s (a) is a modeldependent parameter. The reason for the inclusion of r C (a) is that in the upcoming Sec. II D, we assume a model to describe our periodic crystal where we disembark from static, touching spheres.

D. Construction of the crystal potential
Having conceptualized touching spheres of radii r s (a) around our atoms, we may describe our system using a muffintin-type potential [29,66,67]. Inside of the spheres, a spherically symmetric potential V a (|r − R a |) is assumed, where R a is the position vector of the ath atom. Outside of the spheres, in the interstitial region, the potential becomes a constant, V 0 . Typically, this potential is set to zero [25,66,67,[76][77][78]. However, we will follow Ref. [29] in explicitly calculating V 0 . In Ref. [29], V 0 is the potential value at the Wigner-Seitz radius and denotes the energy value above which the continuum of states starts. In our framework, we do not have a similar physical meaning behind V 0 . What remains true is that a continuum of delocalized states will be present at energies above V 0 in XCRYSTAL. However, there may also be some delocalized states below V 0 , evidently referring to states that are somewhat localized, yet not enough to be considered core orbitals as they do not satisfy Eq. (17). To determine the value of V 0 , we require the value of the potential on the boundaries of the spheres to be equal. We know that the potentials in the spheres are ∼C/r for some C < 0. Therefore, for all our touching spheres, we assign V 0 = min{V a (r s )} and proceed by shrinking all spheres to a radius r V (a) so that V a (r V ) = V 0 . In this framework, this implies that in r C (a) = r V (a) in Eq. (17).
The reason for adopting the muffin-tin potential into our framework is twofold: (1) the implementation of a spherically symmetric potential inside of a sphere can be accurately captured by the available XATOM toolkit, and (2) the muffin-tin potential simplifies the evaluation of W (K) = integral to a one-dimensional one: The inclusion of the muffin-tin crystal potential into Eq. (18) allowed us to split up the integral over the unit cell into two parts: the atomic region and the interstitial region, as schematically illustrated in Fig. 1. In the atomic region, the integral is represented as a sum of spherical integrals (in one dimension), whereas in the interstitial region the integral can be evaluated analytically. The latter is highly beneficial for computational efficiency. The spherically symmetric atomic potential V a (r), used in Eq. (18), is given by Given the limitations on accuracy imposed by the simple muffin-tin approximation, using a higher quality exchangecorrelation functional than the one alluded to in Eqs. (15) and (19) is not warranted. We have made the assumption that the electronic density ρ is also spherically symmetric inside of the spheres. The respective spherical averaging was done as ρ(r) = d r ρ(r)/ d r and may be evaluated analytically. Note that the assumption of a spherically averaged electronic density may be a poor approximation for capturing the electronic behavior at low temperatures for materials that exhibit directional bonding. However, we may expect any directional dependencies to be suppressed once the electronic temperature exceeds the zero-temperature band gap, allowing for thermal excitations. This will make the approximation of a spherically averaged electronic density less crude at these temperatures.

E. Electronic density at finite temperature
Within a finite-temperature HFS framework, the temperature is introduced assuming Fermi-Dirac occupation of the electronic orbitals [79]. Therefore, at temperature T and chemical potential μ, the electron density ρ(r) is calculated as p |ϕ p (r)| 2n p (μ, T ), with p being a spin-orbital index, ϕ p the electron wave function, andn p (μ, T ) a Fermi statistical weight. For the atomic core orbitals, the index p = n C (as defined in Sec. II B), and for the delocalized valence orbitals p = (n, k, s). We assume a degeneracy in the energies with regard to the spin quantum number, so a sum over the spins s will simply lead to a factor of 2. The Fermi weight at T is where ε p is the energy eigenvalue associated with ϕ p , β = 1/T , and μ is the chemical potential of the system. The chemical potential is calculated by imposing charge neutrality on our unit cell. To find μ, we must solve where we used the fact that ϕ p (r) is normalized to 1 in the total volume V . Since |ϕ p (r)| 2 is a periodic function, ϕ p (r) is normalized to 1/N cell in the unit cell .

F. Self-consistent-field method
In the self-consistent-field (SCF) method, we start with an initial guess for the spherically symmetric electronic density ρ(r). The initial guess can be based, for example, on the converged XATOM result for the initial atomic orbitals. This is the typical start for calculations at zero temperature. At nonzero temperatures, we may use the density calculated at a lower temperature, or even extrapolate the results from several lower-temperature runs.
With the initial density, we construct the potentials V a (r) as in Eq. (19), and perform the sphere shrinking in order to calculate V 0 (explained in Sec. II D). At this point, instead of keeping our core orbitals fixed, we use the newly obtained muffin-tin potential in a single XATOM diagonalization in order to update the atomic core orbitals. Therefore, during the SCF iteration, the core orbitals, which are incorporated in our hybrid basis, are always updated to the presence of valence electrons, in contrast to the conventional low-temperature condensed-matter calculations where frozen core orbitals are typically used. We do expect that in dense plasmas the core orbitals strongly react to the presence of plasma electrons, throughout the entire unit cell. After this core-orbital update, we proceed by using our previously obtained muffin-tin potential and Eq. (18) to obtain W (K). It is then utilized in the Fock matrix F in Eq. (12). Therein the quantities referring to the core orbitals, E n C and ψ n C ,k (r), are the results from the updated core orbitals. Subsequently, the equation FC = SC is solved to obtain { , C}.
From the solutions obtained for the core orbitals and the valence orbitals, we calculate the new electronic density ρ(r) = ρ core (r) + ρ val (r), which is spherically averaged analytically. With this electronic density, we construct a new potential V a (r) and update V 0 , thereby closing the loop in the SCF method. The whole procedure is repeated until the result converges. Further discussion on the scheme used to accelerate the SCF convergence can be found in Appendix B. An illustration depicting the SCF algorithm employed is shown in Fig. 2.

G. Numerical parameters for crystal calculations
In addition to the computational parameters related to atomic calculations, there is a couple of numerical parameters utilized in XCRYSTAL, which are related to the crystal structure. They are (1) the cutoff on the number of reciprocal lattice vectors used and (2) the numerical grid defined for the momentum vectors in the first Brillouin zone (BZ).
The sum over reciprocal lattice vectors K i , used in the Fourier series should in principle be infinite and include all possible vectors {K i }. For practical purposes, a cutoff on the number of reciprocal lattice vectors, denoted by N K (in Sec. II B), is indirectly imposed through a cutoff on the norm of the vectors, denoted by |K| max . As can be seen from the coupling term W (K i − K j ) in the Fock matrix F in Eq. (12), the quantity |K| max determines the energy cutoff for the bands. At T = 0 eV we may expect to achieve convergence with a relatively low |K| max , as the absence of thermal excitations implies that higher-lying bands will remain unoccupied. However, by increasing the temperature, the number of thermal excitations will rise, and a higher cutoff on |K| max will be necessary to attain convergence. As one increases |K| max , additional plane waves are added to our hybrid basis with increasingly higher momenta. It may occur that these are high enough so as to describe the atomic core orbitals in our basis, and a linear dependency will ensue. In this case, a simple rank reduction is used on the matrices F, S, and C to remove this linear dependency. For consistency between results at various temperatures, all XCRYSTAL calculations were done using |K| max = 6.0 a −1 0 , which corresponds to N K = 1647 valence bands.
The second parameter of interest defines the numerical grid defined for the momentum vectors in the first Brillouin zone and affects the sum over momentum vectors k. Because of our employed Bloch formalism, the sums over k are restricted to the momentum vectors confined within the first Brillouin zone (BZ) of the reciprocal lattice. However, contrary to the reciprocal lattice vectors {K i }, the vectors k in the first BZ constitute an uncountable set. Therefore, a sum over vectors k is performed as an integral: BZ To perform this integral numerically, a grid is defined inside the first BZ. Through a dedicated analysis, we have found a k grid of 7 × 7 × 7 to be sufficient for achieving convergence. A more elaborate discussion on the employed BZ will be given in Sec. III.

III. RESULTS
In this section we show results obtained with XCRYSTAL for x-ray-excited solid Al, investigated experimentally in Refs. [18,19]. This experiment was conducted at the LCLS [18,19] for solid Al at a density of 2.7 g/cm 3 . The sample reached electronic temperatures in the range of T = 10-80 eV. The Kα fluorescence signal was measured as a function of the incoming photon energy. After comparing the results calculated by XCRYSTAL to experiment, we will show band structure plots of this system. At zero temperature, these can be compared to known, well-established result [80,81], whereas the band diagrams at higher temperatures will be of interest with regard to how they evolve as temperature is increased to values characteristic of WDM conditions (up to 100 eV or 10 6 K). To the best of our knowledge, these temperature-dependent electronic band structures at such high temperatures are the first of their kind for this nonisothermal, transient WDM state.
The Al crystal has a face-centered-cubic (FCC) structure, for which the primitive unit cell can be constructed from the primitive lattice vectors a 1 = a(1/2, 0, 1/2), a 2 = a(0, 1/2, 1/2), a 3 = a(1/2, 1/2, 0), where a = 404.95 × 10 −12 m = 7.652 a 0 is the lattice constant of Al. The first Brillouin zone (BZ) is a truncated octahedron. Due to the implementation of the Bloch formalism, every sum over momentum vectors k can be limited to this first Brillouin zone. However, for modeling this system in XCRYSTAL, it was numerically more efficient to define the unit cell to be cubic, with lattice vectors a 1 = a(1, 0, 0), a 2 = a(0, 1, 0), a 3 = a(0, 0, 1), and four Al atoms placed on the positions a(0, 0, 0), a(1/2, 0, 1/2), a(0, 1/2, 1/2), and a(1/2, 1/2, 0), similar to what was done by Vinko et al. in Ref. [39]. In this manner, any sum over k can be evaluated using a simple linear grid in the cubic BZ. However, our real-space unit cell has a larger volume than the primitive unit cell of FCC. This implies that the volume of our cubic BZ is smaller than the first BZ of FCC. At first glance, this seems to imply we have lost information on k points that lie outside of our cubic BZ but are located inside the first BZ of FCC. However, our employed periodic Bloch formalism ensures that this is not the case: the k-dependent quantities, such as E n,k , that lie in the first BZ of FCC, but outside of our cubic BZ are then shifted to a higher band, n. This effect will reveal itself explicitly in the band structure plots shown in Sec. III B.
The geometry of the defined crystal structure implies a radius of the touching spheres r s = 2.706 a 0 for all atoms. Due to the high level of symmetry, we found that r V (a) = r s = 2.706 a 0 for all atoms, and for all temperatures considered. Thereby, in Eq. (17), r C (a) = r s , which defines the 1s, 2s, and 2p orbitals as the core orbitals. Note that through the use of thermal occupation numbers with a well-defined chemical potential and temperature in xcrystal as was mentioned in Sec. II E, we assume the electrons to be in perfect thermal equilibrium for the conditions in the LCLS experiment, which is merely an approximation.

A. K-shell thresholds and ionization potential depression
The LCLS experiment described in Refs. [18,19] provides the data on the K-shell threshold energy. The K-shell threshold energy refers to the energy that is required for the formation of a K-shell (1s) hole. From the spectrally resolved Kα-emission spectrum in Ref. [18], they could deduce this quantity by observing at which XFEL-photon energy the emission peak started. Additionally, each separate emission peak observed indicates a specific charge state Q of an Al ion, which this peak refers to. Therefore, we will consider the K-shell threshold energy as a function of charge state Q. In addition to a comparison with this experiment, we will compare the result of XCRYSTAL with the results from the AA model and the two-step-HFS model, both shown in Ref. [29]. As such, we will first elaborate on how the K-shell threshold energy and the charge state Q are defined in these three models.
In the AA model, a thermal averaging is employed, so an integer value of Q corresponding to a specific electronic configuration cannot be obtained. Instead, an average chargē is obtained, with Z the charge of the neutral Al atom and N AA b (β, μ) the number of bound electrons. In Ref. [29] a muffin-tin potential is also employed, albeit with Wigner-Seitz spheres, and the flat potential provides a clear energetic distinction between bound orbitals and the continuum. Therefore, N AA b (β, μ) is calculated as the sum over the thermal occupations of orbitals with energies below the flat potential [see Eq. (22a)]. The K-shell threshold energy is calculated as the energy difference between the 1s orbital and the lowest-lying orbital above the flat potential. To be complete, there are two subtleties: (1) if the chemical potential μ is larger than the flat potential, the K-shell threshold energy is calculated as μ − ε 1s to make sure the orbital being excited into is vacant, and (2) at higher temperatures, the 3p orbital drops below the flat potential and a 1s electron can be excited into it. In this case, the K-shell threshold energy is ε 3p − ε 1s . This definition for the K-shell threshold energy is also used in the two-step-HFS method [29]. However, the two-step-HFS model is capable of specifying integer charge states, as it uses an AA calculation in its first step to define a thermal grand-canonical ensemble and proceeds in its second step by working with the most likely electronic configuration for a specific integer charge state Q. In contrast to this, but similar to the AA model, XCRYSTAL also works with a thermally averaged, and atom-specific, charge stateQ XCRY (a). However, our flat potential V 0 does not share the same physical interpretation as the flat potential used in Ref. [29], as delocalized states can be found below V 0 in XCRYSTAL. Therefore, we calculateQ XCRY , μ, a), calculated as the sum of the thermal occupations of the core orbitals (1s, 2s, and 2p). In summary, where V AA 0 is the flat potential in the AA model and p = (n, l, m, s) for a fixed a. Due to the translational symmetry in the crystal, we found the average chargeQ XCRY (a) for each ion to be the same, for a given temperature T . The Kα fluorescence energies employed in Refs. [18,19] to assign atomic charges to spectroscopic features are sensitive primarily to K-and L-shell occupation numbers. The choice made in Eq. (22b) is consistent with this property. As for the K-shell threshold energy, it refers to the energy that is required for the formation of a K-shell (1s) hole by exciting a 1s electron into the lowest-lying, vacant, delocalized valence orbital. In the convention employed here, the K-threshold corresponds to the energy difference between the 1s orbital and the lowest-lying orbital that has a thermal occupation of 0.5 and that is associated with a nonzero bandwidth. Note that this definition is equivalent to the one employed in Ref. [29].
In Fig. 3(a) we plot the K-shell threshold energy as a function of charge state Q calculated with XCRYSTAL and compare it to the experimental data [19], as well as to the theoretical models mentioned previously [29]. For both the AA model and the XCRYSTAL result, the average chargeQ for each method, as defined previously, is portrayed on the x axis. TheQ values shown correspond to a temperature range of 10-90 eV and 10-70 eV for the AA model and XCRYSTAL, respectively. The XCRYSTAL results were shifted by +20.92 eV, similar to what was done in Ref. [29]. This shift corresponds to the difference between the ionization potential calculated at T = 0 eV with XCRYSTAL (1538.68 eV) and the experimentally estimated binding energy (1559.60 eV) [82]. As shown in Fig. 3(a), XCRYSTAL is capable of reproducing the result for the K-shell threshold energy for all charge states considered.
The improvement of the XCRYSTAL result with respect to the AA result alone seemingly implies that this improvement is due to the incorporation of the full crystal structure in XCRYSTAL. However, the two-step-HFS model does not take 033061-8  [29], as well as with the experimental data from Ref. [19]. Calculations using an isolated Al ion are labeled as "Unscreened HFS." (b) K-shell threshold for Al as a function of the charge state, calculated using XCRYSTAL, the two-step-HFS result [29], and the average-atom result withQ defined as in XCRYSTAL [Eq. (22b)]. The experimental data from Ref. [19] is shown as well.
the crystal structure into account and yet shows excellent agreement with experiment as well, through its individual configuration optimization. To shed some light on how large the effects of incorporating a full crystal structure are, we may dismiss their inclusion in XCRYSTAL. This simply amounts to an AA calculation, where the average charge is calculated asQ XCRY , i.e., the thermal occupations of the 1s, 2s, and 2p orbitals subtracted from the nuclear charge of Al. We show this result in Fig. 3(b). The temperature range for the AA result shown in Fig. 3(b) is 10-80 eV. Comparing the AA result between Figs. 3(a) and 3(b), we conclude that the discrepancy between experiment and the AA method shown in Ref. [29] is not so much a limitation of the applicability of the AA model, but rather a consequence of the definition of Q AA used in Ref. [29].
To further strengthen this claim, we may considerQ AA and Q XCRY (a) for a single Al ion as a function of the electronic temperature, depicted in Fig. 4. Initially,Q AA increases with rising temperature as the bound orbitals are being partially thermally vacated. However, as the temperature increases, the previously unbound 3s and 3p orbitals fall below the flat potential V AA 0 [29] and count as bound states. This causes a sudden large contribution to the number of bound electrons, thereby decreasing the average chargeQ AA below what one would expect if such a sudden addition of new bound orbitals had not taken place. The discontinuities at T = 20 eV (Q AA ∼ 3.0) and T = 58 eV (Q AA ∼ 5.5) shown in Fig. 4 directly result in the bumps for the K-shell threshold energy in Fig. 3(a), seen at these same values ofQ. In contrast,Q XCRY does not exhibit these sudden drops. The increase in average charge arises solely from gradually thermally vacating the core orbitals. This analysis, along with the results shown in Fig. 3(b), leads us to the conclusion that both the incorporation of the entire crystal structure, as well as individual configuration optimization, amount to a fairly limited effect overall. The property shared by the AA model, the two-step-HFS approach and XCRYSTAL is the optimization of core orbitals, to which we attribute the success of the three models.
We note that the work performed by Vinko et al. [39] seemingly contradicts this conclusion. Vinko et al. show an excellent agreement with the experiment of Ref. [19]. They used a plane-wave DFT calculation with a frozen-core pseudopotential determined at a fixed configuration and obtained values for the K-shell threshold energy using a SCF scheme. The agreement with experiment was rationalized through their incorporation of a full three-dimensional electronic structure for the valence states and a lack of any spherical averaging. However, as both XCRYSTAL and the two-step-HFS model are able to adequately reproduce the experimental results of Ref. [19], which do incorporate both spherical and thermal averaging, we disagree with this proposed justification. The fact that Vinko et al. employed a frozen-core pseudopotential apparently contradicts our conclusion that core-orbital optimization is of vital importance. However, they do indirectly account for responsive core orbitals through their employment of a SCF scheme to calculate the K-shell threshold energy. Explicitly, they calculate the difference in total free energy between a system with and without a single K-shell hole, accompanied by an additional number of Lshell holes. Therefore, despite the 1s orbital being frozen in both configurations, their energies differ, resulting in an indirect response from the 1s core orbital to the K-shell threshold energy.
In order to illustrate the efficiency of our hybrid-basis approach, we compare numerical parameters of XCRYSTAL with the work done by Vinko et al. [39]. They considered high-energy bands with energies up to at least ten times the considered electronic temperature. Because of the high computational cost, they were incapable of making all temperature calculations with a single fixed N K (see Sec. II B); instead the number of bands ranged from ∼160 bands at T = 10 eV to ∼5000 bands at T = 100 eV. The computational expense in the latter case was reported to exceed ∼190 CPU days for a four-atom supercell containing 25 electrons. One should note that this cell did not contain all the electrons of a chargeneutral system containing four Al atoms, 4 × Z (Al) = 52, because of the employed pseudopotential formalism which used frozen, and not completely filled, core orbitals. In contrast, using XCRYSTAL we performed full all-electron calculations with 52 electrons. All calculations were consistently done with |K| max = 6.0 a −1 0 , corresponding to N K = 1647 bands. We found that this number of bands (1647 < ∼5000) was sufficient to obtain converged XCRYSTAL results. Every temperature run could be completed within 72 CPU days for a four-atom unit cell containing 52 electrons. Additional parallelization of XCRYSTAL further reduced the calculation time to ∼18 days with four CPU cores.
Finally, we present our result for the IPD of this system. For the 1s orbital, this is defined as the difference between the result from the HFS calculation for an unscreened (isolated ion) and the K-shell threshold energy [from Fig. 3(a)]. Defining the ionization potentials for the 2s and 2p orbitals in the analogous way as we did for our K-shell threshold, we can obtain IPD values for these orbitals as well. These values are shown in Fig. 5(a). The XCRYSTAL result was interpolated onto integer values of the charge state Q. Despite the ability of XCRYSTAL to calculate orbital-specific IPD values, Fig. 5(a) illustrates that the differences between the IPDs of different orbitals are minimal, in agreement with the observation of Son et al. [29]. Figure 5(b) shows the IPD values for the 2p orbital calculated with XCRYSTAL, comparing them with the two-step-HFS result [29], and the results from the modified-EK and the SP models, taken from Ref. [19]. We can see that XCRYSTAL reflects a similar trend as the two-step-HFS model, but lies closer to the IPD values determined from the modified-EK model as Q is increased.

B. Band structure at WDM temperatures
We proceed with calculations of the temperaturedependent band structure for Al at WDM conditions. We Comparison between the XCRYSTAL result for the 2p orbital, the two-step HFS model result [29], the result from the modified-EK model, and the result from the SP model, taken from Ref. [19]. The result from XCRYSTAL has been interpolated onto integer values of Q.
start with the discussion of Al at T = 0 eV. In Fig. 6 we present the first 10 energy bands of a nonzero bandwidth at T = 0 eV, along the path -X-W-L--K-X in the BZ. This path was chosen for the comparison of our predictions with Refs. [80,81]. The energy axis on the left is given in Rydbergs (Ry) and was shifted by +1.54 Ry (+20.98 eV), in order to adjust the energy at the point to the energy of the lowest delocalized state, E min,deloc . An unshifted energy axis (in eV) is presented on the right. As was previously mentioned in Sec. III, our choice of a cubic unit cell implies Band structure of aluminum at T = 0 eV along the path -X-W-L--K-X calculated using XCRYSTAL. The blue line represents the XCRYSTAL result, and the red line traces out the lines that would be obtained when using the conventional primitive FCC unit cell, to help guide the eye for comparison with Refs. [80,81]. The horizontal black line indicates the Fermi energy. The energy axis on the left is given in Rydbergs (Ry) and was shifted by +1.54 Ry (+20.98 eV), so that the energy at the point starts at the origin. An unshifted energy axis given in eV is presented on the right. that k-dependent quantities, such as the energies E n,k shown in Fig. 6, that lie in the first BZ of FCC, but outside of our cubic BZ, are shifted to a higher band, n. Therefore, the result from any XCRYSTAL calculation shows many more bands as compared to conventional methods. This is simply a consequence of having chosen a cubic unit cell. Both the blue lines and the red lines presented in Fig. 6 are the result from the same XCRYSTAL calculation. The blue lines show the result as calculated using XCRYSTAL, whereas the red line traces out the bands which would be obtained for the conventional primitive unit cell. One can see excellent agreement with the well-known results of Refs. [80,81]. The solid black line denotes our value for the Fermi energy ε F , calculated as the chemical potential μ at T = 0 eV, with a predicted value of −10.12 eV. Keeping the shift of E min,deloc (=−20.98 eV) in mind, the Fermi energy is calculated relative to the bottom of the conduction band, as ε F = μ − E min,deloc . It is found to be 0.8 Ry = 10.86 eV, whereas the experimental value is 11.7 eV [83]. The value for the constant potential V 0 at T = 0 eV is −18.34 eV, which is above the energy value of the lowestlying delocalized energy band at E min,deloc . This indicates that the bands in Fig. 6 with an (unshifted) energy below V 0 = −18.34 eV correspond to electronic states in which the electron is quasi-bound and tunnels between atomic sites. They are not localized enough to be considered part of the core orbitals [cf. Eq. (17)]. Between each symmetry point ( , X, W, etc.) in Fig. 6, there are 50 k points shown, which is the value we shall retain for the remaining band structure plots. We do not perform an entire XCRYSTAL run with such a fine k grid, instead taking the converged electronic density from the run with the parameters mentioned previously (Sec. II G), and performing a single XCRYSTAL matrix diagonalization in the hybrid basis.

033061-10
In Fig. 7 we show the first 30 energy bands with a nonzero bandwidth at temperatures of T = 25, 50, 75, and 100 eV. For these plots, we traced out the path -X-W-L--K-W-U-X to cover all lines of high symmetry in the BZ. In addition to the interesting observation that at T = 75 eV there appear to be no orbital states present directly above V 0 , we can perceive three general features in Fig. 7 as the electronic temperature rises: (1) all energy bands are lowered, (2) two band gaps start to form and progressively separate, and (3) the bandwidths become smaller. The conclusion is that with the increasing temperature, we observe a formation of quasi-atomic 3s and 3p lines, that was also reported in Ref. [39].
The physical mechanism behind this observation is the following: As the temperature rises, thermal excitations start to generate partial vacancies within the bound 1s, 2s, and 2p orbitals. This causes the nucleus to experience less screening, which in turn makes its potential V (r) more attractive, thereby dragging down all energies. In addition, it is apparent that with increasing electronic temperatures, the bands start to lose their width and their delocalized character, thereby exhibiting features of a more atomic nature.

IV. CONCLUSION AND OUTLOOK
In this work, we have developed an ab initio method for calculating quantum states of hot thermalized electrons confined to a cold ionic crystal lattice. It has been implemented into a new toolkit, XCRYSTAL. Using a mean-field HFS approach in combination with the Bloch formalism, we constructed a hybrid basis consisting of both plane waves and localized atomic orbitals, with respect to which we represented the electronic states in this type of transient WDM system. Allowing for an interwoven optimization between core and valence electrons, we accurately reproduced the experimental data obtained from the LCLS experiment [18,19] on Al plasmas in WDM conditions, in a highly efficient manner. Additionally, our model allowed for the calculation of band structures at various temperatures, T ∼ 0-100 eV. We concluded that the incorporation of the full crystal structure had only a minor effect on the results calculated for comparison with Refs. [18,19] and argued that the role of optimized core orbitals is vital in describing these types of systems. In addition, the band structure of Al was shown to be in good agreement with previous work at zero temperature [80,81]. The computationally efficient scheme in XCRYSTAL facilitated the calculation of the band structure to temperatures and densities typical for WDM conditions.
The new tool XCRYSTAL provides not only the calculation of bands and energy levels but gives full access to the electronic wave functions in a thermalized electron plasma at a certain temperature. This will enable us in the future to model atomic processes in a dense plasma environment, which is critical for a proper description of WDM formation. Until now, the respective cross sections and rates were adapted from isolated-atom models (see, e.g., Ref. [84]). In particular, XCRYSTAL enables access to the evaluation of the electron-impact-ionization cross section in solid-density plasmas, also measured experimentally [85]. It is the ultimate goal of XCRYSTAL to provide fast and accurate data of WDM properties for electronic Monte Carlo simulations of WDM, so as to contribute to the description of matter exposed to high-intensity x-ray pulses, in a similar manner the codes XATOM [71,72], XMDYN [72,86], and XMOLECULE [87,88] do.
The inclusion of molecular dynamics is conceptually possible as was mentioned in Sec II B. For application of this work in an astrophysical context, it is conceivable that the inclusion of the effects of an external magnetic field will be necessary [89]. It would be possible to perturbatively treat the effect of an external magnetic field through an interaction Hamiltonian in XCRYSTAL. For magnetic fields that are nonperturbatively strong, their inclusion would break some of the symmetries being exploited in XCRYSTAL. We believe that the present work will pave the way for a more accurate theoretical description of nonisothermal WDM, and will significantly advance the understanding of this extreme state of matter.

ACKNOWLEDGMENTS
We would like to thank Murali Krishna Ganesa Subramanian for his help with analytical calculations, as well as Zoltan Jurek, Malik Muhammad Abdullah, and Otfried Geffert for their aid in software development. The authors acknowledge funding by the Deutsches Elektronen-Synchrotron (DESY).

APPENDIX A: NORMALIZATION OF CORE BLOCH WAVES
In this appendix, we will derive the normalization of the core Bloch wave functions ψ n C ,k (r) as defined in Eq. (6), from the normalization of the total electronic wave function ϕ n,k (r). The wave functions ϕ n,k (r) are Bloch functions [see Eq. (1)]. These wave functions are orthonormal within the entire crystal volume, ϕ n,k |ϕ n ,k = δ nn δ k,k . Therefore, we have where R denotes a lattice translation vector. We used (1/N cell ) R e −i(k−k )·R = δ k,k and V = N cell . This implies that we normalize the μ n,k (r) within the unit cell volume as For core electrons, we know a priori that our atomic orbital wave functions obey φ n C |φ n C = δ n C n C . However, for consistency, we should also normalize μ n C ,k (r) to , with respect to a unit cell volume , just as was the case for μ n,k (r). This is why the factor factor N in Eq. (6) was split up as N = N a N b . Imposing the preferred normalization, we find revealing that N b = √ . Furthermore, the μ n C ,k (r) are orthogonal to each other, but only for k = k : which is equal to δ n C ,n C only if k = k . In the second line, we used the fact that the n C label on r n C depends only on which atom we are considering, and unless we are considering the same atom, the integrand will be zero for all values of r within the unit cell because we are assuming nonoverlapping spheres. Let us now derive N a for a proper normalization of ψ n C ,k (r): ψ n C ,k (r) = N a e ik·r μ n C ,k (r) and ψ n C ,k ψ n C ,k = δ n C n C δ k,k , 033061-12 which yields ψ n C ,k ψ n C ,k = V d 3 rψ * n C ,k (r)ψ n C ,k (r) = N 2 a V d 3 re −ik·r μ * n C ,k (r)e ik ·r μ n C ,k (r) = N 2 a R d 3 re −ik·(r + R) μ * n C ,k (r + R)e ik ·(r + R) μ n C ,k (r + R) = N 2 a R e −i(k−k )·R d 3 re −i(k−k )·r μ * n C ,k (r)μ n C ,k (r) = N cell N 2 a δ k,k d 3 re −i(k−k )·r μ * n C ,k (r)μ n C ,k (r) = N cell N 2 a δ n C n C δ k,k , showing that N a = 1/ √ V . To summarize, the total wave function has the Bloch-wave form: ϕ n,k (r) = 1 √ V e ik·r μ n,k (r), satisfying ϕ n,k |ϕ n ,k = δ nn δ k,k , and we use the following ansatz with a hybrid basis: Here the core wave functions ψ n C ,k (r) are given by ψ n C ,k (r) = 1 √ V e ik·r μ n C ,k (r), satisfying ψ n C ,k ψ n C ,k = δ n C n C δ k,k , which is also represented as a Bloch wave, with the periodic function μ n C ,k (r) = √ I e −ik·(r−R I,n C ) φ n C r − R I,n C .
Note that μ n,k and μ n C ,k are normalized to with respect to the unit cell .

APPENDIX B: ACCELERATING SCF CONVERGENCE
In this appendix, we will elaborate on a new scheme for improving SCF convergence, which proved very successful for our computational framework.
The convergence of the SCF method is generally rather slow. Therefore, much effort has been invested into methods improving it [81]. A simple example of such a method is linear mixing, described in Ref. [81]. In an attempt to steer the convergence, one uses the information from the previous iterations, i.e., instead of using the potential V i (r) in the ith SCF iteration, one uses V M i (r) = αV i (r) + (1 − α)V M i−1 (r) with a free parameter α ∈ [0, 1]. We distinguish here between the potential V i (r) obtained from the electronic density using Eq. (19), and the "mixed" potential V M i (r) obtained from performing the linear mixing. A different quantity, such as the electron density ρ, may also be used in the mixing [81]. Building on this notion of steering convergence, we have developed a method that adjusts the α parameter per iteration, using the information from all previous iterations. In this way, we are working with α → α i .
Using the information from previous iterations, we impose for the ith SCF iteration that with weights {w m }, which are normalized to 1, i.e., i m=1 w m = 1. It is intuitively understandable that we expect the weight w i associated with the ith iteration to be larger if the error associated with that iteration is smaller. With that in mind, we define c i = ||e i || −1 , for some suitably defined error vector e i . It is related to w i through the normalization condition i m=1 w m = 1: Note that as more SCF iterations are being considered, a specific c m will remain unchanged, whereas k m , and therefore w m , will be updated for every iteration. We can easily combine From this, we can identify α i = k i c i = w i . We can also easily show that k i /k i−1 = 1 − α i , by using (k i−1 ) −1 + c i = (k i ) −1 .
Multiplying both sides of this equation by k i , we get α i + (1 − α i ) = 1, as expected. In the implementation of XCRYSTAL, we defined the error vector to be e i = |V i (r) − V i−1 (r)|/V i−1 (r).
For the first iteration, no mixing is done as the error vector cannot be defined yet. To demonstrate how the new adaptive linear mixing scheme works for our particular framework, we compare the performance of the adaptive linear mixing to static linear mixing, with two different mixing parameters α, in Fig. 8. It illustrates not only that the adaptive linear mixing is useful for speeding up the convergence, but that its application is necessary for high-temperature cases to achieve convergence. The static linear mixing shown in Fig. 8 exhibits an exponential behavior with increasing temperatures and eventually fails to converge within a convergence criterion of 10 −6 % for T > 35 eV (α = 0.6) and T > 15 eV (α = 1.0). We used a relative difference of the total, Fermi-Dirac-weighted energy of the system n,k E n,knn,k between iterations as the SCF-convergence parameter, achieving convergence if this relative difference was smaller than 10 −6 %.