Quantum Dynamics of Charged Fermions in the Wigner Formulation of Quantum Mechanics

: To study the kinetic properties of dense quantum plasma, a new quantum dynamics method in the Wigner representation of quantum mechanics has been developed for extreme conditions, when analytical approximations based on different kinds of perturbation theories cannot be applied. This method combines the Feynman and Wigner formulation of quantum mechanics and uses for calculation the path integral Monte-Carlo (WPIMC) in phase space and quantum generalization of the classical molecular dynamics methods (WMD) allowing to solve the quantum Wigner–Liouville-like equation. The Fermi–Dirac statistical effects are accounted for by the effective pair pseudopotential depending on coordinates and momenta and allowing to avoid the famous “sign problem” due to realization of the Pauli blocking of fermions. Signiﬁcant inﬂuence of the interparticle interaction on the high energy asymptotics of the momentum distribution functions have been observed. According to the quantum Kubo formula, we also study the electron conductivity of dense plasma media.


Introduction
The main advantage of the Wigner formulation of the quantum mechanics in the phase space is that it allows obtaining more information on the quantum system than can be done by any other quantum description. This is important in many fields such as the quantum information processing, quantum electronics, quantum chemistry, signal processing, dynamical properties of many body quantum systems and many others. Moreover, this approach allows treating the entanglement that can occur in quantum systems (for review, see [1]).
In this paper, we apply the Wigner formulation of the quantum simulation of strongly coupled low-temperature two-component plasma (about 10,000-100,000 K). This allows us to avoid the problem that the thermal radiation [2] causes substantial energy loss in weakly coupled high-temperature plasmas [3][4][5] such as thermal nuclear fusion plasmas. For instance, the well-known Lawson condition is derived from the condition that the fusion energy compensates the thermal bremsstrahlung emission loss. During the last decades, many papers have been devoted to studying the thermodynamic and kinetic properties of strongly coupled plasmas. A two component plasma is quantum mechanical even at high temperature and low density since the Heisenberg uncertainty principle is necessary to keep the electrons from collapsing into ions. Some of the most powerful numerical methods for simulation of quantum systems are the Monte-Carlo methods (PIMC), based on the path integral formulation of quantum mechanics [6]. The main difficulty for the PIMC studies of Fermi systems results from the requirement of anti-symmetrization of the density matrix [6]. As a result, all thermodynamic quantities are presented as the sum of terms with alternating sign related to even and odd permutations and are equal to the small difference of two large numbers, which are the sums of positive and negative terms. The numerical calculation in this case is severely hampered. This difficulty is known in the literature as the "sign problem". To overcome this issue, many approaches have been developed (for example, [7,8]).
The second main disadvantage of the PIMC method is that it cannot cope with the problem of calculation of average values of arbitrary quantum operators in phase space and momentum distribution functions, while this problem may be central in treatment of thermodynamic and kinetic properties of matter. The main results of this work can be formulated as follows: to overcome disadvantages of the standard PIMC methods in configuration space, we use a new numerical approach (WPIMC) based on combining the path integral and the Wigner formulation of quantum mechanics [9,10]; to account for the Fermi-Dirac statistical effects in explicit formm the effective pair pseudopotential depending on coordinates and momenta has been used, thus avoiding the famous "sign problem" due to realization of the Pauli blocking of fermions in phase space at the finite temperature; the high energy asymptotics of the momentum distribution functions has been studied under strong interparticle Coulomb interaction; and, according to the quantum Kubo formula, the electron conductivity of strongly-coupled hydrogen plasma has been calculated using the Wigner generalization of the standard classical molecular dynamics method.

Wigner Quantum Dynamics
As example of Coulomb system of particles, we consider a 3D two-component mass asymmetric electron-hole plasma consisting of N e electrons and N h heavier holes in equilibrium (N e = N h = N) [11]. The Hamiltonian of the systemĤ =K +Û c contains the kinetic energyK and the Coulomb interaction energyÛ c =Û c hh +Û c ee +Û c eh contributions. Our starting point is the canonical ensemble-averaged time correlation function [12] whereF andÂ are operators of arbitrary observables, t c = t − ihβ/2 is the complex time, β = 1/k B T and Z = Tr e −βĤ is the partition function. The Wigner representation of Equation (1) in a υ-dimensional space is where ν = 12N, A(pq) and F( pq) denote the Weyl's symbols of the operators and W(pq; pq; t; β) is the spectral density (the Wigner-Liouville-like function) expressed as As has been proved in [13][14][15], W obeys the following integral equation: W (pq; pq; t; β) = dp 0 q 0 d p 0 q 0 G (pq, pq, t; p 0 q 0 , p 0 q 0 , 0) W(p 0 q 0 ; p 0 q 0 ; t = 0, β) ds dp q d p q G pq, pq, t; p q , p q , t with the Green function G pq, pq, t; p q , p q , t = δ p − p(t; p q , t ) δ q − q(t; p q , t ) describing propagation of the spectral density along classical trajectories in positive time direction dp(t; p q , t ) dt and in the reverse time direction where ( q t ) = [ q(t; p q , t )] and, similarly for bared quantities, while vectors F = −∇U c and v symbolize the forces and particle velocities, respectively, ω(s, q) = These equations of motion are supplemented by initial conditions at time t = 0 and by initial conditions at time t = t In fact, Equation (4) are Hamiltonian equations of motion but written for half-time (t/2). Similarly, Equation (5) are half-time Hamiltonian equations of motion reversed in time. The time correlation is taken between instants in the past and the future with the initial conditions fixed in between these instants, i.e., at t = 0 the spectral density is W(pq 0 ; pq 0 ; t = 0, β) = W 0 (pq 0 ; pq 0 ; β). The right-hand sides of Equations (4) and (5) include interparticle interaction, which can be arbitrary strong.
The solution of the integral in Equation (3) can be represented by an iteration series where W t is the initial quantum spectral densities evolving classically during time intervals [0, t], are operators that govern the propagation from time τ i to τ i+1 such as the first term in Equation (3) and momentum jumps due to the convolution structure of Equation (3) (see, e.g., [13][14][15]). Thus, the time correlation function becomes where φ(pq; pq) ≡ F(pq)A( pq) and the parentheses (. . . | . . . ) denote integration over the phase spaces {p 0 q 0 ; p 0 q 0 }, {dp q d p q } and so on. To compute the electron-electron conductivity, we calculate the electron time correlation function C vv (t) and then apply the Kubo formula which contains the Fourier transform of C vv (t) at ω = 0.
The iteration series for C FA (t) can be efficiently computed using MC methods. We have developed a MC scheme which provides domain sampling of the terms giving the main contribution to the iteration series (cf. [13][14][15]). For simplicity, in this work, we take into account only the first term of iteration series, which is related to the propagation of the initial quantum distribution W 0 according to forward and backward in time the Hamiltonian equations of motion. This term, however, does not describe pure classical dynamics but accounts for quantum effects [13] and, in fact, contains arbitrarily high powers of the Planck's constant:

The Initial Wigner Function for Canonical Ensemble
Of course, the exact matrix elements of density matrix of interacting quantum systems is not known (particularly for low temperatures and high densities), but they can be constructed using a path integral approach based on the operator identity The anti-symmetrized Wigner function can be written in the form: Here h and l, t = 1, . . . , N a ) are spin degrees of freedom and the spatial coordinates of electrons and positive particles, respectively. Index m = 0, . . . , M − 1 labels the off-diagonal high-temperature density matrices; each high temperature factor can be presented in the form ) with the error of order 1/M 2 arising from neglecting commutator 2 [K, U c ] /2. In the limit M → ∞, the error of the whole product of high temperature factors is equal to zero (∝ 1/M) and we have the exact path integral representation of the Wigner function.
The spin gives rise to the spin part of the density matrix (S) with exchange effects accounted for by the permutation operatorsP e andP h acting on the electron and hole coordinates q (M) and the spin projections σ . The sum is over all permutations with parity κ P e and κ P h .
Accounting for only identical and pair permutations, the final expression for the Wigner function can be obtained in the form of the path integral over all closed trajectories and can be presented in the form [2]: where v a lt ≈ −kT ln 1 − δ σ l,a σ t,a exp −2π|q l,a − q t,a | 2 exp − |(p l,a −p t,a )| 2 (2π) 2 α 2 p t,a = p t,a + 2 while α is small parameter. The constant C(M) cancels out when we calculate the average values of operators. We imply that momenta and coordinates are dimensionless variables such as p l,aλa /h and q l,a /λ a , whereλ a = 2πhβ m a M . The Pauli blocking of fermions is accounted for in Equation (11) by the exchange pseudopotential in phase space v a lt and provides agreement of the momentum distribution for ideal fermions with analytical Fermi-Dirac distribution in the wide ranges of degeneracy and momenta, where decay of the distribution functions covers at least five orders of magnitude [2].
Here, each particle is represented by a discrete trajectory consisting of a set of M points ("beads") and the whole configuration of the particles is represented by a 3(N e + N h )M-dimensional vector q ≡ {q N h ,h }. When → 0, this multiple integral tends to the exact representation of the Wigner function W(pq; β) in the form of path integral with continuous dimensionless "imaginary time" τ [6], which corresponds to m/M in discrete case. In addition, the set of independent variables q (m) tends into closed trajectory q(τ). This trajectory starts and ends in 0 when τ = 0 and 1. Let us note that integration here is related to the integration over the Wiener measure of all closed trajectories q(τ) [16]. In fact, a particle is presented by the trajectory with characteristic size of order λ a = 2πhβ m a in coordinate space. This is a manifestation of the uncertainty principle.
In Equation (11) , U = U hh + U ee + U eh denotes the sum of all interaction energies, each consisting of the corresponding sum of pair interactions given by the well known Kelbg potential Φ ab defined by the following expression [17,18]: where e a and e b are charges of particles, q m lt = |q m l,a − q m t,b |/λ ab , λ ab = h 2 /(2m ab ), m ab = m a m b /(m a + m b ) is the reduced mass of the (ab)-pair of particles and the error function is defined by erf(x) = 2 √ π x 0 dte −t 2 . Here, coordinates are written in the natural units.

Results of Numerical Calculations
The plasma density is characterized by the Brueckner parameter r s defined as the ratio of the mean interparticle distance of electrons d = ( 3 4πn e ) 1/3 to the Bohr radius a B (n e is the electron densities). To estimate the importance of interaction and degeneracy effects, we use the degeneracy χ = n e λ 3 e and the coupling Γ = e 2 /(r s a B k B T) parameters, where λ e = 2πh 2 /(m e k B T) and m e is the electron mass.
The momentum distribution functions: We define the momentum distribution function for holes (a = h) and electrons (a = e) by the following expressions: where δ is delta function. The momentum distribution function w a (|p|) gives a probability density for particle of type a to have momentum p. For classical systems of particles, due to the commutativity of kinetic and potential energy operators, we have the Maxwell distribution function (MD) proportional to exp(−(pλ a ) 2 /4πh 2 ), even under the strong coupling. Quantum ideal systems of particles, due to the quantum statistics, have the Fermi-Dirac (FD) or Bose-Einstein momentum distribution functions. Interaction of a quantum particle with its surroundings restricts the volume of configuration space, which can affect the shape of momentum distribution function due to the uncertainty relation. Through these effects in frameworks of some models and perturbation theories, the momentum distribution functions have the power-law "tails" (const/p 8 ) even under conditions of thermodynamic equilibrium [19][20][21][22][23][24][25]. In particular, it is shown that the momentum distribution in the asymptotic region of large momenta p may be described by the sum of the Maxwell distribution (MD) and quantum correction proportional to const/p 8 (we use short notation as P8) or sum of the MD and product of const/p 8 and the Maxwell distributions with effective temperature that exceeds the temperature of medium (short notation -P8MDEF) [25]. Reliable analytical calculations of numerical parameters in these asymptotic is problematic.
The WPIMC calculations for the electron and hole momentum distribution functions in electronhole plasma are presented in Figure 1 (left panel) by Lines 1 and 4 for electrons and holes, respectively. Here, the parameter of electronic degeneracy is equal to four (T/E F = 0.261), while the hole mass is two times larger than the electron mass. Figure 1 also presents several related analytical dependencies: FD (Lines 2 and 5), MD (Lines 3 and 6), P8 (Lines 7 and 8), and P8MDEF (Lines 9 and 10).
Here, the constant in P8 and the effective temperature in P8MDEF have been used as the adjustable parameters to fit the WPIMC momentum distribution functions in the large momentum asymptotic regions.
Transport coefficients: A natural way to obtain transport coefficients is making use of the quantum Green-Kubo relations [12]. These relations give the transport coefficients in terms of integrals of equilibrium time-dependent correlation functions. According to Equation (8), the electron conductivity σ is the integral of the velocity-velocity autocorrelation function where the scalar product of 3D-velocities is and the trajectories in positive (bared) and inverse (tilded) time directions are defined by Equations (4) and (5), respectively. Calculations of autocorrelation functions are performed in canonical ensemble and include combination of the Monte-Carlo sampling of initial conditions pq 0 and pq 0 for trajectories and solving the system of dynamic Hamiltonian Equations (4) and (5). The initial conditions pq 0 and pq 0 for the trajectories are sampled by Monte-Carlo method according to the modulus of probability W 0 (pq 0 ; pq 0 ; β), while sign of the W 0 (pq 0 ; pq 0 ; β) is accounted for as a wieght function at calculations average values [2].

(central panel) shows examples of the velocity-velocity autocorrelation functions and its antiderivatives. The initial
Wigner distribution W 0 , due to the path integral representation, accounts for momentum-coordinate principle uncertainty. Consequently, the initial momenta are to a great extent independent form each other and the VVACF approach to zero. Some time later, correlation of the VVACF starts to grow due to strong interaction with surrounding particles. This happens as Γ e ≤ 1, |q 0 − q 0 | ∼ λ e and λ e ≤ βe 2 ≤ r s a B and the virtual trajectories pq and pq are approximately continuations of each other independently from initial conditions. Subsequent decay of the VVACF results from interaction with far particles at distances larger than r s . According to Figure 1, the damping time of the VVACF turns out to be strongly affected by the variations of density and temperature. Let us discuss now the conductivity of a strongly coupled Coulomb system. Figure 1 (right) presents comparison of our conductivity isochores with those obtained from the interpolation formula for conductivity of fully ionized hydrogen plasma derived in [26,27]. This figure demonstrates agreement of both data at low densities and high temperatures. Let us stress that Line 6 restricts approximately from above the region of conductivity sharp drop due to arising bound states of many particle clusters. show results of this work for the fixed r s = 6, 4, 3, 2, 1, respectively, while related lines present results of interpolation formula for conductivity of fully ionized hydrogen plasma [26,27]. Line 6 restricts approximately from above the region of conductivity sharp drop due to arising bound states of many particle clusters.

Discussion
In this work, we use the generalization of the classical molecular dynamics methods allowing to take into account quantum effects. We use the Feynman and Wigner formulations of quantum mechanics combining with Monte Carlo methods for numerical treatment of the kinetic properties of dense quantum plasma media. The Wigner-Liouville-like equation is solved by a combination of Wigner molecular dynamics (WMD) and by Monte-Carlo method in phase space (WPIMC) methods. The initial Wigner-like function has been presented in the form of path integrals. Fermi statistical effects are accounted for by effective pair pseudopotential depending on coordinates, momenta and degeneracy parameter of particles and taking into account Pauli blocking of fermions. WPIMC allows calculating thermodynamic quantities and pair distribution functions in a wide range of density and temperature.
Comparison of the classical Maxwell-Boltzmann and quantum Fermi-Dirac distribution shows the significant influence of the interparticle interaction on the high energy asymptotics of the momentum distribution functions resulting in appearance of the power-law quantum "tails". To study the influence of the interparticle interaction on the dynamic properties of dense plasmas, we also compute the temporal correlation functions of quantum operators. According to the quantum Kubo formula, we determined the electron conductivity and compared obtained results with available theories. Our results show a strong dependence on the plasma parameters and for fully ionized plasma are in a agreement with available theories, simulations, experimental data and interpolation formula obtained by Esser, Redmer and Röpke. Appearance of bound states and many particle clusters in plasma results in sharp drop of electron conductivity.

Abbreviations
The following abbreviations are used in this manuscript:

PIMC
Monte Carlo methods based on the path integral formulation of quantum mechanics WPIMC Path integral Monte Carlo method in phase space WMD Quantum generalization of the classical molecular dynamics methods MD Maxwell distribution function FD Fermi-Dirac momentum distribution function P8 Quantum correction to the momentum distribution function proportional to const/p 8 P8MDEF Quantum correction multiplied on Maxwell distributions with effective temperature VVACF Velocity-velocity auto correlation function