*Article* **Application of the Non-Equilibrium Statistical Operator Method to the Dynamical Conductivity of Metallic and Classical Plasmas**

#### **Mikhail Veysman 1,\*, Gerd Röpke 2,3 and Heidi Reinholz 2,4**


Received: 22 March 2019; Accepted: 23 April 2019; Published: 7 May 2019

**Abstract:** The fruitfulness of the method of a non-equilibrium statistical operator (NSO) and generalized linear response theory is demonstrated calculating the permittivity, dynamical conductivity, absorption coefficient, and dynamical collision frequency of plasmas in the degenerate, metallic state as well as classical plasmas. A wide range of plasma parameters is considered, and a wide range of frequencies of laser radiation acting on such plasmas is treated. New analytical expressions for the plasma response are obtained by this method, and several limiting cases are discussed.

**Keywords:** non-equilibrium statistical operator; linear response theory; permittivity, dynamical conductivity, absorption coefficient, dynamical collision frequency; ordered lattice, disordered lattice; Umklapp process; interband transitions

#### **1. Introduction**

The non-equilibrium quantum statistical operator (NSO) was proposed by D.N. Zubarev [1], and its 100-year anniversary was celebrated recently. The NSO method is an important step in working out a unified, general approach to nonequilibrium phenomena such as transport and relaxation processes. Different approaches to describe nonequilibrium processes, such as kinetic theory, linear response theory, and quantum master equations, are obtained within a very general approach, after specifying the relevant degrees of freedom that characterize the non-equilibrium state of the system. Recent reviews of the NSO method and its applications can be found, e.g., in [2–6].

The NSO method has been successfully applied to many problems related to the transport and optical properties of charged particle systems, such as warm dense matter (WDM). Reviews of the calculations of dynamical response of dense plasmas are found, e.g., in [7,8]. A main advantage of the approach is that special approaches, valid in limiting cases, can be generalized to describe more complex situations. For instance, kinetic theory has been worked out for low-density systems where the single-particle distribution function is relevant to describing the non-equilibrium state of the system, and correlations can be neglected. In contrast, linear response theory has been worked out to describe systems at arbitrary densities, but near thermodynamical equilibrium. A unified theory that covers both limiting cases has been worked out using the NSO method, see [9]. In particular, the correct frequency dependence of the response functions has been found, in contrast to conventional kinetic theory. A further advantage of the approach is that

non-equilibrium properties of quantum as well as classical plasmas can be considered consistently on the same footing. Altogether, the NSO method allows one to study the transport and optical properties of plasmas at arbitrary densities and temperatures, including the region of WDM where correlations are strong, in particular between ions, and where the electrons become degenerate. The response to classical Maxwell fields, represented by a time-dependent external field in the Hamiltonian, can be investigated in a wide range of frequencies of the external field. This refers also to the response of WDM that is irradiated by laser fields, in particular the frequency-dependent absorption coefficient. To include a spontaneous emission of photons, a quantum description of the Maxwell fields is necessary.

In this paper, we present briefly the NSO method and its application to the permittivity of metallic plasmas, which is described by a Hamiltonian accounting for electron–phonon interactions and Umklapp processes, as well as the application to the permittivity of plasmas without a long-range order of ions, which is described by a Hamiltonian accounting for electron–ion and electron–electron interactions. We continue our previous investigations to work out a systematic quantum statistical approach to the response properties of WDM [10–14] with an application to aluminum, considering additional processes of interaction in the strongly coupled Coulomb system. We derive analytical expressions for the dynamical collision frequency and related quantities, in particular the absorption coefficient. The behavior of these response properties, such as temperature dependence and frequency dependence, is discussed in special limiting cases.

These results are of interest for the investigation of material at extremal conditions, e.g., high-particle and high-energy densities. If a solid metallic target is irradiated by powerful laser pulses, the electron component of such a target undergoes modifications of its properties, changing from the state of a metallic "plasma" to a state of a classical non-degenerate plasma. This transition of the electron system from a degenerate state to a classical behavior is consistently described in our approach. In addition, the ion system may change from a solid state, where the ion positions are strongly correlated forming a lattice, to a liquid or a plasma state, where the long-range order of ion positions is destructed. In the present work, we neglect the short-range order of the ion system in the liquid or plasma state, so that the electron–ion interaction is considered as independent scattering at the individual ions. For the electron–ion interaction in the solid state of matter, the coherent part of multiple scattering by the ionic lattice leads to the formation of electron band (Bloch) states. The deviation of ions from the lattice position is described as phonon excitation, and the electron–phonon interaction together with Umklapp processes is considered in this work as responsible for the dynamical conductivity of the solid metal, in addition to electron–electron collisions. Such transitions from collective electron–phonon interactions in a solid to individual electron–ion interactions in the disordered ion configuration directly opens up the question about the switching between the respective Hamiltonians. A more general approach is possible using the concept of the dynamical structure factor which reflects not only the configuration of the ions but also the dynamical behavior, including collective excitations, of the ion system. This problem, however, requires a separate investigation and is not considered in the present work.

With respect to the application to aluminum plasmas, recent experiments in the WDM region [15,16] to measure the dynamical conductivity are also treated by density-functional theory (DFT) for electrons combined with molecular dynamics (MD) calculations for the ions [17]. These calculations have the advantage that optimal single-electron orbitals are calculated, and these orbitals reflect the electron structure of the aluminum ion, improving our effective electron–ion (*e* − *i*) interaction potential. In addition, the static ionic structure factor is calculated so that a short-range order in the liquid or plasma phase is implemented in the calculation of the dynamical conductivity. However, the choice of the density functional for the correlation energy remains an open problem of these DFT-MD calculations. A shortcoming of the DFT-MD approach is that electron–electron (*e* − *e*) collisions are not properly included

in this mean-field theory, in contrast to our generalized linear response approach, where the contribution of *e* − *e* collisions is taken into account, see also [18].

A systematic improvement of our approach is the use of optimal single-electron orbits and *e* − *i* interaction potentials. In addition, the dynamical ion structure factor that causes not only multiple scattering of the electrons, but also the excitation of the ion system, in particular phonons, should be considered in a systematic approach. These improvements in our calculations that establish a closer connection to DFT-MD simulations are a subject of future work.

#### **2. A Brief Description of the Method**

The plasma permittivity *ε*(*ω*) of an electromagnetic field with the frequency *ω* is directly related to the dynamical conductivity *σ*(*ω*) according to *ε*(*ω*) = 1 + *iσ*(*ω*)/(*ε*0*ω*) or, using Gauss units instead of SI units, as follows:

$$
\varepsilon(\omega) = 1 + 4\pi i \sigma(\omega) / \omega. \tag{1}
$$

Here and below we consider isotropic media in the case of weak spacial dispersion or the long-wavelength (with respect to plasma perturbations) limit. In this case, the permittivity does not depend on the wave vector of perturbations. The longitudinal permittivity coincides with the transverse one: *<sup>ε</sup>* = *<sup>ε</sup>*<sup>⊥</sup> = *<sup>ε</sup>* [19].

For the calculation of *σ*(*ω*), one should determine the (quantum) statistically averaged electric current, arising as the system's response to external electromagnetic fields. Such calculation can be done using the NSO *ρ*ˆ = *ρ*ˆrel + *ρ*ˆirrel, which can be constructed as a sum of the so-called *relevant statistical operator ρ*ˆrel describing the quasi-equilibrium and the *irrelevant statistical operator ρ*ˆirrel, which represents the nonequilibrium contribution.

The *relevant statistical operator* is introduced as a generalized Gibbs ensemble, derived from the principle of the maximum of entropy:

$$\begin{aligned} \phi\_{\text{rel}}(t) &= Z\_{\text{rel}}(t)^{-1} \exp\left[ -\beta(\hat{H} - \mu \hat{N}) + \sum\_{n} F\_{n}(t)\hat{\mathcal{B}}\_{n} \right] \\ Z\_{\text{rel}}(t) &= \text{Tr}\left[ -\beta(\hat{H} - \mu \hat{N}) + \sum\_{n} F\_{n}(t)\hat{\mathcal{B}}\_{n} \right] \end{aligned} \tag{2}$$

where *H*ˆ is the Hamiltonian of the system. *N*ˆ = ∑*<sup>p</sup> a*ˆ† *pa*ˆ*<sup>p</sup>* is the particle number operator, where *a*ˆ† *<sup>p</sup>* and *a*ˆ*<sup>p</sup>* are the creation and annihilation operators of the electrons in state *p*. For the chosen set of relevant observables {*Bn*}, *n* = 1... N , we request

$$\text{Tr}\left\{\mathcal{B}\_n\rho(t)\right\} = \langle\mathcal{B}\_n\rangle^t = \text{Tr}\left\{\mathcal{B}\_n\rho\_{\text{rel}}(t)\right\},\tag{3}$$

meaning that the observed statistical averages ...*<sup>t</sup>* at time *<sup>t</sup>* are correctly reproduced by the averages with the relevant statistical operator *ρ*ˆrel(*t*). The Lagrange parameters *Fn*(*t*) are then determined by the set of requested expressions (3) as response parameters. Similar conditions on *N*<sup>ˆ</sup> and *H*<sup>ˆ</sup> determine *μ* as chemical potential and *β* = 1/*Te*, respectively, where *Te* is the electron temperature. (We give the temperature in energy units, *k*<sup>B</sup> = 1, instead of SI units where *β* = 1/*k*B*Te*.)

If N = 0 (i.e., an empty set of relevant observables), the operator *ρ*ˆrel (2) is identical with the statistical operator for the grand canonical ensemble [3]. In that case, the well-known Kubo formula [20] for the electrical conductivity follows immediately [4,5].

The choice of relevant observables can be arbitrary. Different choices lead to the same results if non-perturbative approaches such as infinite diagram summations or numerical calculations of correlation functions are used [5,8,21]. However, the account of only a finite number of terms within a perturbation expansion can lead to even divergent results. Therefore, it is necessary to choose the set of relevant

observables in a way that the relevant statistical operator [5] already ensures a close approximation of the considered system.

According to [9] (see also [4,8,21]), the NSO *ρ*ˆ(*t*) is determined by the dynamical evolution of the system with Hamiltonian *H*ˆtot = *H*ˆ + *H*ˆ ext(*t*):

$$\rho(t) = \lim\_{\eta \to +0} \eta \int\_{-\infty}^{t} dt' e^{-\eta(t-t')} \hat{\mathcal{U}}(t, t') \rho\_{\text{rel}}(t') \hat{\mathcal{U}}^{\dagger}(t, t') \tag{4}$$

where *U*ˆ (*t*, *t* ) is the time evolution operator, which solves the equation

$$i\hbar \partial\_l \hat{\mathcal{U}}(t, t') = \hat{H}\_{\text{tot}} \hat{\mathcal{U}}(t, t') \,, \quad \hat{\mathcal{U}}(t', t') = 1 \tag{5}$$

with *H*ˆ ext(*t*) being the Hamiltonian of the external perturbation. Due to Equation (4), correlations from the initial state are further built up, which is determined by the relevant statistical operator *ρ*ˆrel(*t*) (2).

For a high frequency electromagnetic field with an electric field strength *E*(*t*) acting on matter, the external perturbation *H*ˆ ext(*t*) can be written in dipole approximation as

$$
\hat{H}\_{\text{ext}}(t) = -e\hat{\mathbf{R}} \cdot \mathbf{E}(t), \ \hat{\mathbf{R}} = \sum\_{a} \mathfrak{k}\_{a}, \ \dot{\mathbf{R}} = \hat{\mathbf{P}}\_{1\Sigma}/m\tag{6}
$$

where *m* is the electron mass, and *e* is the electron charge. *P*ˆ <sup>1</sup><sup>Σ</sup> = ∑*νP*ˆ 1,*<sup>ν</sup>* is the operator of the total momentum of electrons, which is the sum of momentums of electrons from different energy zones *ν*. It coincides with the first moment of the density matrix, see Equation (11) below.

In *linear response theory*, which we assume to be applicable, an expansion of the relevant *ρ*ˆrel(*t*), Equation (2), and the irrelevant *ρ*ˆirrel(*t*) = *ρ*ˆ(*t*) − *ρ*ˆrel(*t*) statistical operators with respect to the external perturbation and the response parameters *Fn*(*t*) are considered, see [5,9]. Together with Equation (3) and using the Kubo identity, this gives rise to the following system of equations:

$$
\langle \delta \hat{B}\_n \rangle = \sum\_{m} (\hat{B}\_n ; \delta \hat{B}\_m) F\_m \tag{7}
$$

$$\sum\_{m} \left[ -i\omega \left\{ (\boldsymbol{\mathcal{B}}\_{n}; \boldsymbol{\mathcal{B}}\_{m}) + \langle \boldsymbol{\hat{\mathcal{B}}}\_{m}; \boldsymbol{\mathcal{A}} \boldsymbol{\mathcal{B}}\_{m} \rangle\_{z} \right\} + (\boldsymbol{\mathcal{B}}\_{n}; \boldsymbol{\hat{\mathcal{B}}}\_{m}) + \langle \boldsymbol{\hat{\mathcal{B}}}\_{n}; \boldsymbol{\hat{\mathcal{B}}}\_{m} \rangle\_{z} \right] \boldsymbol{\mathcal{F}}\_{m} = \boldsymbol{\beta} \frac{\boldsymbol{\mathcal{E}}}{m} \left\{ (\boldsymbol{\mathcal{B}}\_{n}; \boldsymbol{\mathcal{P}}\_{\mathrm{1\Sigma}}) + \langle \boldsymbol{\hat{\mathcal{B}}}\_{m}; \boldsymbol{\mathcal{P}}\_{\mathrm{1\Sigma}} \rangle\_{z} \right\} \boldsymbol{\mathcal{E}} \tag{8}$$

where *<sup>z</sup>* <sup>=</sup> *<sup>ω</sup>* <sup>+</sup> *<sup>i</sup>η*; *<sup>δ</sup>B*ˆ*<sup>n</sup>* <sup>=</sup> *<sup>B</sup>*ˆ*<sup>n</sup>* − *B*ˆ*n*0, and *B*ˆ*n*<sup>0</sup> is the statistical average of *<sup>B</sup>*ˆ*<sup>n</sup>* with the *equilibrium statistical operator* of the grand canonical ensemble *ρ*ˆ0 = Z−<sup>1</sup> <sup>0</sup> exp[−(*H*<sup>ˆ</sup> <sup>−</sup> *<sup>μ</sup>N*<sup>ˆ</sup> )/*Te*] with <sup>Z</sup><sup>0</sup> <sup>=</sup> Tr{*e*−(*H*<sup>ˆ</sup> <sup>−</sup>*μN*<sup>ˆ</sup> )/*Te*}.

In Equations (7) and (8), expressions such as (*A*ˆ; *B*ˆ) and 0 *A*ˆ; *B*ˆ 1 *<sup>z</sup>* denote Kubo scalar products of operators *A*ˆ and *B*ˆ and the Laplace transform of the Kubo scalar product of these operators, respectively. The latter are called *equilibrium correlation functions*. They are defined by expressions

$$\hat{A}\left(\hat{A};\mathcal{B}\right) \quad = \int\_0^{\beta} d\tau \,\mathrm{Tr}\left\{\hat{A}(-i\hbar\tau)\mathcal{B}^+\rho\_0\right\} \tag{9}$$

$$\langle \langle \hat{A}; \mathcal{B} \rangle\_z \rangle\_z = \int\_0^\infty dt e^{jzt} \left( \hat{A}(t); \mathcal{B} \right) \tag{10}$$

where the operator *A*ˆ is taken in Heisenberg representation *A*ˆ(*t*) = *eiHt* <sup>ˆ</sup> /¯*hAe*ˆ <sup>−</sup>*iHt* <sup>ˆ</sup> .

For further treatment of Equations (7) and (8), it is necessary to define the set of relevant observables *Bn*. For the description of the permittivity of plasmas, it is convenient to choose the moments of the density matrix as the set of relevant observables:

$$\mathcal{B}\_{\mathbf{n}} = \mathcal{P}\_{\mathbf{n}} = \sum\_{p\_{\mathbf{v}}} \hbar \, p\_{\mathbf{v}} \hbar\_{p, \mathbf{v}} (E\_{p, \mathbf{v}} / T\_{\mathbf{c}})^{(L-1)/2}, \quad \mathbf{n} \equiv \{L, \mathbf{v}\} \tag{11}$$

where *n*ˆ *<sup>p</sup>*,*<sup>ν</sup>* = *a*ˆ<sup>+</sup> *<sup>p</sup>*,*νa*ˆ*p*,*ν*, *p<sup>ν</sup>* = *m ∂Ep*,*ν*/*∂p*, and

$$E\_{\mathcal{V},\mathcal{V}} = \mathfrak{p}^2 / (2\, m \, m\_{\mathcal{V}}) + E\_{\mathcal{O},\mathcal{V}} \tag{12}$$

where *m<sup>ν</sup>* and *E*0,*<sup>ν</sup>* are the effective electron mass (normalized to electron mass *m*) and the energy of the bottom of the *ν*-th zone, respectively. Here, for generality, we consider moments of the density matrix stipulated by different powers of electrons momentum *p* (labeled by the index "*L*") and by different energy zones (labeled by the index "*ν*"). The index "n" of observables comprises both *L* and *ν*; *a*ˆ<sup>+</sup> *<sup>p</sup>*,*<sup>ν</sup>* and *a*ˆ*p*,*<sup>ν</sup>* are creation and annihilation operators for single-electron states with momentum *p* in the *ν*-th band.

The equilibrium correlation functions occurring in Equations (7) and (8) will be evaluated below. Here we only mention that (*B*ˆ*n*; <sup>ˆ</sup> *B*˙ *<sup>m</sup>*) = 0 because the commutator [*B*ˆ*n*, *B*ˆ*m*] vanishes. In the lowest order of interaction, which is proportional to *<sup>e</sup>*2, one can show [22] that the terms in (8) containing only one operator <sup>ˆ</sup> *B*˙ *n* can be omitted in comparison to the leading order. In an isotropic system, we take the electrical field as well as the current density in the *z* direction so that only the absolute value is of relevance. Equations (7) and (8) are rewritten as the following system of equations for the dimensionless response parameters F<sup>m</sup> in terms of the dimensionless correlation functions Nnm and Cnm:

$$
\langle \hat{P}\_{\rm n} \rangle = \frac{e n E}{\omega\_{\rm n.u.}} \sum\_{\rm m} \mathfrak{N}\_{\rm n.m.} \mathcal{F}\_{\rm m} \tag{13}
$$

$$\sum\_{\mathbf{m}} \left[ \mathfrak{C}\_{\text{nm}} - i\omega \mathfrak{N}\_{\text{nm}} \right] \mathcal{F}\_{\text{m}} = \sum\_{\mu} \mathfrak{N}\_{\text{n}\{1\mu\}} \tag{14}$$

where <sup>N</sup>n{1*μ*} has the second index with *<sup>L</sup>* = 1 and band index *<sup>μ</sup>*,

$$\begin{split} \mathfrak{R}\_{\text{nm}} &= \frac{(\hat{\mathbf{P}}\_{\text{n}}; \hat{\mathbf{P}}\_{\text{m}})}{mnT\_{\text{c}}}, \ \mathfrak{C}\_{\text{nm}}(\omega) = \frac{\langle \hat{\mathbf{P}}\_{\text{n}}; \hat{\mathbf{P}}\_{\text{m}} \rangle\_{\omega + i\eta}}{mnT\_{\text{c}}\omega\_{\text{a},\text{u}}}, \\ F\_{\text{m}} &= \mathcal{F}\_{\text{m}} \frac{\varepsilon E}{mT\_{\text{c}}}. \end{split} \tag{15}$$

*n* is the particle number density of electrons in the conduction band (free electron density). *ω*a.u. is the atomic unit of the frequency, so that *<sup>h</sup>*¯ *<sup>ω</sup>*a.u. <sup>=</sup> *EH* <sup>=</sup> *me*4/¯*h*<sup>2</sup> <sup>≈</sup> 27.2 eV is the Hartree energy unit, and *ω*˜ = *ω*/*ω*a.u. is the dimensionless frequency. The indexes m, n = {*L*, *ν*} contain *L*, the power of the momentum, and *ν*, the number of the energy band.

The electrical current ˆ*J* <sup>=</sup> *<sup>σ</sup> <sup>E</sup>* can be calculated using the relation <sup>ˆ</sup>*<sup>J</sup>* <sup>=</sup> *eP*<sup>ˆ</sup> <sup>1</sup>Σ/*m*. Inserting expression (13), we derive the permittivity (1) as a Drude-like formula

$$\lim\_{k \to 0} \varepsilon(k, \omega) = \varepsilon(\omega) = 1 - \frac{\omega\_{\text{pl}}^2}{\omega \left[\omega + i\nu(\omega)\right]'} \tag{16}$$

with *ω*<sup>2</sup> pl = <sup>4</sup>*πne*2/*<sup>m</sup>* and an *effective collision frequency <sup>ν</sup>*(*ω*), which can be expressed in terms of dimensionless response parameters and correlation functions Nnm and Cnm, as defined in (15), according to

$$\boldsymbol{\nu}(\omega) = \boldsymbol{\nu}\_1(\omega)\boldsymbol{r}\_\omega(\omega) \tag{17}$$

$$\nu\_1(\omega) = \omega\_{\text{a.u.}} \frac{\mathfrak{C}\_{\text{11}}}{\mathfrak{N}\_{\text{11}}},\\r\_{\omega}(\omega) = \frac{\mathfrak{N}\_{\text{11}}}{\mathfrak{C}\_{\text{11}}} \frac{1 + i\omega^\* \sum\_{\mathbf{m},\nu} \mathfrak{N}\_{\{1,\nu\}\text{m}} \mathcal{F}\_{\text{m}}}{\sum\_{\mathbf{m},\nu} \mathfrak{N}\_{\{1,\nu\}\text{m}} \mathcal{F}\_{\text{m}}} \tag{18}$$

where the index "1" means n = 1 = {1, 1}, i.e., *L* = 1 in Equation (11), and *ν* = 1 is the conduction band. *ν*1(*ω*) is the complex *effective collision frequency* of electrons for Drude-like transitions (within the conduction band) and in single-moment approximation, while *rω* (18) is the so-called *renormalization factor*, which takes into account the influence of higher moments of the density operator [8,12,22] and, in the case considered here, the influence of transitions between different bands. Inserting the solution of the system of Equations (14) for F<sup>m</sup> into Expressions (18), we obtain the complex *effective collision frequency* in terms of the correlation functions Nnm and Cnm.

In our further considerations we look at the two following cases:

(a) *L* = 1, *ν* = 1, 2 ..., i.e., the case of different energy bands, but only the 1st moment of the density operator. Then <sup>n</sup> = {*L*, *<sup>ν</sup>*} = {1, *<sup>ν</sup>*} = *<sup>ν</sup>*. Inserting the definition of moments (13) for *<sup>P</sup>*{1,*ν*} into Expression (15) for N*νμ*, using the Kubo identity and expressing the electron momentum *p* via *r*˙, one can show [3] that

$$\mathfrak{N}\_{\nu\mu} = \delta\_{\nu,\mu} n\_{\nu} / n \tag{19}$$

where *nν* is the number of electrons in the *ν*-th band.

(b) *L* = 1, 2 ... , *ν* = 1, i.e., the case of a single conduction band and different moments of the density operator. Then n = {*L*, *ν*} = {*L*, 1} = *L*.

Expressions for the correlation functions N*lm*, *l*, *m* ≥ 1 can be found elsewhere (see, e.g., [3,18]):

$$\mathfrak{N}\_{lm} = \frac{\Gamma[(l+m+3)/2]}{\Gamma(5/2)} \frac{I\_{(l+m-1)/2}(\mathfrak{e}\_{\mu})}{I\_{1/2}(\mathfrak{e}\_{\mu})}, \ l, m \ge 1 \tag{20}$$

where *μ* = *μ*/*Te* is dimensionless chemical potential,

$$
\epsilon\_{\mu} = X\_{1/2} \left( 2 \epsilon\_{\text{F}}^{3/2} / 3 \right),
\tag{21}
$$

with the Fermi integrals *Iν*(*y*) = Γ(*ν* + 1) <sup>∞</sup> <sup>0</sup> *<sup>x</sup>ν*[*ex*−*<sup>y</sup>* <sup>+</sup> <sup>1</sup>] <sup>−</sup>1*dx*; <sup>F</sup> = *EF*/*Te* = Θ−1, where Θ is the degeneracy parameter, and *E*<sup>F</sup> = *h*¯ 2/(2*m*)(3*π*2*n*)2/3 is the Fermi energy; the dimensionless chemical potential *μ* in (21) is expressed via the inverse Fermi integral *X*1/2(*x*), which refers to the Fermi integral *I*1/2(*x*). Particularly, from (20) one has

$$\mathfrak{N}\_{11} = 1, \quad \mathfrak{N}\_{31} = \frac{5}{2} \frac{I\_{3/2}(\mathfrak{e}\_{\mu})}{I\_{1/2}(\mathfrak{e}\_{\mu})}, \quad \mathfrak{N}\_{33} = \frac{35}{4} \frac{I\_{5/2}(\mathfrak{e}\_{\mu})}{I\_{1/2}(\mathfrak{e}\_{\mu})}. \tag{22}$$

In the non-degenerate case, *Iν*(*μ*) = *e<sup>μ</sup>* for all *ν* and N<sup>31</sup> = 5/2, N<sup>33</sup> = 35/4, see [9].

In previous works [9,22], correlation functions with only the first and third moments of the density matrix (11) in the sum (18) were considered. It was shown that this leads to an accuracy of a few % for the calculation of the renormalization factor. Therefore, restricting our approximation to only two bands

or two moments of the density matrix in Cases (a) and (b), respectively, and using the solution of Equation (14), the *renormalization factor* can be written in terms of respective correlation functions:

$$\begin{split} r(\omega) &= \frac{1}{\mathfrak{C}\_{11}} \frac{1 + i\tilde{\omega}Q\_{\omega}}{Q\_{\omega}}\\ Q\_{\omega} &= \frac{\mathfrak{A}\_{22} - 2\mathfrak{A}\_{22}\mathfrak{A}\_{21} + \mathfrak{A}\_{22}^{2}\mathfrak{A}\_{11}}{\mathfrak{A}\_{11}\mathfrak{A}\_{22} - \mathfrak{A}\_{21}^{2}} \quad \text{for Case (a)}\\ Q\_{\omega} &= \frac{\mathfrak{A}\_{33} - 2\mathfrak{A}\_{31}\mathfrak{A}\_{31} + \mathfrak{A}\_{31}^{2}\mathfrak{A}\_{11}}{\mathfrak{A}\_{11}\mathfrak{A}\_{33} - \mathfrak{A}\_{31}^{2}} \quad \text{for Case (b)}\\ \mathfrak{A}\_{lm} &= \mathfrak{C}\_{lm} - i\tilde{\omega}\mathfrak{A}\_{lm} \; l, m \ge 1. \end{split} \tag{23}$$

Here, Expressions (19) and (22), respectively, have been used.

According to the definitions (10), (11) and (15), the correlation functions C*nm*(*ω*) can be expressed in terms of correlation functions of electron creation and annihilation operators as

$$\mathcal{C}\_{\text{nm}}(\omega) = \frac{-\beta}{\text{mm}\omega\_{\text{a.u.}}} \sum\_{p\_\nu, p\_\mu} p\_{\nu, z} p\_{\mu, z} \times \left(\mathbb{E}\_{p, \nu} / T\_e\right)^{\frac{1}{2} - 1} \left(\mathbb{E}\_{p, \mu} / T\_e\right)^{\frac{M - 1}{2}} \langle [\![\![\![\!n\_{\nu}\!] \!p\_{\nu}\!] \!], [\![\![\!n\_{\nu}\!] \!p\_{\mu}\!] \!] \rangle\_{\omega + i\eta} \tag{24}$$

where n = {*L*, *ν*}, and m = {*M*, *μ*}. The further calculation of the correlation functions Cnm requires the specific Hamiltonian *H*ˆ .

In the general case, the four-particle correlation function of creation and annihilation operators, arising after respective elementary transformations of the correlation function [*H*<sup>ˆ</sup> , *<sup>n</sup>*<sup>ˆ</sup> *<sup>p</sup>*,*ν*]; [*H*<sup>ˆ</sup> , *<sup>n</sup>*<sup>ˆ</sup> *<sup>p</sup>*,*μ*]*ω*+*i<sup>η</sup>* with known Hamiltonian *H*ˆ , can be expressed via thermodynamic Green functions [3,8]. Green function techniques allow one, in principle, to take into account all orders of interactions via the summation of respective Feynman diagrams [3,8,22]. Below we consider the first Born approximation, which follows at the lowest order with respect to the interaction either directly from the definition of the correlation functions in (24), or from the four-particle Green functions expressed as a product of single-particle Green functions, which are related to the correlation functions. As a result, we obtain concise and simple analytic results. It shall be noted that it can be reasonable to calculate the renormalization factor (18) in the Born or screened Born (see below) approximation and take into account strong-coupling effects only in the calculation of the collision frequency *ν*1(*ω*) (17), (18), see [18].

Below we look at two different cases:

(A) We consider individual interactions of electrons with each other and with randomly distributed ions. This is the case for ion temperatures higher than the melting one, where the long-range order of the ion lattice is destructed. The short-range order is described by the ionic structure factor. As an example, we can consider a metallic solid sample irradiated by intense short-pulse laser beams. The ion lattice disappears after a heating process longer than the characteristic melting time *τm*, which is of the order of the time between ion collisions. This time is given by the interatomic distance *ra* <sup>∼</sup> (4*πnat*/3)−1/3 divided by the sound velocity *vs* <sup>∼</sup> <sup>√</sup>*ZTe*/*mat*, where *nat*, *mat* are the concentration and the mass of the heavy particles, respectively, so that

$$
\tau\_m = \frac{r\_a}{\upsilon\_s} = k\_m A\_{at}^{5/6} \sqrt{\frac{T\_1}{ZT\_c}} \varrho^{-1/3}
$$

$$
k\_m \approx \quad 7.5 \text{ fs}, \quad T\_1 = 1 \text{ eV}. \tag{25}
$$

Here, *Aat* is the atomic number, and  is the mass density of matter in g/cm3. For example, for aluminum, *τ<sup>m</sup>* ≈ 11 fs with *Z* = 3, *-* = 2.71 g/cm3 and *Te* = 20 eV. We will not discuss the phenomenon of melting here, but we consider it as an example of a disordered configuration of ions.

(B) In the solid state, the interaction of the electrons with a perfect ionic lattice is already taken into account, introducing the band structure with the dispersion relation *Ep*,*ν*, Equation (12). Only deviations from the perfect lattice lead to scattering of the electron quasiparticles. We consider here the collective interactions of electrons with ion lattice vibrations (electron–phonon interaction). In addition, we have electron–electron collisions that will not contribute to the conductivity because of the conservation of the total momentum at the Coulomb interaction. However, Umklapp processes are possible, and they lead to a transfer of momentum from the electron system to the ion lattice.

#### *2.1. The Case of Individual Electron–Electron and Electron–Ion Interactions*

For this case, the Hamiltonian accounting for only the electronic degrees of freedom can be written as

$$H = \sum\_{p} E\_{p} \mathfrak{a}\_{p}^{\dagger} \mathfrak{a}\_{p} + \sum\_{pk} V\_{\text{el}}(k) \mathfrak{a}\_{p+k}^{\dagger} \mathfrak{a}\_{p}$$

$$\frac{1}{2} + \frac{1}{2} \sum\_{p\_1 p\_2 k} V\_{\text{ev}}(k) \mathfrak{a}\_{p\_1 + k}^{\dagger} \mathfrak{a}\_{p\_2 - k}^{\dagger} \mathfrak{a}\_{p\_2} \mathfrak{a}\_{p\_1 \text{\textquotedblleft}p\_1 \text{\textquotedblright}p} \tag{26}$$

with *Ep* = *h*¯ <sup>2</sup> *p*2/(2*m*). Only a single conduction band is considered in this subsection. The interactions between ions and electrons *V*ei(*k*) = *V*(*k*) are given by the Coulomb potential. *V*(*k*) = <sup>−</sup>*Zv*(*k*) and *<sup>V</sup>*ee(*k*) = *<sup>v</sup>*(*k*) = <sup>4</sup>*πe*2/*k*<sup>2</sup> is the potential of the *<sup>e</sup>* <sup>−</sup> *<sup>e</sup>* interaction. The ions can be treated in adiabatic approximation via the static ion structure factor *S*ii [23], which reflects the ion configuration. The ion component will be described in terms of an average charge number *Z* with the particle density *n*<sup>i</sup> = *n*/*Z* due to charge neutrality. The ion temperature is denoted as *T*i.

In a more general case, pseudo-potentials *V*ps ei (*k*) should be considered to take into account the structure of complex ions, the screening of the Coulomb potential, and the influence of bound and free electrons on their interaction with free electrons [12,24]. In particular, the expression for |*V*ei| 2, which appears in Born approximation can be rewritten as

$$\left|\left|V\_{\rm el}(\mathbf{k})\right|\right|^{2} = \left|V\_{\rm el}^{\rm p\*}(\mathbf{k})\right|^{2} S\_{\rm il}(\mathbf{k}).\tag{27}$$

One can use the following expression for the individual *<sup>e</sup>* <sup>−</sup> *<sup>i</sup>* interaction *<sup>V</sup>*ps ei (*k*) [12]:

$$V\_{\rm eli}^{\rm Fe}(k) = V(k)[1 + \varkappa\_{\rm ei}^2/k^2]^{-1}\cos(kr\_{\rm cut}),\tag{28}$$

where the prefactor *V*(*k*) is the electron–ion Coulomb potential. The second term on the r.h.s. is owing to the account of statical screening by the free electrons (the contribution of ions to screening is already taken into account by the ion structure factor *S*ii). It can be derived [12,25] in the low-frequency limit from the more sophisticated Lennard–Balescu approximation by summing up ring diagrams [8,22], leading to a characteristic inverse screening radius. It is κee = (4*πne*2/*Te*)1/2 in the classical limit. The third term is due to an empty-core model pseudo-potential [24,26], where *r*cut is a free parameter that can be fitted to match experimental data on transport and optical properties.

Similarly, in calculations within Born approximation, it is reasonable to replace the electron–electron Coulomb potential by a screened Coulomb potential:

$$V\_{\rm ce}(k) = v(k)[1 + \varkappa\_{\rm ce}^2/k^2]^{-1} \tag{29}$$

where the electron–electron Coulomb potential *v*(*k*) is statically screened with the screening parameter κee.

Substituting (26)– (29) into the correlation function [*H*<sup>ˆ</sup> , *<sup>n</sup>*<sup>ˆ</sup> *<sup>p</sup>*,*ν*]; [*H*<sup>ˆ</sup> , *<sup>n</sup>*<sup>ˆ</sup> *<sup>p</sup>*,*μ*]*ω*+*iη*, one gets from (24), the following expressions for the real and imaginary parts of the correlation functions C*eq* nm (where superscript "eq" denotes electron–electron ("q"="e") or electron–ion ("q"="i") contributions to the full correlation function C*nm* = Cee *nm* + Cei *nm*), see [12]:

$$\mathfrak{C}\_{nm}^{eq} = a\_{\neq} / (\mathfrak{Z}\pi w) \times \int\_0^{\infty} f^{eq}(y) dy \mathcal{R}\_{nm}^{eq} \left(\frac{w}{y}, y\right) \ln\left[\frac{1 + e^{c\_{\neq} - (w/y - y)^2}}{1 + e^{c\_{\neq} - (w/y + y)^2}}\right] \tag{30}$$

$$\mathfrak{C}\_{nm}^{\prime \prime \text{eq}} = \frac{\mathfrak{a}\_q}{3\pi^2 w} \int\_0^\infty f^{\text{eq}}(y) dy \left[ \sum\_{l=\pm 1} \mathcal{Z}\_{nm}^{\text{eq},l}(y) - 2\mathcal{Z}\_{nm}^{\text{eq},0}(y) \right] \tag{31}$$

where *α*<sup>i</sup> = *Z*, *α*<sup>e</sup> = 1/ <sup>√</sup>2, and <sup>C</sup>e*<sup>q</sup> nm* and <sup>C</sup>e*<sup>q</sup> nm* denote real and imaginary parts of correlation functions, respectively;

$$\mathcal{X}\_{nm}^{eq,l} = \int\_0^\infty \frac{d\xi}{\xi} \sum\_{\sigma=\pm 1} \sigma \mathcal{R}\_{nm}^{eq} \left(\xi + \sigma \frac{lw}{\mathcal{Y}}, \mathcal{Y}\right)$$

$$\times \ln\left[1 + e^{c\_\mu - \left[\frac{\sigma}{\xi} + \sigma(y + lw/y)\right]^2}\right],\tag{32}$$

with *<sup>l</sup>* <sup>=</sup> 0, <sup>±</sup>1. The factors *<sup>R</sup>*e*<sup>q</sup> nm* are the following polynomials for *n*, *m* = 1, 3, see [9]:

$$\begin{aligned} R\_{11}^{\text{ei}} &= 1\\ R\_{31}^{\text{ei}}(\mathbf{x}, \mathbf{y}) &= R\_{13}^{\text{ei}}(\mathbf{x}, \mathbf{y}) = 1 + \mathbf{y}^2 + 3\mathbf{x}^2\\ R\_{33}^{\text{ei}}(\mathbf{x}, \mathbf{y}) &= 2 + 2\mathbf{y}^2 + \mathbf{y}^4 + 2\mathbf{x}^2(5 + 3\mathbf{y}^2) + 9\mathbf{x}^4\\ R\_{11}^{\text{eve}} &= R\_{31}^{\text{eve}} = R\_{13}^{\text{eve}} = 0\\ R\_{33}^{\text{eve}}(\mathbf{x}, \mathbf{y}) &= 1 + 19\mathbf{x}^2/4. \end{aligned} \tag{33}$$

Similar expressions can be given for the higher order polynomials, see [7,23]. The following dimensionless units are used hereafter:

$$\begin{aligned} \check{k} &= k/k\_{\hat{\mathsf{R}}} \ y = \check{k}/(2\sqrt{2}), \ k\_{\hat{\mathsf{A}}}^{-1} = \Lambda = \hbar/\sqrt{mT\_{\varepsilon}};\\ \check{r} &= k\_{\hat{\mathsf{A}}}r; \ w = \hbar\omega/(4T\_{\varepsilon}) \end{aligned} \tag{34}$$

(*r* is any quantity having the dimension of a coordinate).

The functions *f* eq(*y*) occurring in Equation (30) are

$$f^{\rm el}(y) = f^{\rm el}\_{\rm scr}(y) \cos^2(2\sqrt{2}y\mathbb{F}\_{\rm cut})\mathcal{S}\_{\vec{n}}(y), \quad f^{\rm oc}(y) = f^{\rm oc}\_{\rm scr}(y) \tag{35}$$

where *f* ei scr(*y*) and *f* ee scr(*y*) are screening functions,

$$f\_{\rm scr}^{\rm el}(y) = y^3 / [y^2 + \kappa\_{\rm el}^2 / 8] \quad f\_{\rm scr}^{\rm oc}(y) = y^3 / [y^2 + \kappa\_{\rm oc}^2 / 4]. \tag{36}$$

Here,

$$\begin{aligned} \kappa\_{\rm el}^2 &= \min \left\{ \tilde{k}\_{\rm D}^2, \tilde{k}\_{\rm max}^2 \right\} \quad \kappa\_{\rm oc} = \tilde{k}\_{\rm D} \\ \tilde{k}\_{\rm max}^2 &= \mathcal{C}\_{\tilde{k}\_{\rm max}} 8\varepsilon\_{\rm F} / (18\pi Z)^{2/3}, \quad \tilde{k}\_{\rm D}^2 = \left[ \mathcal{R}\_{\rm D}^2 (1 + 2\varepsilon\_{\rm F}/3) \right]^{-1} \end{aligned} \tag{37}$$

where *<sup>R</sup>*˜ <sup>D</sup> <sup>=</sup> *<sup>R</sup>*D/*λ*¯ is the dimensionless Debye radius, *<sup>R</sup>*<sup>D</sup> <sup>=</sup> *<sup>V</sup>*th/*ω*pl, *<sup>V</sup>*th <sup>=</sup> <sup>√</sup>*Te*/*m*, and *<sup>C</sup>*˜ *<sup>k</sup>*max ≈ 1 is constant. Expressions for the ion–ion structure factor can be found in [12,27,28].

The value of ˜ *k*<sup>D</sup> leads to a proper interpolation between the degenerate and non-degenerate limits of screening and ensures a good agreement between Lennard–Balescu calculations and screened Born calculations everywhere, except for some frequency regions in the vicinity of the plasma frequency *ω*pl [12]. The parameter ˜ *k*max gives a respective restriction of screening for strongly coupled plasmas [12,29] at distances of the order of the interatomic distance

$$R\_0 = \left(4\pi n\_i/3\right)^{-1/3}.\tag{38}$$

In the vicinity *<sup>ω</sup>* <sup>∼</sup> *<sup>ω</sup>*pl, a more precise expression for the screening function *<sup>f</sup>* ei scr can be derived by comparing it with the Lennard–Balescu expression for the dynamical collision frequency *ν*1(*ω*), see [12]:

$$f\_{\mathsf{scr}}^{\mathsf{el}} = f\_{\mathsf{dyn}}(y, w) = \varepsilon\_{\mathsf{RPA}}^{\ast}(y, w) / \left[ y \varepsilon\_{\mathsf{RPA}}^{\prime}(y, 0) / \varepsilon\_{\mathsf{RPA}}(y, w) /^{2} \right] \tag{39}$$

where *ε*RPA is the RPA (random phase approximation) permittivity, *ε*∗ RPA is its complex conjugate, and *ε* RPA is its real part;

$$\varepsilon\_{\rm RPA}(y, w) = 1 + \frac{\sqrt{\omega\_{\rm uu}}}{8\sqrt{2}\pi} \frac{1}{y^3} \left[ -\sum\_{l=\pm 1} \mathcal{Z}\_{11}^l(y) \right.\\ \left. + i\pi \ln \left( \frac{1 + \exp[\varepsilon\_{\mu} - (w/y - y)^2]}{1 + \exp[\varepsilon\_{\mu} - (w/y + y)^2]} \right) \right] \tag{40}$$

$$\mathcal{X}\_{11}^{\mathbb{I}} = \int\_0^\infty \frac{d\xi^\mathbb{v}}{\xi^\mathbb{v}} \sum\_{\sigma=\pm 1} \sigma \ln\left[1 + e^{c\_\mu - \left[\frac{\pi}{2} + \sigma(y + lw/y)\right]^2}\right].\tag{41}$$

*l* = 0, ±1.

Taking in mind that in the considered case of dynamical screening the screening function *f*scr(*y*, *w*) (39) is a complex function (note that in the case of statical screening *f* ei scr is dependent only on *y*, rather than on *y* and *w*), one can rewrite the expressions for real and imaginary parts of the correlation functions stipulated by electron–ion interactions as

$$\begin{aligned} \mathfrak{C}\_{nm}^{\prime \mathrm{el}} &= \mathfrak{C}\_{nm}^{\prime \mathrm{el}}(f\_{\mathrm{dyn}}^{\prime}) - \mathfrak{C}\_{nm}^{\prime \mathrm{el}}(f\_{\mathrm{dyn}}^{\prime \prime}) \\ \mathfrak{C}\_{nm}^{\prime \prime \mathrm{el}} &= \mathfrak{C}\_{nm}^{\prime \prime \mathrm{el}}(f\_{\mathrm{dyn}}^{\prime}) + \mathfrak{C}\_{nm}^{\prime \mathrm{el}}(f\_{\mathrm{dyn}}^{\prime \prime}) \end{aligned} \tag{42}$$

where designations Cei *nm*(*f* dyn) and <sup>C</sup>ei *nm*(*f* dyn) (where superscripts "ei" designate *e* − *i* interactions) mean that in the respective Expressions (30) and (31) for the real and imaginary part of the correlation function, respectively, the real part of the screening function (39) is substituted by *f* ei scr in (36), and similarly designations Cei *nm*(*f* dyn) and <sup>C</sup>ei *nm*(*f* dyn) indicate the substitution of the imaginary part of *<sup>f</sup>*dyn (39) by *<sup>f</sup>* ei scr in (36).

With the same arguments as for static screening (36), one can suppose that *f* dyn should be replaced by

$$f\_{\text{scr,min}} = y^3 / \left[ y^2 + \vec{k}\_{\text{max}}^2 / 8 \right] \tag{43}$$

where ˜ *k*max is given by (37) if *f* dyn <sup>&</sup>lt; *<sup>f</sup>*scr,min.

It should be noted that the restriction of screening at distances <sup>≈</sup> *<sup>R</sup>*<sup>0</sup> occur at <sup>Γ</sup>ei,g <sup>&</sup>gt; 1/3, where

$$\Gamma\_{\rm el,g} = \frac{\Gamma\_{\rm el}}{\epsilon\_{\rm F} \left[ 4(1 + e^{-\epsilon\_{\mu}(\epsilon\_{\rm F})}) / (3\sqrt{\pi t}) \right]^{2/3}} \tag{44}$$

is the generalized electron–ion coupling parameter for plasma at arbitrary degeneracy [12], and Γei = *Ze*2/(*R*0*Te*) is the classical plasma coupling parameter. In the non-degenerate case (<sup>F</sup> 1), the relation (44) gives the conventional expression <sup>Γ</sup>ei,g <sup>≈</sup> <sup>Γ</sup>ei <sup>=</sup> *Ze*2/(*R*0*Te*), while in the case of strongly degenerate (<sup>F</sup> 1) plasmas, the generalized coupling parameter depends on the Fermi energy, <sup>Γ</sup>ei,g <sup>≈</sup> (9*π*/16)1/3*Ze*2/(*R*0*EF*) <sup>≈</sup> 1.21*Ze*2/(*R*0*EF*).

#### *2.2. Collective Electron–Phonon Interactions and Electron–Electron Interactions via Umklapp Processes*

For ion temperatures *<sup>T</sup>*<sup>i</sup> <sup>&</sup>lt; *<sup>T</sup>*<sup>m</sup> or when a lattice is heated to *<sup>T</sup>*<sup>i</sup> <sup>&</sup>gt; *<sup>T</sup>*<sup>m</sup> during times *<sup>t</sup>* <sup>&</sup>lt; *<sup>τ</sup>m*, where *<sup>τ</sup><sup>m</sup>* is given by (25), one can assume that the electrons interact with lattice vibrations (phonons). In addition, electron–electron interactions can contribute to the change in the energy of the electron gas through electron momentum transfer to the lattice (Umklapp processes). This situation can be described by the Hamiltonian

$$\begin{split} \hat{H} &= \sum\_{k,i,\sigma} E\_{k,i} \hat{a}\_{k,i,\sigma}^{+} \hat{a}\_{k,i,\sigma} + \sum\_{q,\lambda} \hbar \omega\_{q,\lambda} \hat{b}\_{q,\lambda}^{+} \hat{b}\_{q,\lambda} \\ &+ \sum\_{k,q,i,i',\lambda,\sigma} g\_{k}(q,i,i',\lambda) \hat{a}\_{k+q,i,\sigma}^{+} \hat{a}\_{k,i',\sigma} (\hat{b}\_{q,\lambda}^{+} + \hat{b}\_{-q,\lambda}) \\ &+ \frac{\mathcal{U}}{2n} \sum\_{k,k',q,j,i,\sigma} \hat{a}\_{k+q-g,i,\sigma}^{+} \hat{a}\_{k'-q,i,-\sigma}^{+} \hat{a}\_{k',i,-\sigma} \hat{a}\_{k,i,\sigma} \end{split} \tag{45}$$

where the first two terms on the r.h.s. of (45) represent electron and phonon kinetic energies, respectively. The third term represents the electron–phonon interaction in the Fröhlich form [30]. The fourth term represents the electron–electron interaction accounting for Umklapp processes in the Hubbard [31] form, where *g* is the wave vector of the inverse lattice. *i*, *i* are electron band numbers, *λ* is the phonon mode number, *<sup>σ</sup>* is the spin quantum number, "*a*<sup>+</sup> *k*,*i* , "*ak*,*i*, "*b*<sup>+</sup> *<sup>q</sup>*,*λ*,"*b*−*q*,*<sup>λ</sup>* are the creation and annihilation operators of electrons and phonons, respectively, *Ek*,*<sup>i</sup>* and *ωq*,*<sup>λ</sup>* are the energy of electrons in the *i*-th band, given by (12) and the frequency of phonons in *λ*-th mode, respectively, and *gk*(*q*, *i*, *i* , *λ*) is the coefficient of the electron–phonon interaction. *U* is the single-site approximation to the Coulomb interaction of electrons with opposite spin orientations.

This effective Hamiltonian (45) can be derived from a more fundamental Hamiltonian describing the Coulomb system consisting of electrons and ions. The phonon excitations are obtained from the dynamical structure factor, and the Hubbard term is obtained from the *e* − *e* interaction. For our exploratory calculations, we consider the case of electrons interacting with a single phonon mode of longitudinal optical (LO) phonons with a frequency independent of the wave vector [30]. LO phonons have been considered for simplicity. One can expect that a consideration of acoustical phonons (see, e.g., [32]) gives similar results for ion temperatures greater than the Debye one.

We also disregard in this subsection *interband* transitions between different bands *n*, *i* and consider only the contribution due to free–free transitions of electrons with effective mass *mi* = *m*<sup>∗</sup> in the conduction band (the consideration of a two-level system for the case of electron–phonon interactions is given in the subsequent subsection). Besides that, we disregard cross terms from contributions of different types of interactions (normal and Umklapp processes) when calculating the correlation functions of the first

moment of the density operator. Calculating commutators in (24) and making respective transformations, one obtains in the Born approximation

$$\mathfrak{C}\_{11} = \mathfrak{C}\_{\mathfrak{e}-ph} + \mathfrak{C}\_{\text{U}} \tag{46}$$

where

$$\begin{aligned} \mathfrak{C}\_{\text{e-ph}} &= \ \frac{i\epsilon\_{\text{F}}^{-3/2}}{2\pi^{5/2}} \frac{w\_{\text{lo}}}{\widehat{w}\_{\text{a}\text{in}}^{1/2}} m\_{\text{e}}^{2} \mathbb{C}\_{\text{eph}} \int ^{\infty} y dy \ \int\_{-\infty}^{\infty} dx \\ & \times \left[ \frac{1}{w - x + w\_{\text{lo}} + i\eta} + \frac{1}{w + x - w\_{\text{lo}} + i\eta} \right] \\ & \times \frac{(e^{4x} - 1)^{-1} - (e^{4ww\_{\text{IO}}} - 1)^{-1}}{x - aw\_{\text{lo}}} \\ & \times \ln\left[ \frac{1 + \exp[\epsilon\_{\mu} - (y - x/y)^{2}]}{1 + \exp[\epsilon\_{\mu} - (y + x/y)^{2}]} \right] \end{aligned} \tag{47}$$

is the contribution due to electron–phonon interactions, derived earlier in [13], where *C*eph is a constant of the order of unity, *<sup>α</sup>* <sup>=</sup> *Te*/*Ti*; *<sup>w</sup>* <sup>=</sup> *<sup>h</sup>*¯ *<sup>ω</sup>*/(4*Te*), *<sup>ω</sup>*au <sup>=</sup> *<sup>h</sup>*¯ *<sup>ω</sup>*au/*Te*, *<sup>w</sup>*LO <sup>=</sup> *<sup>h</sup>*¯ *<sup>ω</sup>*LO/(4*Te*) with the frequency of the longitudinal optical phonons *ω*LO, and *η* is an infinitesimal small value.

The contribution of Umklapp processes to the correlation function is given by the second term in Equation (46),

$$\mathfrak{C}\_{\rm U} = \frac{9im\_{\ast}}{4} \frac{\mathcal{U}^2 T\_{\varepsilon}^2}{E\_F^3 E\_H} \sum\_{\mathbf{g}} \mathcal{g}^2 I\_{\Omega}(\mathbf{g}) I\_{\mathbb{E}}(\mathcal{W}, \mathfrak{e}\_{B\prime} \mathfrak{e}\_{\mu}) \tag{48}$$

$$\mathcal{J}\_{\Omega}(\mathbf{g}) = \frac{1}{(4\pi)^4} \iiint d\Omega \, d\Omega' d\Omega\_1 d\Omega'\_1 \times \delta\left(\mathbf{k}\_1 + \mathbf{k}'\_1 - \mathbf{k} - \mathbf{k}' - \mathbf{g}\right) \tag{49}$$

$$\mathrm{J}\_{\mathrm{E}}(\mathcal{W}, \boldsymbol{\varepsilon}\_{\mathrm{B}}, \boldsymbol{\varepsilon}\_{\mathrm{\mu}}) = \iiint\limits\_{-\boldsymbol{\varepsilon}\_{\mathrm{\mu}}} \frac{n\_{\mathrm{x}\_{1}} n\_{\mathrm{x}\_{1}'} n\_{\mathrm{x}} n\_{\mathrm{x}'} \mathrm{d} \mathbf{x}\_{1} d \mathbf{x}\_{1}' d \mathbf{x} \, d \mathbf{x}'}{\mathrm{x}\_{1} + \mathrm{x}\_{1}' - \mathrm{x} - \mathrm{x}' + \mathcal{W} + i\eta} \times \frac{\mathbf{e}^{\mathbf{x}\_{1} + \mathbf{x}\_{1}'} - \mathbf{e}^{\mathbf{x} + \mathbf{x}'}}{\mathrm{x}\_{1} + \mathrm{x}\_{1}' - \mathrm{x} - \mathrm{x}'}. \tag{50}$$

In (49), the vectors **g**, **k**1, **k** <sup>1</sup>, **k**, and **k** are taken dimensionless, as the ratio to the absolute value of the Fermi wave vector *kF* = *kF* = (3*π*2*n*)1/3. The integration in (49) is performed on solid angles for each of the respective wavevectors **k**1, **k** <sup>1</sup>, **k**, **k** . Therefore, the integral (49) is dependent only on **g**. Furthermore, *W* = *h*¯ *ω*/*Te* = 4*w*;  *<sup>B</sup>* = *EB*/*Te*, and *EB* is the energy of electrons on the surface of the 1st Brillouin zone boundary. In (50), *nx* = 1/[1 + *ex*] is the Bose distribution function.

The integral over *dx* in Equations (47) and (48) can be performed applying the Sokhotski–Plemelj formula

$$\int\_{-\infty}^{\infty} \frac{f(\mathbf{x})}{\mathbf{x} + i\eta} d\mathbf{x} = -i\pi f(0) + \mathcal{P} \int\_{-\infty}^{\infty} \frac{f(\mathbf{x})}{\mathbf{x}} d\mathbf{x} \Big|$$

where P denotes the principal value of the integral. In this case, real and imaginary parts of the correlation function can easily be derived. Below we consider the real parts of the correlation functions (46):

$$\begin{split} \mathbf{c}\_{\mathbf{e}-ph}^{\prime} &= \frac{\mathbf{c}\_{\mathbf{F}}^{-3/2}}{2\pi^{3/2}} \frac{m\_{\ast}^{2}\mathbf{C}\_{\mathbf{e}\text{ph},\mathbf{w}\_{\text{LO}}}}{\hat{w}^{1/2}} \\ &\times \sum\_{\boldsymbol{\sigma}=\pm 1} \frac{(\epsilon^{4}[\boldsymbol{w}\_{\text{LO}} + \boldsymbol{\sigma}\boldsymbol{w}] - 1)^{-1} - (\epsilon^{4}\boldsymbol{w}\_{\text{LO}}\boldsymbol{a} - 1)^{-1}}{w\_{\text{LO}}(\boldsymbol{a} - 1) - \boldsymbol{\sigma}\boldsymbol{w}} \\ &\times \int \mathbf{y} dy \ln \left[ \frac{1 + \exp\left\{\epsilon\_{\mu} - \left[y - \frac{\boldsymbol{w}\_{\text{LO}} + \boldsymbol{\sigma}\boldsymbol{w}}{y}\right]^{2}\right\}}{1 + \exp\left\{\epsilon\_{\mu} - \left[y + \frac{\boldsymbol{w}\_{\text{LO}} + \boldsymbol{\sigma}\boldsymbol{w}}{y}\right]^{2}\right\}} \right] \end{split} \tag{51}$$

and

$$\begin{array}{rcl} \mathfrak{E}'\_{\mathbb{U}} &=& \frac{9}{4} \mathbb{C}\_{\mathbb{U}} N\_{\mathbb{X}} \bar{g}^2 m\_\* \frac{\underline{\mathbb{U}}^2 T\_{\varepsilon}^2}{E\_F^3 E\_H} l'\_{\mathbb{E}} \\\ j'\_{\mathbb{E}} &=& \frac{4}{W} \int\_{-2\varepsilon\_{\mu}}^{2\varepsilon\_{\Lambda} - 2\varepsilon\_{\mu}} dt \left[ \frac{1}{\varepsilon^{t-W} - 1} - \frac{1}{\varepsilon^t - 1} \right] \\\ & \times \ln \left[ \frac{\varepsilon^{t/2} + \varepsilon^{-B/2}}{\varepsilon^{t/2 - B/2} + 1} \right] \ln \left[ \frac{\varepsilon^{t/2 - W/2} + \varepsilon^{-B/2}}{\varepsilon^{t/2 - B/2 - W/2} + 1} \right] \end{array} \tag{52}$$

where <sup>Δ</sup> = Δ*E*/*Te*, where Δ*<sup>E</sup>* is the energetic distance between the Fermi surface and the Brillouin zone boundary, *B* = *μ* + <sup>Δ</sup>.

While deriving Equation (52) from Equation (48), the substitution *x* = (*t* + *r*)/2, *x* = (*t* − *r*)/2, *x*<sup>1</sup> = (*t*<sup>1</sup> + *r*1)/2, *x* <sup>1</sup> = (*t*<sup>1</sup> − *r*1)/2 was done, the integrals over *r* and *r*<sup>1</sup> were calculated explicitly, and the *δ*-function was accounted for while integrating over *t*1. To derive Expression (49), the approximation

$$\sum\_{\mathfrak{X}} \mathfrak{g}^2 J\_{\Omega}(\mathfrak{g}) \approx \mathcal{C}\_{\mathbf{U}\mathfrak{J}} \mathfrak{Z} N\_{\mathfrak{X}}.$$

is made. *C*<sup>U</sup> is a constant of the order of 1, which can be found, e.g., from optical measurements, as in [10,33,34]. *Ng* is the number of different wavevectors of the inverse lattice in the first Brillouin zone, which coincides with the number of nearest neighbors in the inverse lattice for the point g = 0. ¯ *g*<sup>2</sup> is the average value of *g*<sup>2</sup> in the first Brillouin zone.

#### *2.3. Interband Transitions in a Two-Level System*

In this subsection, we consider expressions for the force–force correlation functions of first order moments <sup>C</sup>{1,*ν*}{1,*μ*} = <sup>C</sup>*νμ* for different energy bands *<sup>ν</sup>*, *<sup>μ</sup>* = 1, 2. The expressions will be derived for the case of collective electron–phonon interactions, when the electronic system Hamiltonian is described by (45) with the account of only the first three terms, i.e., without Umklapp processes. With this Hamiltonian and with an electron–phonon coupling function for longitudinal optical phonons, one can obtain from Equation (24) the following expression for correlation functions C*νμ* similarly, as was done above:

(i) for the correlation functions related to intra-band electron transitions:

$$\mathfrak{C}\_{\mathsf{VV}} = -im\_{\bullet}^{2}\mathfrak{x}\int\_{0}^{\infty}y\,dy\int\_{-\infty}^{\infty}\mathrm{d}x\mathrm{X}\_{\mathsf{VV}}(\mathbf{x})\,\Delta\mathbf{F}\_{\mathsf{VV}}(\mathbf{x},\mathsf{y}) - \sum\_{\mu\neq\nu}\mathfrak{C}\_{\mu\nu}^{0}\tag{53}$$

where *<sup>κ</sup>* <sup>=</sup> <sup>−</sup>3/2 F 2*π*5/2 *w*LO *<sup>ω</sup>*1/2 *a*.*u*. *C*eph; C<sup>0</sup> *μν* stands for the indirect influence of interband transitions:

$$\begin{cases} \mathfrak{C}^{0}\_{\mu\nu} = \inf\_{\omega \succsim \mathsf{X}} \lambda / 4 \bigcap\_{0 \precsim \infty}^{\infty} \mathsf{Ax} \Big\{ \mathsf{X}\_{\mu\nu}(\mathbf{x}) \\ \times \left[ \left( \frac{\mathbf{x}^{2}}{y^{3}} + \frac{2\mathbf{x}}{y} + y \right) \Delta F\_{\mu\nu}(\mathbf{x}, y) + \frac{m\_{\star}^{-1}}{y} \Delta \widetilde{F}\_{\mu\nu}(\mathbf{x}, y) \right] \\ + \mathsf{X}\_{\nu\mu}(\mathbf{x}) \\ \times \left[ \left( \frac{\mathbf{x}^{2}}{y^{3}} - \frac{2\mathbf{x}}{y} + y \right) \Delta F\_{\nu\mu}(\mathbf{x}, y) + \frac{m\_{\star}^{-1}}{y} \Delta \widetilde{F}\_{\nu\mu}(\mathbf{x}, y) \right] \end{cases} \tag{54}$$

where

$$X\_{\mu\nu}(\mathbf{x}) = \left[ \frac{1/q\_{\nu\mu}(\mathbf{x},a)}{w + i\eta + q\_{\nu\mu}(\mathbf{x},1)} + \frac{1/q\_{\nu\mu}(\mathbf{x},a)}{w + i\eta - q\_{\nu\mu}(\mathbf{x},1)} \right] \times \left[ \frac{1}{e^{4x + 4w\_{\mu\nu}} - 1} - \frac{1}{e^{4w\_{\text{LO}}} - 1} \right] \tag{55}$$

where *ϕνμ*(*x*, *<sup>t</sup>*) = *<sup>x</sup>* <sup>+</sup> *<sup>w</sup>νμ* <sup>−</sup> *tw*LO, *<sup>w</sup>νμ* <sup>=</sup> <sup>1</sup> <sup>4</sup> (*E<sup>ν</sup>* − *Eμ*)/*Te*;

$$\begin{aligned} \Delta F\_{\mu\nu}(\mathbf{x}, \boldsymbol{y}) &= F\_1(A\_{-}^{\mu}(\mathbf{x}, \boldsymbol{y})) - F\_1(A\_{+}^{\nu}(\mathbf{x}, \boldsymbol{y})), \\ \Delta \tilde{F}\_{\mu\nu}(\mathbf{x}, \boldsymbol{y}) &= F\_2(A\_{-}^{\mu}(\mathbf{x}, \boldsymbol{y})) - F\_2(A\_{+}^{\nu}(\mathbf{x}, \boldsymbol{y})) \\ A\_{\pm}^{\mu}(\mathbf{x}, \boldsymbol{y}) &= \exp\left[\left(\mathbf{x}/\boldsymbol{y} \pm \boldsymbol{y}\right)^2 + \left(E\_{\mu} - \mu\right)/T\_{\varepsilon}\right] \\ F\_1(t) &= \ln\left(1 + \frac{1}{T}\right), \quad F\_2(t) = \ln^2(t) + 2\operatorname{Li}\_2(-t). \end{aligned} \tag{56}$$

(Note that in the above expression, in the term *E<sup>μ</sup>* − *μ*, the second value *μ* denotes the chemical potential, and the index *μ* in the first term *Eμ* denotes the energy band).

(ii) for the correlation functions related to interband electrons transitions:

$$\mathcal{C}\_{\mu\nu}^{\mu\rho\tau\nu} = \mathrm{im}\_{\ast}^{2}\mathrm{x} / 4 \int\_{0}^{\infty} dy \int\_{-\infty}^{y} \mathrm{d}x \left\{ X\_{\mu\nu}(\mathbf{x}) \times \left[ \left( \frac{\mathbf{x}^{2}}{y^{3}} - y \right) \Delta F\_{\mu\nu}(\mathbf{x}, y) + \frac{m\_{\ast}^{-1}}{y} \Delta F\_{\mu\nu}(\mathbf{x}, y) \right] + \text{the same with } \mu \leftrightarrow \nu \right\}. \tag{57}$$

#### **3. Calculations and Discussion**

We perform exploratory calculations of the absorption coefficient of aluminum plasmas with the solid density *-* = 2.71 g/cm<sup>3</sup> and the average ion charge *Z* = 3 (step-like density profile), irradiated by lasers with different wavelengths shown in Figures 1 and 2. The dielectric function (16) was determined for two cases: an ordered lattice with an account of collective electron–phonon interactions and Umklapp processes by means of Expressions (46), (51), and (52), assuming *Ti* <sup>&</sup>lt; *<sup>T</sup>*melt or *<sup>t</sup>* <sup>&</sup>lt; *<sup>τ</sup><sup>m</sup>* (25), and a disordered lattice, with an account of individual electron–ion and electron–electron interactions by means of Expressions (23), (30), and (31), with respect to a Percus–Yevick-like model [16] for *Sii*, with the restrictions of screening (37) (see also [12]).

The effective mass *m*<sup>∗</sup> was calculated according to the Huttner model [35], see also [11]. The value of *C*eph ≈ 5.73 in (47) was chosen to reproduce the low-frequency cold reflectivity of aluminum [11,33] for the laser wavelength *λ* = 0.4 μm. The value of *ω*LO was determined by the position of the maximum of the phonon spectrum for aluminum [36] as *h*¯ *ω*LO ≈ 0.006 eV. Furthermore, the following parameters have been used: *<sup>C</sup>*<sup>U</sup> = 1.5, *<sup>U</sup>* = 2 eV, <sup>Δ</sup>*<sup>E</sup>* = 7.3 eV, and ¯ *g*<sup>2</sup> = 2. The value *Ng* = 8 is chosen, assuming the fcc lattice structure and the bcc inverse lattice structure for aluminum.

**Figure 1.** Absorption coefficient (first row) as well as real (second row) and imaginary (third row) parts of the complex collision frequency of electrons (46), as a function of the electron temperature *Te*, for a solid-density aluminum plasma (*-* = 2.71 g/cm3) and laser radiation of different wavelengths *λ* (left column: 80 nm; center: 0.4 μm; right: 10 μm). Curves are shown for different ion temperatures (see legend), calculated for the case of an ordered lattice by Models (46), (51), and (52) (with an account of electron–phonon interactions and Umklapp processes, labeled by "e-ph"), as well as for the case of a disordered lattice (electron–ion and electron–electron interactions, labeled by "e-i"). Curves with the markers "\*" and "" are according to the models cited in the text (see legend).

Figure 1 shows the dependence of the absorption coefficient *A*, Equation (2), the real and imaginary parts of the effective collision frequency *ν*(*ω*) (17) on the electron temperature. Ordered (electron–phonon) and disordered ion lattices (*e* − *i*) are considered at different wavelengths of laser radiation (*λ* = 0.08, 0.4, 10 μm). The results for *λ* = 0.4 μm are compared with a semi-empirical model by Povarnitsyn et al. [33,37] and with calculations of the absorption coefficient taken from the work of Cauble et al. [38].

The most essential feature of the dependence of the effective collision frequency on electron temperature, for the case of ordered ion lattice with electron–phonon interactions and Umklapp processes, is that for *Te* <sup>&</sup>gt; *EF* (*EF* <sup>≈</sup> 11.6 eV for solid-density aluminum) the real part *<sup>ν</sup>* is decreasing as

$$\vee \sim \tau\_e^{-4}$$

according to the asymptotic behavior of Expressions (52), see [14]. This is much faster than the scaling

$$\nu' \sim T\_e^{-3/2}$$

for the case of electron–ion interaction in a high-temperature plasma [12]. This is shown by the solid and marked curves in Figure 1.

It should be noted that, in the case of an ordered ion lattice, unlike the case of a disordered ion lattice, the temperature dependence of the real part of the effective collision frequency as well as the absorption coefficient shows a clear peak in the vicinity of the Fermi temperature for all wavelengths. The maximum values of *A*, calculated by Models (46), (51), and (52) for an ordered ion lattice, are relatively close to that of the semi-empirical model [33] and to the values of *A*, calculated by Models (23), (30), and (31) for a disordered lattice.

It should also be noted that under the conditions of Figure 1 the contribution of Umklapp processes (52) exceeds the contribution of the electron–phonon contribution (51). The influence of different ion temperatures is mainly owing to the dependence of the effective electron mass *m*<sup>∗</sup> on the ion temperature *Ti*, as calculated according to the Huttner model [35] (compare the thin and thick solid curves in Figure 1).

**Figure 2.** The absorption coefficient and the real part of the complex collision frequency of electrons (46) in a solid-density aluminum plasma are shown as a function of the laser frequency *h*¯ *ω* for ion temperatures *Ti* = 0.04 eV and different electron temperatures *Te* (see legend). Thick curves (left column) refer to the case of an ordered lattice with electron–phonon interactions and Umklapp processes, and thin curves (right column) refer to a disordered lattice with electron–ion and electron–electron interactions.

The imaginary part of *ν*(*ω*) calculated by a consequent quantum statistical model (for semi-empirical models, it can be calculated by equating the respective expressions for permittivity to the Drude-like expression (16), see [12]) is small for laser frequencies *h*¯ *ω EF* or for *h*¯ *ω EF*, while for *h*¯ *ω* ∼ *EF* its value can be considerable, see the subfigure at the left corner of Figure 1.

For *<sup>h</sup>*¯ *<sup>ω</sup> EF* and *<sup>h</sup>*¯ *<sup>ω</sup>* <sup>&</sup>lt; *Te*, Re *<sup>ν</sup>*(*ω*) (and hence the value of the density of absorbed power of laser radiation, which is proportional to Im *ε* ∼ Re *ν*(*ω*)) is weakly dependent on *ω*, see Figure 2. For *h*¯ *ω Te*

and *h*¯ *ω* - *EF*, the value of Re *ν*(*ω*) is increasing with *ω* (the rate of the increase is proportional to *ω*2) for Umklapp processes in an ordered ion lattice, see [13]). For *h*¯ *ω Te* and *h*¯ *ω EF*, the value of Re *ν*(*ω*) is decreasing with *ω*, see Figure 2. The rate of such a decrease is proportional to *ω*−3/2 for the case of a disordered ion lattice [12] and to *ω*−<sup>1</sup> for the case of an ordered ion lattice [13]. One should note here that the contribution of interband transitions to the dielectric function should be taken into account for quantum energies exceeding the gap between the conductivity band and inner electron energy levels [17]. As described above, the respective theoretical approach can be elaborated on the basis of the nonequilibrium statistical operator method. Note also that in the case of a disordered ion lattice (individual electron–ion and electron–electron interactions), Expressions (30) and (31) for the correlation functions give rise to the Ziman–Evans formula [39,40] for the electric conductivity in the limit *ω* → 0 [12].

#### **4. Conclusions**

The NSO method allows one to describe consistently the dielectric function of WDM in a wide range of frequencies and plasma parameters. Different kinds of electron–ion interactions, in particular in systems with a disordered distribution of ions (individual electron–ion and electron–electron interactions) and systems with an ordered ion configuration, described as an ion lattice (with collective electron–phonon interactions and Umklapp processes), are successfully treated using this method.

A main peculiarity of the absorption of irradiated laser energy in the case of an ordered ion lattice, owing to collective electron–phonon interactions and Umklapp processes, is a much stronger (<sup>∼</sup> *<sup>T</sup>*−<sup>4</sup> *<sup>e</sup>* ) decrease of the real part of the effective collision frequency at electron temperatures exceeding the Fermi temperature, if compared to the case of individual electron–ion and electron–electron interactions (where Re *<sup>ν</sup>* <sup>∼</sup> *<sup>T</sup>*−3/2 *<sup>e</sup>* for *Te* <sup>&</sup>gt; *EF*). An essential feature of the electron–phonon interaction and the Umklapp process is that they show, as a function of the electron temperature, a clear peak structure of the absorption coefficient and of the real part of the effective collision frequency at *Te* ∼ *EF*. In addition, in both cases of ordered and disordered ion lattices, the real part of the effective collision frequency as a function of the photon energy shows a peak structure: Re *<sup>ν</sup>*(*ω*) is nearly independent of *<sup>ω</sup>* for *<sup>h</sup>*¯ *<sup>ω</sup> EF* and *<sup>h</sup>*¯ *<sup>ω</sup>* <sup>&</sup>lt; *Te*, but Re *<sup>ν</sup>*(*ω*) is rising with *<sup>ω</sup>* (<sup>∼</sup> *<sup>ω</sup>*<sup>2</sup> for the case of ordered ion lattice) for *<sup>h</sup>*¯ *<sup>ω</sup> Te*, *<sup>h</sup>*¯ *<sup>ω</sup>* - *EF*, and Re *ν*(*ω*) is decreasing with *<sup>ω</sup>* (<sup>∼</sup> *<sup>ω</sup>*−<sup>1</sup> for the case of ordered ion lattice and <sup>∼</sup> *<sup>ω</sup>*−3/2 for the case of disordered ion lattice) for ¯*hω Te* and ¯*hω EF*.

The NSO method can be also applied for the description of interband contributions to the dielectric function, which are essential for photon energies exceeding the energy gaps between the conduction band and electron bands corresponding to excited electron energy levels. The relation to DFT calculation [17], which provides us with ab initio calculations of optimal single-electron states and replaces the approximation for the electron–ion potential used in our calculations, is of high interest. In addition, our approach allows one to take into account the contribution of electron–electron collisions to the dielectric function, which is not possible in mean-field approaches such as DFT. More detailed investigations of this case, as well as the calculation of the imaginary part of the effective collision frequency for the case of an ordered ion lattice with electron–phonon interactions and Umklapp processes, will be the subject of our future work.

**Author Contributions:** H.R. and G.R. have been working on the general concepts of the approach. M.V. has conducted the particular derivations and calculations.

**Funding:** The work of M.V. was partly supported by the Presidium of the RAS program "Thermo-physics of High Energy Density", No I.13P. H.R. and G.R. acknowledge support by the DFG (German research foundation) via RE 975/4-1 and RO 905/32-1, respectively. G.R. acknowledges support from the MEPhI Academic Excellence Program under contract number 02.a03.21.0005.

**Acknowledgments:** We are greatfull to N.E. Andreev for valuable comments

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**



c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
