Open access peer-reviewed chapter

# Nonequilibrium Statistical Operator

Written By

Gerd Röpke

Submitted: 27 May 2018 Reviewed: 24 January 2019 Published: 10 May 2019

DOI: 10.5772/intechopen.84707

From the Edited Volume

## Non-Equilibrium Particle Dynamics

Edited by Albert S. Kim

Chapter metrics overview

View Full Metrics

## Abstract

Nonequilibrium statistical physics is concerned with a fundamental problem in physics, the phenomenon of irreversibility, which is not rigorously solved yet. Different approaches to the statistical mechanics of nonequilibrium processes are based on empirical assumptions, but a rigorous, first principle theory is missing. An important contribution to describe irreversible behavior starting from reversible Hamiltonian dynamics was given by Zubarev, who invented the method of the nonequilibrium statistical operator (NSO). We discuss, in particular, the extended von Neumann equation and the entropy concept in this approach. The method of NSO proved to be a general and universal approach to different nonequilibrium phenomena. Typical applications are the quantum master equation, kinetic theory, and linear response theory which are outlined and illustrated solving standard examples for reaction and transport processes. Some open questions are emphasized.

### Keywords

• Zubarev method
• NSO
• quantum master equation
• kinetic equations
• linear response theory

## 1. Introduction: irreversible processes

### 1.1 Irreversibility and arrow of time

Irreversibility belongs to the unsolved fundamental problems in recent physics. Nonequilibrium processes are omnipresent in our daily experience. However, a fundamental, microscopic description of such processes is missing yet.

Our microscopic description of physical phenomena is expressed by equations of motion, well known in mechanics, electrodynamics, quantum mechanics, and field theory. We model a physical system, we determine the degrees of freedom and the forces, and we introduce a Lagrangian. The equations of motion are differential equations. If we know the initial state, the future of the system can be predicted solving the equations of motion. Anything is determined. The equations of motion are invariant with respect to time reversion. The time evolution is reversible. No arrow of time is selected out, nothing happens what is not prescribed by the initial state.

This picture was created by celestial dynamics. It is very successful, very presumptuous, and many processes are described with high precision. However, it is in contradiction to daily experience. We know birth and death, decay, destruction, and many other phenomena that are irreversible, selecting out the arrow of time.

A qualitative new discipline in physics is thermodynamics. It considers not a model but any real system. The laws of thermodynamics define new quantities, the state variables. The second law determines the entropy S as state variable (and the temperature T ) via

dS = 1 T δ Q reversible E1

where δQ is the heat imposed to the system within a reversible process, together with the third law which fixes the value S T = 0 = 0 independent on other state variables. For irreversible processes holds

dS dt > δQ T . E2

In particular, for isolated system, δQ = 0 , irreversible processes are possible so that dS / dt > 0 . Typical examples are friction that transforms mechanical energy to thermal energy, temperature equilibration without production of mechanical work, diffusion processes to balance concentration gradients. An arrow of time is selected out, time reversion describes a phenomenon which is not possible. How can irreversible evolution in time be obtained from the fundamental microscopic equations of motion which are reversible in time?

For equilibrium thermodynamics, a microscopic approach is given by statistical physics. Additional concepts are introduced such as probability and distribution function, ensembles in thermodynamic equilibrium, and information theory. New thermodynamic quantities are introduced, basically the entropy, which have no direct relation to mechanical quantities describing the equation of motion. However, nonequilibrium processes are described in a phenomenological way, and no fundamental solution of the problem of irreversibility is found until now. A substantial step to develop a theory of irreversible evolution is the Zubarev method of the nonequilibrium statistical operator (NSO) [1, 2, 3, 4, 5, 6] to be described in the following section. It is a consistent theory to describe different nonequilibrium processes what is indispensable for a basic approach.

### 1.2 Langevin equation

To give an example for a microscopic approach to a nonequilibrium process, let us consider the Brownian motion. A particle suspended in a liquid, moving with velocity v medium , experiences a friction force F fric t ,

d dt v t = 1 m F fric t = γ v t v medium , E3

with the coefficient of friction γ . The solution

v t = v t 0 e γ t t 0 + v medium 1 e γ t t 0 E4

describes the relaxation from the initial state v t 0 at t 0 to the final state v medium for t t 0 . Independent of the initial state, the particle rests in equilibrium with the medium. In the general case not considered here, an external force can be added.

As it is well known, this simple relaxation behavior cannot be correct because it does not describe the Brownian motion, showing the existence of fluctuations also in thermal equilibrium. This problem was solved with the Langevin equation: instead of the trajectory v t as solution of a differential equation, the stochastic process V t is considered. It obeys the stochastic differential equation

d dt V t = γ V t v rel t + R t . E5

The random acceleration R t (or the stochastic force m R t ) is a stochastic process, which is characterized by special properties. For instance, white noise is a Gaussian process that is characterized by the mean value R t = 0 and the auto-correlation function

R i t 1 R j t 2 = φ ij t 2 t 1 = 2 D δ ij δ t 2 t 1 . E6

D is the diffusion coefficient. An interesting result is the Einstein relation (fluctuation-dissipation theorem, FDT)

D γ = k B T m E7

which relates the friction coefficient γ (dissipation) to the fluctuations φ in the system (stochastic forces), characterized by the parameter D ; see [5] for more details.

### 1.3 Von Neumann equation

Within statistical mechanics, the thermodynamic state of an ensemble of many-particle systems at time t is described by the statistical operator ρ t . We assume that the time evolution of the quantum state of the system is given by the Hamiltonian H t which may contain time-dependent external fields. The von Neumann equation follows as equation of motion for the statistical operator,

t ρ t + i H t ρ t = 0 . E8

The von Neumann equation describes reversible dynamics. The equation of motion is based on the Schrödinger equation. Time inversion and conjugate complex means that the first term on the left-hand side as well as the second one change the sign, since i i and both the Hamiltonian and the statistical operator are Hermitean. However, the von Neumann equation is not sufficient to determine ρ t because it is a first-order differential equation, and an initial value ρ t 0 at time t 0 is necessary to specify a solution. This problem emerges clearly in equilibrium.

### 1.4 Thermodynamic equilibrium and entropy

By definition, in thermodynamic equilibrium, the thermodynamic state of the system is not changing with time. Both, H t and ρ t , are not depending on t so that

t ρ eq t = 0 . E9

The solution of the von Neumann equation in thermodynamic equilibrium becomes trivial, i H ρ eq = 0 . The time-independent statistical operator ρ eq commutes with the Hamiltonian. We conclude that ρ eq depends only on constants of motion C n that commute with H . But, the von Neumann equation is not sufficient to determine how ρ eq depends on constants of motion C n . We need a new additional principle, not included in Hamiltonian dynamics.

Equilibrium statistical mechanics is based on the following principle to determine the statistical operator ρ eq : consider the functional (information entropy)

S inf ρ = Tr ρ ln ρ E10

for arbitrary ρ that is consistent with the given conditions T r ρ = 1 (normalization) and

Tr ρ C n = C n E11

(self-consistency conditions). Respecting these conditions, we vary ρ and determine the maximum of the information entropy for the optimal distribution ρ eq so that δ S inf [   ρ eq ] = 0 . As it is well-known, the method of Lagrange multipliers can be used to account for the self-consistency conditions (11). The corresponding maximum value for S inf ρ

S eq ρ eq = k B Tr ρ eq ln ρ eq E12

is the equilibrium entropy of the system at given constraints C n and k B is the Boltzmann constant. The solution of this variational principle leads to the Gibbs ensembles for thermodynamic equilibrium, see also Section 4.

As an example, we consider an open system which is in thermal contact and particle exchange with reservoirs. The corresponding equilibrium statistical operator has to obey the given constraints: normalization Tr ρ = 1 , thermal contact with the bath so that Tr ρ H = U (internal energy), particle exchange with a reservoir so that for the particle number operator N c of species c , the average is given by Tr ρ N c = n c Ω , where Ω denotes the volume of the system (we do not use V to avoid confusion with the potential), and n c is the particle density of species c . Looking for the maximum of the information entropy functional with these constraints, one obtains the grand canonical distribution

ρ eq = e β H c μ c N c Tr e β H c μ c N c . E13

The normalization is explicitly accounted for by the denominator (partition function). The second condition means that the energy of a system in heat contact with a thermostat fluctuates around an averaged value H = U = u Ω with the given density of internal energy u . This condition is taken into account by the Lagrange multiplier β that must be related to the temperature, a more detailed discussion leads to β = 1 / k B T . Similarly, the contact with the particle reservoir fixes the particle density n c , introduced by the Lagrange multiplier μ c , which has the meaning of the chemical potential of species c .

Within the variational approach, the Lagrange parameters β , μ c have to be eliminated. This leads to the equations of state ( eq = Tr { ρ eq } ) which relate, e.g., the chemical potentials μ c to the particle densities n c ,

H eq = U Ω β μ c , N c eq = Ω n c T μ c . E14

The entropy S eq Ω β μ follows from Eq. (12). The dependence of extensive quantities on the volume Ω is trivial for homogeneous systems. After a thermodynamic potential is calculated, all thermodynamic variables are derived in a consistent manner. The method to construct statistical ensembles from the maximum of entropy at given conditions, which take into account the different contacts with the surrounding bath, is well accepted in equilibrium statistical mechanics and is applied successfully to different phenomena, including phase transitions.

Can we extend the definition of equilibrium entropy (12) also for ρ t that describes the evolution in nonequilibrium? Time evolution is given by an unitary transformation that leaves the trace invariant. Thus, the expression Tr ρ t ln ρ t is constant for a solution ρ t of the von Neumann equation

d d t Tr ρ t ln ρ t = 0 . E15

The entropy for a system in nonequilibrium, however, may increase with time, according to the second law of thermodynamics. The equations of motion, including the Schrödinger equation and the Liouville-von Neumann equation, describe reversible motion and are not appropriate to describe irreversible processes. Therefore, the entropy concept (12) elaborated in equilibrium statistical physics together with the Liouville-von Neumann equation cannot be used as fundamental approach to nonequilibrium statistical physics.

## 2. The method of nonequilibrium statistical operator (NSO)

After the laws of thermodynamics have been formulated in the nineteenth century, in particular, the definition of entropy for systems in thermodynamic equilibrium and the increase of intrinsic entropy in nonequilibrium processes, a microscopic approach to nonequilibrium evolution was first given by Ludwig Boltzmann who formulated the kinetic theory of gases [7] using the famous Stoßzahlansatz. The question how irreversible evolution in time can be obtained from reversible microscopic equations has been arisen immediately and was discussed controversially.

The rigorous derivation of the kinetic equations from a microscopic description of a system was given only a long time afterward by Bogoliubov [8] introducing a new additional theorem, the principle of weakening of initial correlation.

### 2.1 Construction of the Zubarev NSO

A generalization of this principle has been given by Zubarev [1, 2], who invented the method of the nonequilibrium statistical operator (NSO). This approach has been applied to various problems in nonequilibrium statistical physics, see [3, 4] and may be considered as a unified, fundamental approach to nonequilibrium systems which includes different theories such as quantum master equations, kinetic theory, and linear response theory to be outlined below. An exhaustive review of the Zubarev NSO method and its manifold applications cannot be given here, see [1, 2, 3, 4, 5].

In the first step, we interrogate the concept of thermodynamic equilibrium. This is an idealization, because slow processes are always possible. As example, we may take the nuclear decay of long-living isotopes, hindered chemical reactions, or the long-time relaxation of glasses. Concepts introduced for equilibrium have to be generalized to nonequilibrium. An example is the thermodynamics of irreversible processes.

#### 2.1.1 The relevant statistical operator

A solution of the problem to combine equilibrium thermodynamics and nonequilibrium processes was proposed by Zubarev [1, 2]. To characterize the nonequilibrium state of a system, we introduce the set of relevant observables B n extending the set of conserved quantities C n . At time t , the observed values B n t have to be reproduced by the statistical operator ρ t , i.e.,

Tr ρ t B n = B n t . E16

However, these conditions are not sufficient to fix ρ t , and we need an additional principle to find the correct one in between many possible distributions which all fulfill the conditions (16). We ask for the most probable distribution at time t , where the information entropy has a maximum value (see Section 4)

δ Tr ρ rel t ln ρ rel t = 0 E17

with the self-consistency conditions

Tr ρ rel t B n = B n t E18

and Tr ρ rel t = 1 . Once more, we use Lagrange multipliers λ n t to account for the self-consistency conditions (18). Since the averages are, in general, time dependent, the corresponding Lagrange multipliers are now time-dependent functions as well. We find the generalized Gibbs distribution

ρ rel t = e Φ t n λ n t B n , Φ t = ln Tr e n λ n t B n , E19

where the Lagrange multipliers λ n t (thermodynamic parameters) are determined by the self-consistency conditions (18). Φ t is the Massieux-Planck function, needed for normalization purposes and playing the role of a thermodynamic potential. Generalizing the equilibrium case, Eq. (12), we can consider the relevant entropy in nonequilibrium

S rel t = k B Tr ρ rel t ln ρ rel t . E20

Relations similar to the relations known from equilibrium thermodynamics can be derived. In particular, the production of entropy results as

S rel t t = n λ n t B ̇ n t E21

as known from the thermodynamics of irreversible processes. In contrast to Eq. (15), this expression can have a positive value so that S rel t can increase with time.

The relevant statistical operator ρ rel t is not the wanted nonequilibrium statistical operator ρ t because it does not obey the Liouville-von Neumann equation. Also, S rel t is not the thermodynamic entropy because it is based on the arbitrary choice of the set B n of relevant observables, and not all possible variables are correctly reproduced. As example, we consider below the famous Boltzmann entropy that is based on the single-particle distribution function, but does not take higher order correlation functions into account.

### 2.2 The Zubarev solution of the initial value problem

The solution of the problem how to find the missing signatures of ρ t not already described by ρ rel t was found by Zubarev [1, 2] generalizing the Bogoliubov principle of weakening of initial correlations [8]. He proposed to use the relevant statistical operator ρ rel t 0 at some initial time t 0 as initial condition to construct ρ t ,

ρ t 0 t = U t t 0 ρ rel t 0 U t t 0 . E22

The unitary time evolution operator U t t 0 is the solution of the differential equation

iℏ t U t t 0 = H t U t t 0 , E23

with the initial condition U t 0 t 0 = 1 . This unitary operator is known from the solution of the Schrödinger equation. If the Hamiltonian is not time dependent, we have

U t t 0 = e i H t t 0 . E24

If the Hamiltonian is time dependent, the solution is given by a time-ordered exponent.

Now, it is easily shown that ρ t 0 t is a solution of the von Neumann equation. All missing correlations not contained in ρ rel t 0 are formed dynamically during the time evolution of the system. However, incorrect initial correlations contained in ρ rel t 0 may survive for a finite time interval t t 0 , and the self-consistency conditions (18) valid at t 0 are not automatically valid also at t .

To get rid of these incorrect initial correlations, according to the Bogoliubov principle of weakening of initial correlations, one can consider the limit t 0 . According to Zubarev, it is more efficient to average over the initial time so that no special time instant t 0 is singled out. This is of importance, for instance, if there are long-living oscillations determined by the initial state. According to Abel’s theorem, see [1, 2, 3, 4], the limit t 0 can be replaced by the limit ϵ + 0 in the expression

ρ ϵ t = ϵ t e ϵ t 1 t U t t 1 ρ rel t 1 U t t 1 d t 1 . E25

This averaging over different initial time instants means a mixing of phases so that long-living oscillations are damped out. Finally, we obtain the nonequilibrium statistical operator as

ρ NSO t = lim ϵ 0 ρ ϵ t . E26

This way, ρ rel t 1 for all times < t 1 < t serves as initial condition to solve the Liouville-von Neumann equation, according to the Bogoliubov principle of weakening of initial correlations. The missing correlations are formed dynamically during the time evolution of the system. The more information about the nonequilibrium state are used to construct the relevant statistical operator, the less dynamical formation of the correct correlations in ρ t is needed. The limit t 0 is less active to produce the remaining missing correlating. The past that is of relevance, given by the relaxation time τ , becomes shorter, if the relevant (long-living) correlations are already correctly implemented. The limit ε + 0 has to be performed after the thermodynamic limit, see below.

### 2.3 Discussion of the Zubarev NSO approach

#### 2.3.1 The extended Liouville-von Neumann equation

The nonequilibrium statistical operator ρ ϵ t , Eq. (25), obeys the extended von Neumann equation

ρ ϵ t t + i H t ρ ϵ t = ϵ ρ ϵ t ρ rel t . E27

as can be seen after simple derivation with respect to time. In contrast to the von Neumann equation (8), a source term arises on the right-hand side that becomes infinitesimal small in the limit ϵ + 0 . This source term breaks the time inversion symmetry, so that for any finite value of ϵ , the solution ρ ϵ t describes in general an irreversible evolution with time.

The source term can be interpreted in the following way:

1. The source term implements the “initial condition” in the equation of motion as expressed by ρ rel t . Formally, the source term looks like a relaxation process. In addition to the internal dynamics, the system evolves toward the relevant distribution.

2. The construction of the source term is such that the time evolution of the relevant variables is not affected by the source term (we use the invariance of the trace with respect to cyclic permutations),

t B n t = Tr ρ ϵ t t B n = Tr i [ H t ρ ϵ t ] B n = i H t B n t = B ̇ n t . E28

The source term cancels because of the self-consistency conditions (18). Thus, the time evolution of the relevant observables satisfies the dynamical equations of motion according to the Hamiltonian H t .

1. The value of ϵ has to be small enough, ϵ 1 / τ , so that all relaxation processes to establish the correct correlations, i.e., the correct distribution of the irrelevant observables, can be performed. However, ℏϵ has to be large compared to the energy difference of neighbored energy eigenstates of the system so that mixing is possible. For a system of many particles, the density of energy eigenvalues is high so that we can assume a quasi-continuum. This is necessary to allow for dissipation. The van Hove limit means that the limit ϵ + 0 has to be performed after the thermodynamic limit.

2. Differential equations can have degenerated solutions. For instance, we know the retarded and advanced solution of the wave equation that describes the emission of electromagnetic radiation. An infinitesimal small perturbation can destroy this degeneracy and select out a special solution, here the retarded one. Similar problems are known for systems (magnetism) where the ground state has a lower symmetry than the Hamiltonian.

3. Any real system is in contact with the surroundings. The intrinsic dynamics described by the Hamiltonian H t is modified due to the coupling of the open system to the bath. Within the quantum master equation approach, we can approximate the influence term describing the coupling to the bath by a relaxation term such as the source term. At present, we consider the source term as a purely mathematical tool to select the retarded solution of the von Neumann equation, and physical results are obtained only after performing the limit ϵ 0 .

#### 2.3.2 Selection of the set of relevant observables

The Zubarev method to solve the initial value problem for the Liouville-von Neumann equation is based on the selection of the set B n of relevant observables which characterize the nonequilibrium state. The corresponding relevant statistical operator ρ rel t is some approximation to ρ t . According to the Bogoliubov principle of weakening of initial correlations, the missing correlations to get ρ t are produced dynamically. This process, the dynamical formation of the missing correlations, needs some relaxation time τ . If we would take instead of ρ rel t the exact (but unknown) solution ρ t , the relaxation time τ is zero. The Liouville-von Neumann equation, which is a first-order differential equation with respect to time, describes a Markov process.

There is no rigorous prescription how to select the set of relevant observables B n . The more relevant observables are selected so that their averages with ρ rel t reproduce already the correctly known averages B n t , see Eq. (18), the less the effort to produce the missing correlations dynamically, and the less relaxation time τ is needed. Taking into account that usually perturbation theory is used to treat the dynamical time evolution (23), a lower order of perturbation theory is then sufficient. We discuss this issue in Section 3.

In conclusion, the selection of the set of relevant observables is arbitrary, as a minimum the constants of motion C n have to be included because their relaxation time is infinite, their averages cannot be produced dynamically. The resulting ρ NSO t (26) should not depend on the (arbitrary) choice of relevant observables B n if the limit ε 0 is correctly performed. However, usually perturbation theory is applied, so that the result will depend on the selection of the set of relevant observables. The inclusion of long-living correlations into B n allows to use lower order perturbation expansions to obtain acceptable results.

#### 2.3.3 Entropy of the nonequilibrium state

An intricate problem is the definition of entropy for the nonequilibrium state. In nonequilibrium, entropy is produced, as investigated in the phenomenological approach to the thermodynamics of irreversible processes, considering currents induced by the generalized forces.

Such a behavior occurs for the relevant entropy defined by the relevant distribution (20),

S rel t = k B Tr ρ rel t ln ρ rel t . E29

A famous example that shows the increase of the relevant entropy with time is the Boltzmann H theorem where the relevant observables to define the nonequilibrium state are the occupation numbers of the single-particle states, i.e., the distribution function, see Section 3.2 for discussion.

Note that the increase of entropy cannot be solved this way. It is related to so-called coarse graining. The information about the state is reduced because the degrees of freedom to describe the system are reduced. This may be an averaging in phase space over small cells. The loss of information then gives the increase of entropy. This procedure is artificial, anthropomorphic, depending on our way to describe the details of a process.

The method of nonequilibrium statistical operator ρ NSO t allows to extend the set of relevant observables arbitrarily so that the choice of the set of relevant observables seems to be irrelevant. All missing correlations are produced dynamically. We can start with any set of relevant operators, but have to wait for a sufficient long time to get the correct statistical operator, or to go to very small ϵ . A possible definition of the entropy would be

S NSO t = k B Tr ρ NSO t ln ρ NSO t . E30

The destruction of the reversibility of the von Neumann equation (27) is connected with the source term on the right-hand side that produces the mixing by averaging over the past in Eq. (25). This source term is responsible for the entropy production. At present, there is no proof that the entropy S NSO t will increase also in the limit ϵ + 0 .

## 3. Applications

The NSO method is a fundamental step in deriving equations of evolution to describe nonequilibrium phenomena. It can be shown that any currently used description can be deduced from this approach. We give three typical examples, the quantum master equations, see [9, 10], kinetic theory, see [11], and linear response theory, see [12]. In all of these applications, we have to define the set of relevant observables, and to eliminate the Lagrange parameters determined by the self-consistency conditions. We shortly outline these applications, for a more exhaustive presentation see [3, 4, 5].

### 3.1 Quantum master equation

#### 3.1.1 Open systems

The main issue is that any physical system cannot be completely separated from the surroundings, so that the isolated system is only a limiting case of the open system which is in contact with a bath. More general, we subdivide the degrees of freedom of the total system into the relevant degrees of freedom which describe the system S under consideration, and the irrelevant part describing the bath B. Examples are a harmonic oscillator coupled to a bath consisting of harmonic oscillators, such as an oscillating molecule interacting with phonons or photons, or radiation from a single atom embedded in the bath consisting of photons, see below.

The Hamiltonian H of the open system can be decomposed

H = H S + H B + H int . E31

The system Hamiltonian acts only in the Hibert space of the system states leaving the bath states unchanged. It is expressed in terms of the system observables A ν . The bath Hamiltonian acts only in the Hilbert space of the bath states leaving the system states unchanged. It is expressed in terms of the bath observables B μ . Both sets of operators are assumed to be hermitean and independent so that A ν B μ = 0 .

We project out the relevant part of the nonequilibrium statistical operator ρ t

ρ s t = Tr B ρ t E32

where the trace over the bath can be performed after the eigenstates of the bath are introduced. The operator Tr B means the trace over the quantum states of the heat bath. If we have no further information, we construct the relevant statistical operator taking the equilibrium distribution ρ B = ρ eq (13) for the irrelevant degrees of freedom,

ρ rel t = ρ s t ρ B . E33

#### 3.1.2 Born-Markov approximation

Starting with the extended Liouville-von Neumann equation (27), we perform the trace T r B over the variables of the bath (see Eq. (32)),

t ρ s t 1 i H s ρ s t = 1 i Tr B H int ρ t E34

since the remaining terms disappear and 1 i Tr B H B ρ t ρ t H B = 0 because of cyclic invariance of the trace Tr B . To obtain a closed equation for ρ s t , the full nonequilibrium statistical operator ρ t occurring on the right-hand side has to be eliminated.

For this, we calculate the time evolution of the irrelevant part of the statistical operator Δ ρ t = ρ t ρ r el t ,

t Δ ρ t = t ρ t t ρ s t ρ B E35

inserting the time evolution for ρ t (8) and ρ s t (34) given above:

t + ε Δ ρ t = 1 i H ρ t 1 i H s ρ s t ρ B ρ B 1 i Tr B H int ρ t . E36

We eliminate ρ t = Δ ρ t + ρ s t ρ B and collect all terms with Δ ρ t on the left-hand side. We can assume that H int B = Tr B H int ρ B = 0 because the heat bath do not exert external forces on the system (if not, replace H s by H s + H int B and H int by H int H int B ) so that also Tr B H int ρ B ρ s t = 0 and the last term ρ B 1 i Tr B H int ρ B ρ s t + ρ B ρ s t 1 i Tr B ρ B H int vanishes. Also, the term 1 i H B ρ s t ρ B disappears since H B ρ B = 0 .

We obtain

t + ε Δ ρ t 1 i H s + H int + H B Δ ρ t + ρ B 1 i Tr B H int Δ ρ t = 1 i H int ρ s t ρ B . E37

The deviation Δ ρ t vanishes when H int 0 . In lowest order with respect to H int , the solution is found as

Δ ρ t = t d t e ε t t e 1 i t t H s + H B 1 i H int ρ s t ρ B e 1 i t t H s + H B . E38

Inserting the solution (38) into the equation of motion of ρ s t (34), a closed equation of evolution is obtained eliminating ρ t . In the lowest (second) order with respect to the interaction considered here, memory effects are neglected. We can use the unperturbed dynamics to replace ρ s t = e 1 i t t H s ρ s t e 1 i t t H s and H int τ = e 1 i τ H s + H B H int e 1 i τ H s + H B so that after a shift of the integration variable

t ρ s t 1 i H s ρ s t = 1 2 0 d τ e ετ Tr B H int H int τ ρ s t ρ B = D ρ s t . E39

This result is described as quantum master equation in Born approximation. For higher orders of H int see [4, 5].

#### 3.1.3 Rotating wave approximation and Lindblad form

We assume that the interaction has the form

H int = α A α B α . E40

We use the interaction picture that coincides at t 0 with the Schrödinger picture,

O int t t 0 = e i H S + H B t t 0 / Oe i H S + H B t t 0 / E41

for any operator O . In particular, we denote

D int ρ S t t t 0 = e i H S + H B t t 0 / D ρ S t e i H S + H B t t 0 / , ρ S int t t t 0 = e iH S t t 0 / ρ S t e iH S t t 0 / E42

(note that H B commutes with ρ S t which is defined in the Hilbert space H S ).

Then, the dynamical evolution of the system is given by

t ρ S int t t t 0 = D int ρ S t t t 0 . E43

On the left-hand side, we cancel H B because it commutes with the system variables. The right-hand side, the influence term, has the form (note that ρ B commutes with H B )

D int ρ S t t t 0 = 1 2 t d t e ε t t Tr B H int int t t 0 [ H int int t t 0 ρ S int ( t t t 0 ) ] ρ B . E44

In zeroth order of interaction, ρ S int t t t 0 = e iH S t t 0 / ρ S t e iH S t t 0 / is not depending on t because the derivative with respect to t vanishes. This fact has already been used when in the Markov approximation ρ S t is replaced by ρ S t . This corresponds to the Heisenberg picture where the state of the system does not change with time. The time dependence of averages is attributed to the temporal changes of the observables.

To include the interaction, we characterize the dynamics of the system observable A introducing the spectral decomposition with respect to the (discrete) eigenstates ϕ n of H S . We introduce the eigen-energies E n s of the system S according to H S ϕ n = E n s ϕ n , and with exp ikx dx = 2 πδ k ,

A ω = d t e t t 0 e iH S t t 0 / A e iH S t t 0 / = A ω = 2 π nm ϕ n ϕ n A ϕ m ϕ m δ E n s E m s + ω E45

(the index α in (40) is dropped). In interaction picture ( A commutes with the bath observables) we have e iH S t t 0 / A e iH S t t 0 / = d ω / 2 π exp t t 0 × A ω . Now, we find for the influence term

D int ρ S t t t 0 = 1 2 t d t d ω 2 π d ω 2 π e ε t t e i ω t t e i ω + ω t t 0 × B t t B B [ A ω ρ S int ( t t t 0 ) A ω ] + BB t t B [ A ω ρ S int ( t t t 0 ) A ω ] E46

with the time-dependent bath operators B t t = exp iH B t t / B exp iH B t t / .

We can perform the integral over t that concerns the bath observables. The bath enters via equilibrium auto-correlation functions of the time-dependent bath operators B α τ . We introduce the Laplace transform of the bath correlation function (the response function of the bath)

Γ ω = 0 e i ω + i ϵ τ / Tr B ρ B B τ B = 1 2 γ ω + i 1 S ω E47

that is a matrix Γ αβ ω if the observable B has several components. We find in short notation

D int ρ S t t t 0 = d ω 2 π d ω 2 π e i ω ω t t 0 × Γ 2 ω A ω ρ S int ( t t t 0 ) A ω + Γ 1 ω A ω ρ S int ( t t t 0 ) A ω E48

after the transformation ω ω and using Eq. (45). Note that this expression for the influence term is real because the second contribution is the Hermitean conjugated of the first contribution. Using symmetry properties, all correlation functions of bath variables are related to Γ ω .

The expression ρ S int t t t 0 = e iH S t t 0 / ρ S t e iH S t t 0 / is not depending on time t because in the Heisenberg picture (we consider the lowest order of interaction) the state of the system does not depend on time. Oscillations with e i ω ω t t 0 occur that vanish for ω = ω . The rotating wave approximation (RWA) takes into account only contributions with ω = ω that are not depending on t 0 . Oscillations with e i ω ω t t 0 , ω ω 0 exhibit a phase, depending on t 0 . Any process of dephasing will damp down these oscillations.

In the case of a discrete spectrum, the spectral function (45) can be used, and the integrals over ω , ω are replaced by sums over the eigenstates ϕ n of the system S :

D int ρ S t t t 0 = 1 2 n n , m m e E n s E n s E m s + E m s t t 0 / Γ E n s E m s / × ϕ n ϕ n A ϕ m ϕ m e E m s E m s t t 0 / ρ S t ϕ m ϕ m A ϕ n ϕ n ϕ m ϕ m A ϕ n ϕ n ϕ n ϕ n A ϕ m ϕ m e E m s E m s t t 0 / ρ S t + h . c . E49

The rotating wave approximation means that n = n , m = m so that

D int ρ S t t t 0 = 1 2 n , m Γ E n s E m s / × ϕ n ϕ n A ϕ m ϕ m ρ S t ϕ m ϕ m A ϕ n ϕ n ϕ m ϕ m A ϕ n ϕ n A ϕ m ϕ m ρ S t + h . c . E50

The generalization to a more complex coupling to a bath (40) is straightforward, leading to matrices. More difficult is the discussion if the spectral function A ω is continuous, see [5]. Going back to the Schrödinger picture we have

D ρ S t = αβ Γ αβ ω A β ω ρ S t A α ω A α ω A β ω ρ S t + h . c . E51

The influence term D ρ S t cannot be given in the form of a commutator of an effective Hamiltonian with ρ S t that characterizes the Hamiltonian dynamics. Only a part can be separated that contributes to the reversible Hamiltonian dynamics, whereas the remaining part describes irreversible evolution in time and is denoted as dissipator D ρ S t .

With Γ αβ ω = γ αβ ω / 2 + iS αβ ω , we introduce the Hermitian operator H infl = αβ S αβ ω A α ω A β ω and obtain the quantum master equation

t ρ S t 1 i H S ρ S t 1 i H infl ρ S t = D ' ρ S t . E52

The dissipator has the form

D ρ S t = αβ γ αβ ω A β ω ρ S t A α ω 1 2 { A α ω A β ω ρ S t } E53

where A B = AB + BA denotes the anticommutator. The influence Hamiltonian H infl commutes with the system Hamiltonian, H S H infl = 0 , because the operator A α ω A β ω commutes with H S . It is often called the Lamb shift Hamiltonian since it leads to a shift of the unperturbed energy levels influenced by the coupling of the system to the reservoir, similar to the Lamb shift in QED. The Lindblad form follows by diagonalization of the matrices γ αβ ω ,

D ρ S t = k γ k A k ρ S t A k 1 2 { A k A k ρ S t } . E54

#### 3.1.4 Example: harmonic oscillator in a bath

A typical example is the absorption or emission of light. An isolated atom (e.g., hydrogen) is usually treated with the Schrödinger equation which gives the well-known energy eigenvalues and the corresponding eigenstates. However, this is not correct, and the finite (natural) linewidth indicate that the energetically sharp eigenstates have not an infinite life-time. The coupling to the environment, the electromagnetic field (even in the vacuum at T = 0 ) leads to transitions and a finite life-time. The electromagnetic field which is considered as bath can be represented as a system of harmonic oscillators (for each mode of the field), and the interaction with the atomic system is (dipole approximation, dipole moment D = e r )

H int = e r E = D E . E55

We discuss this phenomenon of radiation in a simplified version [5]. We consider a one-dimensional harmonic oscillator with the eigen-frequency ω 0 ,

H S = 1 2 m p 2 + m ω 0 2 2 x 2 = ω 0 a a + 1 2 , E56

with the creation a = m ω 0 / 2 1 / 2 x i / 2 m ω 0 1 / 2 p and destruction operator a = m ω 0 / 2 1 / 2 x + i / 2 m ω 0 1 / 2 p ( a a = 1 ). The discrete eigenstates ϕ n of H S are the well-known harmonic oscillator states, with eigen-energies E n s = ω 0 n + 1 / 2 . The matrix elements of the construction operators are ϕ n a ϕ n = n δ n 1 , n and its adjoint complex. In interaction picture, the equations of motion are d a t / d t = i ω 0 a t , d a t / d t = i ω 0 a t . The spectral representation reads

a ω = 2 π n n + 1 ϕ n + 1 ϕ n | δ ω + ω 0 , a ω = 2 π n n | ϕ n 1 ϕ n δ ω ω 0 . E57

At this moment, we do not specify the bath any more in detail. Suppose we have the solutions n of the energy eigenvalue problem H B m = E B , m m , then we can construct the statistical operator for the canonical distribution as

ρ B , m m 0 = m ρ B m = δ m m 1 Z e E B , m / k B T , Z = m e E B , m / k B T . E58

We introduce a weak coupling between the system and the bath

H int = exE = λ a + a B , E59

where the operator B acts only on the variables of the bath and commutes with a   and   a . In interaction picture we have

H int int t t 0 = λ a e i ω 0 t t 0 + ae i ω 0 t t 0 B t t 0 . E60

The influence term is calculated as given above. With the response function of the bath Γ ω (47) we find

t ρ S t 1 i H S ρ S t 1 i [ ( S ω 0 a a + S ω 0 aa , ρ S t ] = γ ω 0 a ρ S t a 1 2 a a ρ S t + γ ω 0 a ρ S t a 1 2 aa ρ S t . E61

The curly brackets in the dissipator denote the anticommutator. There are eight additional terms containing aa or a a . In interaction picture, they are proportional to e ± 2 i ω 0 t t 0 and are dropped within the rotating wave approximation. For a bath in thermal equilibrium, using eigenstates the detailed balance relation is easily proven,

γ ω 0 = γ ω 0 e ω 0 / k B T . E62

The evolution equations for the averages a t = Tr S ρ S a , a a t = Tr S ρ S a a are immediately calculated as

d dt a t = Tr S t ρ S t a = i ω 0 1 2 γ ω 0 γ ω 0 a t E63

with the renormalized frequency ω 0 = ω 0 + S ω 0 + S ω 0 / . The solution is

a t = a t 0 e i ω 0 γ ω 0 / 2 + γ ω 0 / 2 t t 0 . E64

Similar expressions are obtained for a t . We find for the occupation number n t = a a t = p n t

d dt a a t = γ ω 0 γ ω 0 γ ω 0 a a t E65

with the solution

a a t = a a t 0 e γ ω 0 γ ω 0 t t 0 + γ ω 0 γ ω 0 γ ω 0 1 e γ ω 0 γ ω 0 t t 0 . E66

The asymptotic behavior t t 0 is determined by the properties of the bath,

γ ω 0 γ ω 0 γ ω 0 = 1 e ω 0 / k B T 1 = n B ω 0 , E67

the system relaxes to the thermal equilibrium distribution that is independent on the initial distribution a a t 0 .

#### 3.1.5 Electromagnetic field

As example for the response function of the bath, we give the result for the blackbody radiation (Maxwell field)

Γ ij ω = 0 e i ω + τ E i τ E j 0 B = δ ij 1 2 γ ω + i S ω E68

with

γ ω = 4 ω 3 3 c 3 1 + n B ω , S ω = 2 3 π c 3 P 0 d ω k ω k 3 1 + n B ω k ω ω k + n B ω k ω + ω k . E69

Note that the Planck distribution satisfies n B ω = 1 + n B ω such that γ ω = 4 ω 3 1 + n B ω / 3 c 3 for ω > 0 and γ ω = 4 ω 3 n B ω / 3 c 3 for ω < 0 .

The resulting quantum optical master equation which, e.g., describes the coupling of atoms to the radiation field H int = D E in dipole approximation,

t ρ S t 1 i H S ρ S t 1 i H infl ρ S t = D ρ S t , E70

has the Lindblad form. The influence Hamiltonian H infl = S ω D ω D ω leads to a renormalization of the system Hamiltonian H S that is induced by the vacuum fluctuations of the radiation field (Lamb shift) and by the thermally induced processes (Stark shift). The dissipator of the quantum master equation reads

D ρ S t = 0 4 ω 3 3 c 3 1 + n B ω D ω ρ S t D ω 1 2 { D ω D ω ρ S t } + 0 4 ω 3 3 c 3 n B ω D ω ρ S t D ω 1 2 { D ω D ω ρ S t } , E71

where the integral over the negative frequencies has been transformed into positive frequencies. This result can be interpreted in a simple way. The application of the destruction operator D ω on a state of the system lowers its energy by the amount ω and describes the emission of a photon. The transition rate 4 ω 3 3 c 3 1 + n B ω contains the spontaneous emission as well as the thermal emission of photons. The term D ω gives the creation of excitations with transition rate 4 ω 3 3 c 3 n B ω describing the absorption of photons.

#### 3.1.6 The Pauli equation

We consider a system whose state is described by the observable A , and which takes the value a . This can be a set of numbers in the classical case that describe the degrees of freedom we use as relevant variables. In the quantum case, this is a set of relevant observables that describe the state of the system. The eigenvalue a corresponds to a state vector a in the Hilbert space.

At time t , we expect a probability distribution p 1 a t to find the system in state a , if the property A is measured. The change of the probability p 1 a t with time is described by a master equation or balance equation

d dt p 1 a t = a a w a a p 1 a t w a a p 1 ( a t ) . E72

In the context of the time evolution of a physical system, this master equation is also denoted as Pauli equation. We derive it from a microscopical approach using perturbation theory. The statistical operator ρ t follows the von Neumann equation of motion (8) with the Hamiltonian

H = H 0 + λ H E73

where the solution of the eigenvalue problem for H 0 is known, H 0 n = E n n . The probabilities to find the system in the state n are given by the diagonal elements of ρ t in this representation,

p 1 n t = n ρ t n . E74

First, we consider the special case λ = 0 , where the von Neumann equation is easily solved:

ρ nm t = n ρ t m = e i ω nm t t 0 ρ nm t 0 , ω nm = E n E m E75

if ρ nm t 0 is given. The nondiagonal elements ρ nm t , n m are oscillating. The periodic time dependence of the density matrix that arises in the nondiagonal elements has nothing to do with any time evolution or irreversibility. It expresses the coherences in the system. The diagonal elements

ρ nn t = p 1 n t = n ρ t n E76

do not change with time and can be considered as conserved quantities if λ = 0 .

To find the initial distribution, we consider the probabilities as relevant observables that describe the nonequilibrium state at t 0 . If there are no further information on coherence, the relevant statistical operator is diagonal,

ρ rel t 0 = n p 1 n t 0 n n = n p 1 n t 0 P n . E77

We introduced the projection operator P n = n n . The solution is ρ t = ρ rel t 0 . The case λ = 0 is a trivial case, nothing happens.

Now, we consider a small perturbation as expressed by the parameter λ . As before, we consider the probabilities as relevant observables that describe the system in nonequilibrium. We project the diagonal part of the statistical operator,

ρ rel t = diag ρ t = D n ρ t = n P n ρ t P n . E78

The difference ρ irrel t = ρ t ρ rel t = 1 D n ρ is the irrelevant part of the full statistical operator.

The problem to obtain the time evolution of the probabilities p 1 n t is solved if we find an equation of evolution for ρ rel t . We use the method of the nonequilibrium statistical operator and start with the extended von Neumann equation (27). For the projection, we obtain ( D n is linear and commutes with / t )

t ρ rel t = 1 i D n λ H ρ irrel t . E79

We assumed that H 0 is diagonal with ρ rel t so that the commutator vanishes. Furthermore, the diagonal elements of the commutator of a diagonal matrix with an arbitrary matrix disappear. For the irrelevant part we have

t ρ irrel t + ϵ ρ irrel t 1 i 1 D n H ρ irrel t = 1 i 1 D n λ H ρ rel t . E80

On the right-hand side, we can drop the projector D n . Its action disappears because ρ rel is diagonal. It is seen that ρ irrel t is of the order λ .

In the remaining projection 1 D n H 0 ρ irrel t + 1 D n H ρ irrel t , the second contribution is of second order in λ and will be dropped here because we consider only the lowest order in λ ( ρ irrel t is also of the order λ ). This is denoted as Born approximation. We have

t ρ irrel t + ερ irrel t 1 i H 0 ρ irrel t = 1 i λ H ρ rel t . E81

The solution is simple by integration,

ρ irrel t = 1 i t e ε t 1 t e i H 0 t 1 t λ H ρ rel t 1 e i H 0 t 1 t dt 1 . E82

The proof is given by insertion.

With this expression for ρ irrel t , we find a closed equation for ρ rel t ,

t ρ rel t = λ 2 2 D n t e ε t 1 t H e i H 0 t 1 t H ρ rel t 1 e i H 0 t 1 t dt 1 . E83

This result describes a memory effect. The change of ρ rel t is determined by the values ρ rel t 1 at all previous times t 1 t . In the Markov approximation, we replace ρ rel t 1 by ρ rel t so that memory effects are neglected. This is justified in the limit λ 0 because then ρ rel t changes only slowly with time. Then

t ρ rel t = λ 2 2 D n t e ε t 1 t H e i H 0 t 1 t H e i H 0 t 1 t ρ rel t dt 1 . E84

This expression has similar structure as the QME (39) an can be treated in the same way. The right-hand side D ρ rel t is related to the dissipator after subtracting the Lamb shift contribution.

Explicit expressions for the time evolution of the density matrix are obtained by projection on the basis n . With the matrix elements n ρ rel t m = δ n , m p 1 n t as well as n H 0 m = δ n , m E n and n H m = H nm we have

d dt p 1 n t = λ 2 2 m H nm H mn p 1 n t p 1 m t × t e ε t 1 t e i E m E n t 1 t + e i E m E n t 1 t dt 1 . E85

Performing the integral over t 1 , we find [with the Dirac identity lim ϵ + 0 1 x + i ϵ P 1 x iπδ x ] the Pauli equation

d dt p 1 n t = n n w n n p 1 n t w n n p 1 ( n t ) . E86

The transition rates are given by Fermi’s Golden rule,

w nm = lim ϵ 0 λ 2 2 H nm 2 1 i ω nm + ϵ + 1 i ω nm + ϵ = 2 πλ 2 H nm 2 δ E n E m . E87

#### 3.1.7 Properties of the Pauli equation

The transition rate w nm obeys the condition of detailed balance, w mn = w nm , the inverse transition has the same rate. This follows because H is hermitean,

n H m = m H ′+ n = m H n . E88

An important property is that it describes irreversible evolution with time. For the relevant entropy S r el t = k B n p 1 n t ln p 1 n t we find

dS rel t dt = k B n m w nm p 1 m t p 1 n t ln p 1 n t k B n p 1 n t p 1 n t p 1 n t t = 1 2 k B n m w nm p 1 n t p 1 m t ln p 1 n t ln p 1 m t 0 . E89

We used d dt n p 1 n t = d dt 1 = 0 and interchanged n with m in the half of the expression. Since ln x is a monotonic function of x , the relation x 1 x 2 ln x 1 ln x 2 0 holds. Considering states n , m where transitions are possible, equilibrium ( dS rel t / dt = 0 ) occurs if p 1 m t = p 1 n t ; else S rel t increases with time. Equipartition corresponds to the microcanonical ensemble in equilibrium.

#### 3.1.8 Example: transition rates

We consider transitions between eigenstates of H 0 owing to interaction. A typical case is the collisions expressed by a k 1 a k 2 a k 2 a k 1 between the (momentum) eigenstates k of H 0 . This is discussed in the following section on kinetic theory. Another example is minimal coupling known from QFT between a Dirac fermionic field (electron) and the Maxwell bosonic field (photons), with ( E k = 2 k 2 / 2 m , ω q = c q )

H 0 = k E k a k a k + q ω q b q b q E90

(spin and polarization variables are not indicated separately), and the interaction

H int = k , k , q v k k q a k a k b q + h . c . E91

The transition rates (87) are calculated between the initial state n = k , energy E n = E k , and the final state m = k , q , energy E m = E k + ω q for emission in the vacuum state. For absorption, the corresponding process can be given. For free particles k = k , σ , the matrix element v k σ k σ q δ k ' + q , k must fulfill momentum conservation. Together with the conservation of energy in Eq. (87), the second-order transition rate vanishes. Only in fourth order, different contributions (Compton scattering, pair creation) are possible. If considering an radiating atom, the electrons are moving in orbits around the nucleus, k = nlm , σ . Momentum conservation is not required, and the standard expressions (Fermi’s Golden rule) for absorption and emission of light by an atom are obtained. The corresponding rate equation (86) describes natural line width, detailed balance, and thermal equilibrium as stationary solution.

#### 3.1.9 Conclusions

Quantum master equations and the Pauli equation are fundamental expressions to describe nonequilibrium phenomena, such as one-step processes of excitation and deexcitation, two-level systems, nuclear decay, chemical reactions, and also conductivity where electrons are scattered by ions, etc. A basic assumption is the subdivision into a system and a bath. In Born-Markov approximation, correlations between system and bath (back-reactions) are neglected. Projection to diagonal elements of the reduced density matrix or the Rotating wave approximation lead to irreversible equations of evolution (dissipator) as derived by Zwanzig, Lindblad, Kossakowski, and others. Further developments of the theory are, e.g., the Nakajima-Zwanzig equation or the Quantum Fokker-Planck equation [4]. A fundamental problem is the subdivision in relevant (system) and irrelevant (bath) degrees of freedom. If correlations between the system and bath become relevant, the corresponding degrees of freedom of the bath must be included in the set of system variables.

### 3.2 Kinetic theory

Historically, nonequilibrium statistical physics was first developed as the kinetic theory of gases [7] by Boltzmann. We start with classical systems to explain the problem to be solved in kinetic theory. The more general case of quantum systems contains no additional complications, but the concepts become more evident in the classical limit. We give results for both cases, the general quantum case and the classical limit. Reduced distribution functions are considered as the relevant observables. Closed equations of evolution are obtained describing irreversible processes.

#### 3.2.1 The Liouville equation

The standard treatment of a classical dynamical system can be given in terms of the Hamilton canonical equations. In classical mechanics, we have generalized coordinates and canonic conjugated momenta describing the state of the system, e.g., a point in the 6 N -dimensional phase space ( Γ -space) in the case of N point masses. The 6 N degrees of freedom r 1 p 1 r N p N define the microstate of the system. The evolution of a particular system with time is given by a trajectory in the phase space. Depending on the initial conditions different trajectories are taken.

Within statistical physics, instead of a special system, an ensemble of identical systems is considered, consisting of the same constituents and described by the same Hamiltonian, but at different initial conditions (microstates), which are compatible with the values of a given set of relevant observables characterizing the macrostate of the system. The probability of the realization of a macrostate by a special microstate, i.e., a point in the 6 N -dimensional phase space ( Γ -space), is given by the N -particle distribution function f N r i p i t which is normalized,

f N r i p i t = 1 ; d Γ = d N r d N p N ! h 3 N = d 3 N x d 3 N p N ! h 3 N . E92

In nonequilibrium, the N -particle distribution function depends on the time t .

The macroscopic properties can be evaluated as averages of the microscopic quantities a r i p i with respect to the distribution function f N r i p i t :

A t = d Γ a r i p i f N r i p i t . E93

In addition to these so-called mechanical properties there exist also thermal properties, such as entropy, temperature, and chemical potential. Instead of a dynamical variable, they are related to the distribution function. For example, the equilibrium entropy is given by

S eq = k B f N r i p i t ln f N r i p i t E94

We derive an equation of motion for the distribution function f N r i p i t , the Liouville equation, see [5]:

d f N d t = f N t + i = 1 N f N r i r ̇ i + f N p i p ̇ i = 0 . E95

We shortly remember the quantum case. Instead of the N -particle distribution function f N t , the statistical operator ρ t is used to indicate the probability of a microstate in a given macrostate. The equation of motion is the von Neumann equation (8). Both equations are closely related and denoted as Liouville-von Neumann equation.

#### 3.2.2 Classical reduced distribution functions

To evaluate averages, instead of the N -particle distribution function f N r 1 r N p 1 p N t often reduced s -particle distribution functions

f s r 1 p s t = d 3 r s + 1 d 3 p N N s ! h 3 N s f N r 1 p N t E96

are sufficient. Examples are the particle density, the Maxwell distribution of the particle velocities, and the pair correlation function.

We are interested in the equations of motion for the reduced distribution functions. For classical systems, one finds a hierarchy of equations. From the Liouville equation, Eq. (95) without external potential,

d f N d t = f N t + i N v i f N r i i j N V ij r i f N p i = 0 E97

we obtain the equation of motion for the reduced distribution function f s through integration over the 3 N s other variables:

d f s d t = f s t + i = 1 s v i f s r i i j s V ij r i f s p i = i = 1 s d 3 r s + 1 d 3 p s + 1 h 3 V i , s + 1 r i f s + 1 r 1 p s + 1 t p i . E98

This hierarchy of equations is called BBGKY hierarchy, standing for Bogoliubov, Born, Green, Kirkwood, and Young.

The equation of motion (98) for the reduced distribution function f s is not closed because on the right-hand side the higher order distribution function f s + 1 appears. In its turn, f s + 1 obeys a similar equation that contains f s + 2 , etc. This structure of a system of equations is denoted as hierarchy. To obtain a kinetic equation that is a closed equation for the reduced distribution function, one has to truncate the BBGKY hierarchy, expressing the higher order distribution function f s + 1 by the lower order distribution functions f 1 f s .

#### 3.2.3 Quantum statistical reduced distributions

In the quantum case, the distribution function f N is replaced by the statistical operator ρ that describes the state of the system, and the equation of motion is the von Neumann equation (8). The quantum statistical reduced density matrix is defined as average over creation and annihilation operators,

ρ s r 1 r s t = Tr ρ t ψ r 1 ψ r s ψ r s ψ r 1 . E99

It is related to correlation functions, the Wigner function, Green functions, dynamical structure factor, and others.

We consider the equations of motion for reduced distribution functions. For the single-particle density matrix in momentum representation, we have

ρ 1 p p t = Tr ρ t ψ p ψ p . E100

Derivation with respect to time gives

t ρ 1 p p t = 1 i Tr H ρ ψ p ψ p = 1 i Tr ρ ψ p ψ p H . E101

Similar as for the BBGKY hierarchy, we obtain in general a hierarchy of equations of the form

ρ s t t = function of ρ s t ρ s + 1 t . E102

Like in the classical case, we have to truncate this chain of equations. For example, in the Boltzmann equation for f 1 t , the higher order distribution function f 2 t is replaced by a product of single-particle distribution functions f 1 t .

#### 3.2.4 Stoßzahlansatz and Boltzmann equation

To evaluate the averages of single-particle properties such as particle current or kinetic energy, only the single-particle distribution must be known. Then, the single-particle distribution contains the relevant information, the higher distributions are irrelevant and will be integrated over.

We are looking for an equation of motion for the single-particle distribution function f 1 r p t , taking into account short range interactions and binary collisions. For the total derivative with respect to time we find, see Eq. (95),

d f 1 d t = t f 1 + r ̇ r f 1 + p ̇ p f 1 = t f 1 + v r f 1 + F p f 1 = 0 .

The crucial point in this equation is the force F . It is the sum of external forces F ext acting on the system under consideration and all forces resulting from the interaction V ij r i r j between the constituents of the system.

Before discussing the derivation of kinetic equations using the method of the nonequilibrium statistical operator, we give a phenomenological approach using empirical arguments. To describe the change in the distribution function f 1 due to collisions among particles, we write

t f 1 = t f 1 D + t f 1 St , E103

where the drift term contains the external force,

t f 1 D = v r f 1 F ext p f 1 E104

and the internal interactions are contained in the collision term t f 1 St for which, from the BBGKY hierarchy (98), an exact expression has already been given:

t f 1 St = d 3 r d 3 p h 3 V r r r p f 2 rp r p t . E105

Collisions or interactions among particles occur due to the interaction potential V r r , which depends on the coordinates of the two colliding partners. For every particle, one has to sum over collision with all partners in the system. In this way, we have an equation for the single-particle distribution function, but it is not closed because the right-hand side contains the two-particle distribution function f 2 rp r p t .

As an approximation, similar to the master equation, we assume a balance between gain and loss:

f 1 t St = G L . E106

With some phenomenological considerations [5], we can find the collision term as

f 1 t St = d 3 v 2 d σ v 1 v 2 f 1 r v 1 t f 1 r v 2 t f 1 r v 1 t f 1 r v 2 t , E107

where we have introduced the differential cross section

d σ d Ω = b ϑ sin ϑ d b ϑ d ϑ . E108

Inserting expression (108) into Eq. (103), we obtain a kinetic equation only for the single-particle distribution, the Boltzmann equation.

#### 3.2.5 Derivation of the Boltzmann equation from the nonequilibrium statistical operator

The relevant observable to describe the nonequilibrium state of the system is the single-particle distribution function. First, we consider classical mechanics where the single-particle distribution function is f 1 r p t .

We can write the single-particle distribution as an average (93) of a microscopic (dynamic) variable, the single-particle density

f 1 r p t = n 1 r 1 p N r p t , n 1 r 1 p 1 r N p N r p = 3 i = 1 N δ 3 r r i δ 3 p p i . E109

The self-consistency conditions (18) are realized with the Lagrange parameter F 1 r p t . The relevant distribution F rel reads (see (19) and replace n by d 3 rd 3 p / h 3 )

F rel r 1 p N t = exp Φ t i = 1 N F 1 ( r i p i t ) , Φ t = ln exp i = 1 N F 1 r i p i t . E110

The constraints f 1 r p t F rel r 1 p N t n 1 r 1 p N r p are solved according to

f 1 r p t = h 3 N e F 1 r p t e F 1 r p t d 3 r d 3 p 1 , F 1 r p t = ln f 1 r p t . E111

This means, we can eliminate the Lagrange parameters F 1 r p t that are expressed in terms of the given distribution function f 1 r p t . The relevant distribution is

F rel r 1 p N t = 1 Z rel j f 1 r j p j t , Z rel = j f 1 r j p j t d Γ N = N N N ! e N . E112

The Boltzmann entropy is then

S rel t = k B ln F rel t = k B f 1 r p t ln f 1 r p t e d 3 r d 3 p h 3 . E113

Below, we show that it increases with time for nonequilibrium distributions.

The relevant distribution can be used to derive the collision term (107), for details see [3]. We will switch over to the quantum case where the presentation is more transparent.

In the quantum case, we consider the single-particle density matrix. In the case of a homogeneous system ( n 1 r = n ), ρ 1 p p is diagonal. The set of relevant observables are the occupation number operators n p ,

n p t = f 1 p t . E114

Considering these mean values as given, we construct the relevant statistical operator as

ρ rel t = e Φ t p F 1 p t n p , Φ t = ln Tr e p F 1 p t n p . E115

The Lagrange parameters F 1 p t are obtained from the self-consistency conditions (114) similar to Eq. (111);

f 1 p t = Tr e p F 1 p t n p n p Tr e p F 1 p t n p = i n i e F 1 p i t n i 1 + δ p i , p n i 1 i n i e F 1 p i t n i E116

so that

f 1 p t = 1 e F 1 p t ± 1 + : Fermions : Bosons , F 1 p t = ln 1 f 1 p t ln f 1 p t . E117

As in the classical case, also in the quantum case, the Lagrange parameters can be eliminated explicitly.

We now derive the Boltzmann equation for the quantum case, see [3]. With the statistical operator (Eq. (25) after integration by parts)

ρ t = ρ rel t t e ϵ t 1 t d d t 1 e i H ( t t 1 ρ rel t 1 e i H ( t t 1 d t 1 , E118

With n ̇ p = i H n p , we get the time derivative of the single-particle distribution function

t f 1 p t = Tr ρ rel t n ̇ p 0 e ϵ t Tr n ̇ p d d t e i H t ρ rel t + t e i H t d t . E119

Because the trace is invariant with respect to cyclic permutations and ρ rel t commutes with n p , see (115),

Tr ρ rel t n ̇ p = i Tr ρ rel H n p = i Tr H n p ρ rel = 0 , E120

and Eq. (119) can be written as

f 1 t = 1 2 0 d t e ϵ t T r H n p e i H t [ H ρ rel ] e i H t , E121

if we neglect the explicit time dependence of ρ rel t (no memory effects, the collision term is local in space and time). Next, we introduce two more integrations via delta functions to get rid of the time dependence in the trace:

f 1 t = 1 2 d E d E 0 d t e ϵ + i E E t Tr V n p δ E H V ρ rel δ E H . E122

(We take into account that the kinetic energy in H commutes with n p so that only the potential energy V remains.) This equation can be expressed by so-called T matrices, T = V + V 1 E H T ,

f 1 t = π d E Tr T n p δ E H 0 T ρ rel δ E H 0 , E123

For further treatment, we choose the approximation of binary collisions, that means that only two particles change their momentums during a collision. In second quantization, the T matrix is then

T p 1 , p 2 , p 1 p 2 a p 1 a p 2 t p 1 p 2 p 1 p 2 a p 2 a p 1 δ p 1 + p 2 p 1 p 2 , E124

with the two-particle T matrix t p 1 p 2 p 1 p 2 . With this T matrix, we find the collision term (time t is dropped)

f 1 p 1 t St = p 2 p 1 p 2 w p 1 p 2 p 1 p 2 f 1 p 1 f 1 p 2 1 f 1 p 1 1 f 1 p 2 f 1 p 1 f 1 p 2 1 f 1 p 1 1 f 1 p 2 E125

with the transition probability rate

w p 1 p 2 p 1 p 2 = 2 π t p 1 p 2 p 1 p 2 t p 1 p 2 p 2 p 1 2 δ E p 1 + E p 2 E p 1 E p 2 δ p 1 + p 2 p 1 p 2 , E126

which leads to the quantum statistical Boltzmann equation.

#### 3.2.6 Properties of the Boltzmann equation

The Boltzmann equation is a nonlinear integro-differential equation for the single-particle distribution function in the classical case. In the quantum case, we can use the density matrix or the Wigner function to characterize the nonequilibrium state of the system. The Boltzmann equation is valid in low-density limit (only binary collisions). At higher densities also three-body collisions, etc., must be taken into account. Further density effects such as the formation of quasi particles and bound states have to be considered. The collision term is approximated to be local in space and time, no gradients in the density and no memory in time is considered. The assumption of molecular chaos means that correlations are neglected, the two-particle distribution function is replaced by the product of single-particle distribution functions.

The increase of entropy (Boltzmann H theorem) can be proven. In terms of the relevant statistical operator, the entropy is

S rel = k B p 1 + f 1 p ln 1 f 1 p f 1 p ln f 1 p . E127

The change with time follows from

d S rel d t = k B p f 1 t ln f 1 k B p f 1 t + k B p f 1 t ln 1 f 1 + k B p f 1 t = k B p 1 p 2 p 1 p 2 w p 1 p 2 p 1 p 2 ln 1 f 1 p 1 1 1 f 1 p 1 1 1 f 1 p 2 1 1 f 1 p 1 1 1 f 1 p 2 1 × f 1 p 1 f 1 p 1 f 1 p 2 f 1 p 2 . E128

We interchange indices 1 2 , 1 2 ; furthermore 1 1 , 2 2 ; and 1 2 , 2 1 , use the symmetries of w p 1 p 2 p 1 p 2 and x 1 x 2 ln x 1 ln x 2 0 because ln x is a monotonous function of x . We obtain 4 d S rel d t 0 , the Boltzmann (relevant) entropy can increase.

The collision integral guarantees conservation of total momentum, particle number, and kinetic energy. However, the total energy including the interaction part is not conserved. The equilibrium solution f 1 0 p follows from d S rel d t = 0 :

1 f 1 0 p 1 1 f 1 0 p 1 1 1 f 1 0 p 1 1 f 1 0 p 1 1 = 0 . E129

If f 1 0 p depends only on energy, we find the well-known result for ideal quantum gases,

1 f 1 0 p 1 = e β E p μ , f 1 0 p = e β E p μ ± 1 1 . E130

In the classical limit, we have f 1 0 p = e β E p μ with e βμ = N Ω 2 π 2 mk B T 3 / 2 1 2 s + 1 , where s denotes the spin of the particle.

#### 3.2.7 Beyond the Boltzmann kinetic equation

In deriving the Boltzmann equation, different approximations have been performed: only binary collisions are considered, three-particle, and higher order collisions are neglected. Memory effects and spatial inhomogeneities have been neglected. The single-particle distribution was considered as relevant observable in the Markov approximation. These approximations can be compared with the Born-Markov approximation discussed in context with the quantum master equation. Instead of the Born approximation that is possible for weak interactions, the binary collision approximation is possible in the low-density limit, where three- and higher order collisions are improbable.

In the case of thermal equilibrium, the Boltzmann entropy S r el (127) coincides with the entropy of the ideal (classical or quantum) gas. The equilibrium solution of the Boltzmann equation leads to the entropy of the ideal gas and gives not the correct equation of state for an interacting system that are derived from the Gibbs entropy ( Φ = ln Z is the Matthieu-Planck function)

S eq = k B Φ + βH exp Φ βH , E131

see Eq. (13). This deficit of the Boltzmann equation arises because binary collisions are considered where the kinetic energy of the asymptotic states is conserved. Only the single-particle distribution is a relevant observable and is correctly reproduced. It can be improved if the total energy, which is conserved, is considered as a relevant observable. Alternatively, we can also include the two-particle distribution function in the set of relevant observables. An important example is the formation of bound states as a signature of strong correlations in the system. Then, the momentum distribution of bound states has to be included in the set of relevant observables.

#### 3.2.8 The linearized Boltzmann equation

Different approximations are known to obtain solutions of the Boltzmann equation, see [4, 5]. A serious problem in solving the Boltzmann equation is its nonlinearity as we have terms of the form f 1 p 1 t f 1 p 2 t . Special cases that allow for linearization are two-component systems with a large difference in the masses or concentration. Linearization is also possible in the case where the deviation from some equilibrium distribution is small. As an application, we consider the calculation of electrical conductivity in plasmas.

We investigate a plasma of ions and electrons under the influence of an external electric field E ext . For simplicity, we assume E ext to be homogeneous and independent of time (statical conductivity σ ). For moderate fields, we await a linear behavior of the plasma following Ohm’s Law:

j el = σ E . E132

[Note that in Eq. (132) E is not the external field, but the effective electric field in the medium (the plasma), being the superposition of the external field E ext and the polarization field ε P ]. j el is the average electric current defined via the single-particle distribution function f 1

j e l = 1 Ω i N e i v i = s e s d 3 v v f 1 v s = s e s m s d 3 p 2 π 3 p f 1 p s . E133

Here, we have kept the index s for the different sorts. In the following, we will skip this index as we only consider electrons being responsible for the electric current.

We recall the Boltzmann equation

p m r f 1 + e E p f 1 + t f 1 St = 0 , E134

m is the electron mass and e the electron charge. The first term in this equation vanishes because of homogeneity of the system. For the collision term, we take the expression Eq. (125) in the generalized form for quantum systems. After the distribution function of the collision partner has been replaced by the equilibrium distribution, we have

t f 1 St = d 3 p Ω 2 π 3 f 1 p w p p 1 f 1 p f 1 p w p p 1 f 1 p , E135

where w p p is the transition rate from the momentum state p to the state p . The quantum behavior of the collisions is taken into account via the Pauli blocking factors 1 f 1 p .

#### 3.2.9 Example: conductivity of the Lorentz plasma

In the Lorentz plasma model, the electron-electron collisions are neglected, and only electron-ion collisions are considered, interaction potential V ei r . In the adiabatic approximation where the ions are regarded as fixed at positions R i (elastic collisions), the interaction part of the Hamiltonian reads

H = i V ei r R i . E136

In Born approximation (or time-dependent perturbation theory), the transition rate is given by Fermi’s Golden rule:

w p p = 2 π H p p 2 δ E p E p = w p p ; E p = p 2 / 2 m . E137

To solve the Boltzmann equation Eq. (135), we make use of the ansatz

f 1 p = f 1 0 E p + Φ p d f 1 0 E p d E p k B T = f 1 0 E p 1 + Φ p 1 f 1 0 E p . E138

For equilibrium distributions, we have the detailed balance condition

w p p f 1 0 E p 1 f 1 0 E p = w p p f 1 0 E p 1 f 1 0 E p . E139

Insertion of Eq. (138) into the Boltzmann equation Eq. (135) yields with Eq. (139)

e mk B T E p f 1 0 E p 1 f 1 0 E p = d 3 p Ω 2 π 3 w p p f 1 0 E p 1 f 1 0 E p Φ p Φ p , E140

where we have neglect terms with higher order of E and have used the fact that Φ p E . With the definition of the relaxation time tensor τ ̂ p , according to Φ p = e / mk B T E τ ̂ p p , the equation reads

e E p = d 3 p Ω 2 π 3 w p p f 1 0 E p f 1 0 E p e E τ ̂ p p τ ̂ p p , E141

e E = E / E . The electric current density Eq. (133) depends only on the deviation of the distribution function since f 1 0 is an even function in p (isotropy). We obtain by insertion of Eq. (138) into Eq. (133)

j el = e Ω 2 d 3 p Ω 2 π 3 p m Φ p f 1 0 E p 1 f 1 0 E p . E142

The conductivity σ is the proportionality factor between the current density and the effective field E :

σ = e 2 m 2 k B T 2 d 3 p 2 π 3 p z τ ̂ p p z f 1 0 E p 1 f 1 0 E p . E143

We have derived an analytical expression for the conductivity of a Lorentz plasma in terms of the relaxation time tensor τ ̂ p . For isotropic systems, τ ̂ ij = τδ ij , the well-known Ziman formula σ τ = τ ne 2 / m for the conductivity results.

The solution of Eq. (141) for a momentum-dependent relaxation time is

τ E p = d 3 p Ω 2 π 3 w p p 1 cos ϑ 1 E144

as can be verified by insertion. Now, the conductivity reads with Eq. (137)

σ = 2 e 2 m 2 k B T Ω d 3 p p z 2 f 1 0 E p 1 f 1 0 E p 2 π d 3 p H p p 2 δ E p E p 1 cos ϑ 1 . E145

Considering the screened interaction potential (Debye potential) V ei D r = e 2 4 π ϵ 0 r e κ r with the Debye screening parameter κ 2 = e 2 N / ϵ 0 k B T Ω , the evaluation can be performed. With

Λ p = 0 2 p / 1 q 2 + κ 2 2 q 3 d q = ln 1 + b 1 2 b 1 + b , b = 4 p 2 k B T Ω ϵ 0 e 2 2 N , E146

we finally obtain for the conductivity [5].

σ = 2 5 / 2 k B T 3 / 2 4 π ϵ 0 2 π 3 / 2 m 1 / 2 e 2 Λ ; Λ Λ p 2 / 2 m = 3 k B T / 2 . E147

#### 3.2.10 Conclusions

The method of the nonequilibrium statistical operator gives not only the derivation of the Boltzmann equation (quantal and classical), but indicates also possible improvements such as conservation of total energy, inclusion of bound state formation, hydrodynamic equations, etc.

The solution of the general Boltzmann equation is not simple, in addition to numerical simulations different approximations have been worked out. For the linearized Boltzmann equation, the relaxation time approximation can be used for elastic scattering, but for the general case (inclusion of electron-electron collisions in a plasma), the Kohler variational principle [11] can be applied. Landau-Vlasov equations for mean-field effects as well as Fokker-Planck equations for the collision term have been investigated.

The basic assumption to derive the Boltzmann equation is the selection of the single-particle distribution as relevant observable. Correlations are neglected and have to be built up in higher orders of approximation or extending the set of relevant observables. The most appropriate systems for kinetic theory are dilute gases where the collision time is short compared with the time of free flight. Irreversibility is owing to the Stoßzahlansatz for the intrinsic interaction.

### 3.3 Linear response theory

A third example, which allows the explicit elimination of the Lagrange multipliers to fulfill the self-consistency conditions, is a system near to thermodynamic equilibrium which is under the influence of mechanical (external forces) or thermodynamic (gradients of temperature, pressure, chemical potentials, etc.) perturbations. As response, currents appear in the system. Assuming linearity for small perturbations, transport coefficients are defined. Fluctuations in equilibrium are considered as a nonequilibrium state which relaxes to equilibrium, see Eq. (7).

#### 3.3.1 Response to an external field

We consider a system under the influence of external (time dependent) fields acting on the particles, see [4, 11, 12, 13, 14, 15, 16],

H t = H S + H F t , E148

where H S denotes the system Hamiltonian, containing all kinetic energies of the particles as well as the full interaction part. The second part H F t describes the coupling of the system to the external fields h j :

H F t = j h j e i ωt A j . E149

We characterize the nonequilibrium state by the set B n of relevant observables. In the following, we assume that the equilibrium expectation values of the nonequilibrium fluctuations disappear, B n e q = 0 (else, we have to subtract the equilibrium values).

Treating the conserved observables explicitly, we write the relevant statistical operator ρ rel in the form ( H = H S c μ c N c )

ρ rel t = e Φ t β H n F n t B n , Φ t = ln Tr e β H n F n t B n , E150

where the Lagrange multipliers are divided into the equilibrium parameters β , μ and the generalized response parameters F n t , coupled to the corresponding observables. All Lagrange parameters are determined by the given mean values of these observables. In particular, we have the self-consistency conditions (18)

B n rel t = Tr ρ rel t B n = Tr ρ t B n = B n t E151

or

Tr ρ irrel t B n = 0 , ρ irrel t = ρ t ρ rel t . E152

The corresponding self-consistency condition for N and H S lead to the well-known equations of state for the temperature 1 / β and the chemical potential μ . Φ t is the Massieu-Planck functional that normalizes ρ rel t .

We consider the limit of weak external fields. Compared with the equilibrium distribution (13), we expect that the changes of the state of the system are also weak. We characterize the nonequilibrium state by the set B n of relevant observables and assume that the averages

B n t = Tr ρ t B n h j e i ωt E153

are proportional to the external fields (linear response).

The basic assumption of LRT is that the average values B n t of the additional observables, which characterize the response of the system, are proportional to the external fields. Because these external fields are arbitrarily weak, we expand all quantities with respect to the fields up to first order. If the fluctuations B n t are proportional to these fields, we have also F n h j . Below, we derive linear equations that relate the response of the system to the causing external fields.

In the linear regime, we await the response parameters F n t to exhibit the same time dependence as the external fields:

F n t = F n e i ωt . E154

Here, we have harmonic fields h j e i ωt , but the formulation rests general as we can always express arbitrary time dependences by means of a Fourier transformation. Within the linear regime, the superposition of different components of the field gives the superposition of the corresponding responses. The treatment of spatial dependent external forces is also possible. As a specific advantage of the Zubarev method, thermodynamic forces such as gradients of temperature or chemical potentials can be treated [4, 5, 15, 16].

#### 3.3.2 Elimination of the Lagrange multipliers

The main problem is to eliminate the Lagrange multipliers, the generalized response parameters F n t . As in the case of kinetic theory, this is also possible explicitly in the case of linear response theory (LRT). With the operator relation e A + B = e A + 0 1 d λ e λ A + B B e 1 λ A , we get for the relevant statistical operator (150) up to first order of the nonequilibrium fluctuations B n

ρ rel t = ρ eq + β 0 1 d λ n F n t B n i βλ ρ eq . E155

Here, we made use of the modified-Heisenberg picture O τ = exp i H τ / O exp i H τ / with τ i βλ replacing in the exponents H S by H = H S c μ c N c . We want to calculate expectation values of macroscopic relevant variables that commute with the particle number operator N c so that we can use both H and H S synonymously. (Mention that also the Massieu-Planck functional Φ t has to be expanded so that the fluctuations around the equilibrium averages B n B n eq appear).

#### 3.3.3 Linearization of the NSO

All terms have to be evaluated in such a way, that the total expression rests of order O h . For expressions (25) and (26), we find after integration by parts

ρ ϵ t = ρ rel t t d t 1 e ϵ t 1 t U t t 1 i H S + H F t 1 ρ rel t 1 + t 1 ρ rel t 1 U t t 1 . E156

Since H S commutes with ρ eq (equilibrium!), the curly bracket is of order O h . In particular, we have for the first term the time derivative in the Heisenberg picture,

i H S β 0 1 d λ n F n t 1 B n i λβ ρ eq = β 0 1 d λ n F n t i B ̇ n i λβ ρ eq . E157

For the second term of the integral in Eq. (156), we use Kubo’s identity

B e A = 0 1 d λ e λ A B A e 1 λ A . E158

so that

i H F t 1 ρ eq = β e i ω t 1 0 1 d λ j h j A ̇ j i λβ ρ eq . E159

The last term in the curly bracket can be rewritten as

t 1 ρ rel = β 0 1 d λ n F ̇ n t 1 B n i λβ ρ eq . E160

Because we restrict ourselves to the order O h , for the time evolution operator we have U t t 1 e iH S t t 1 / .

After linearization with respect to the external fields h j and the response parameters F n , finally we have

ρ ϵ t = ρ rel t β e i ωt 0 d t 1 e i zt 1 0 1 d λ j h j A ̇ j i λβ + t 1 ρ eq + n F n B ̇ n i λβ + t 1 ρ eq i ω F n B n i λβ + t 1 ρ eq E161

( z = ω + ). Here, we used that h j t and F n t , Eq. (154), are proportional to e i ωt .

We multiply this equation by B m , take the trace and use the self-consistency relation (151). We obtain a set of linear equations for the thermodynamically conjugated parameters F n (response parameters):

n B m B ̇ n z i ω B m B n z F n = j B m A ̇ j z h j , E162

with the Kubo scalar product (the particle number commutes with the observables)

A B = 0 1 d λ Tr A e λβ H B e λβ H ρ eq = 0 1 d λ Tr A B i λβ ρ eq , E163

and its Laplace transform, the thermodynamic correlation function

A B z = 0 d t e i zt A B t = 0 d t e i zt A t B . E164

The linear system of equations (162) has the form

n P mn F n = j D mj h j E165

to determine the response parameters F n , the number of equations coincides with the number of variables to be determined. The coefficients of this linear system of equations are given by equilibrium correlation functions. We emphasize that in the classical limit the relations become more simple because the variables commute, and we have not additional integrals expanding the exponential.

We can solve this linear system of equations (162) using Cramers rule. The response parameters F n are found to be proportional to the external fields h j with coefficients that are ratios of two determinants. The matrix elements are given by equilibrium correlation functions. This way, the self-consistency conditions are solved, and the Lagrange multipliers can be eliminated. The nonequilibrium problem is formally solved. The second problem, the evaluation of equilibrium correlation functions, can be solved by different methods such as numerical simulations, quantum statistical perturbation theories such as thermodynamic Green functions and Feynman diagrams, path integral methods, etc. Using partial integration, we show the relation

i z A B z = A B + A ̇ B z = A B A B ̇ z . E166

Then, the generalized linear response equations (162) can be rewritten in the short form (165) with the matrix elements

P mn = B m B ̇ n + B ̇ m B ̇ n ω + i ω B m B n i ω B ̇ m B n ω + , E167
D mj = B m A ̇ j + B ̇ m A ̇ j ω + . E168

that can be interpreted as generalized transition rates (collision integral, left-hand side) and the influence of external forces (drift term, right-hand side of Eq. (165)).

Having the response parameters F n to our disposal, we can evaluate averages of the relevant observables, see Eq. (151),

B n t = B n rel t = β m F m e i ωt N mn , N mn = B m B n . E169

Eliminating F m , these average fluctuations B n t are proportional to the fields h j e i ωt .

#### 3.3.4 Force-force correlation function and static (dc) conductivity

As an example for the generalized linear response theory, we calculate the conductivity of a plasma of charged particles (electrons and ions) that is exposed to a static homogeneous electric field in x -direction: ω = 0 , E = E e x ,

H F = eE X , X = i N e x i . E170

Instead of h j , we have only one constant external field E . For the treatment of arbitrary ω to obtain the dynamical (optical) conductivity see [11, 13, 16, 17]. The conjugated variable A from Eq. (149) that couples the system to the external field is A = e X . The time derivative follows as A ̇ = e / m P , with P = i N e p x , i denoting the total momentum in x direction.

For simplicity, the ions are considered here as fixed in space because of the large mass ratio (adiabatic approximation). Then, the transport of charge is owing to the motion of the electrons. In general, the ions can also be treated as moving charged particles that contribute to the current.

A stationary state will be established in the plasma where the electrons are accelerated by the external field, but loose energy (and momentum) due to collisions with the ions. This nonequilibrium state is characterized by an electrical current that is absent in thermal equilibrium. We can take the electric current density j e l = e / m Ω P = e / Ω X ̇ as a relevant observable that characterizes the nonequilibrium state. Instead, we take the total momentum B = P = m X ̇ . The generalized linear response equations (165) and (167) read

F P ̇ P + P ̇ P ̇ = e m E P P + P P ̇ , E171

The term P ̇ P = P P e q vanishes as can be shown with Kubo’s identity, see Eq. (158). With the Kubo identity, we also evaluate the Kubo scalar product

P P = m 0 1 d λ X ̇ i βλ P eq = i m β Tr ρ eq X P = mN β . E172

The solution for response parameter F is

F = e m E mN / β + P P ̇ P ̇ P ̇ . E173

With Eq. (169) we have

j el = e m Ω P rel = m Ω F P P = σ dc E . E174

The resistance R in the static limit follows as

R = 1 σ dc = Ω β e 2 N 2 P ̇ P ̇ 1 + P P ̇ β / mN . E175

#### 3.3.5 Ziman formula for the Lorentz plasma

To evaluate the resistance R , we have to calculate the correlation functions P ̇ P ̇ and P P ̇ . For this, we have to specify the system Hamiltonian H S , which reads for the Lorentz plasma model (136)

H S = H 0 + H int = p E p a p a p + p , q V q a p + q a p , E p = 2 p 2 2 m . E176

We consider the ions at fixed positions R i so that V r = i V ei r R i . The Fourier transform V q depends for isotropic systems only on the modulus q = q and will be specified below. A realistic plasma Hamiltonian should consider also moving ions and the electron-electron interaction so that we have a two-component plasma Hamiltonian with pure Coulomb interaction between all constituents. This has been worked out [14], but is not the subject of our present work so that we restrict ourselves mainly to the simple Lorentz model.

The force P ̇ on the electrons follows from the x component of the total momentum ( p is the wave-number vector)

P = p p x a p a p , H S P = p , q V q q x a p + q a p . E177

We calculate the force-force correlation function (only x component)

P ̇ P ̇ = 0 d t e ϵ t 0 1 d λ i [ H S P t i λβ ] i [ H S P ] eq E178

in Born approximation with respect to V q . In lowest order, the force-force correlation function is of second order so that in the time evolution exp i / H S t i λβ the contribution H int of interaction to H S , Eq. (176), can be dropped as well as in the statistical operator. The averages are performed with the noninteracting ρ 0 . The product of the two commutators is evaluated using Wick’s theorem. One obtains

P ̇ P ̇ = p , p , q , q 0 d t e ϵ t 0 1 d λ e i E p E p + q t i βλ V q V q q x q x a p + q a p a p + q a p eq = p , q V q 2 δ E p E p + q f p 1 f p π q x 2 . E179

Because the x direction can be arbitrarily chosen in an isotropic system, we replace q x 2 = q x 2 + q y 2 + q z 2 / 3 = q 2 / 3 if the remaining contributions to the integrand are not depending on the direction in space.

Evaluating Eq. (175) in Born approximation, the correlation function P P ̇ β / mN can be neglected in relation to 1 because it contains the interaction strength. For the resistance, this term contributes only to higher orders of the interaction.

The force-force correlation function (179) is further evaluated using the relations 1 β d f E p d E p = f p 1 f p and δ E p E p + q = m 2 qp δ cos θ q 2 p . The q integration has to be performed in the limits 0 q 2 p . Finally the resistance can be calculated by inserting the previous expressions Eqs. (172) and (179) into Eq. (175) so that the Ziman-Faber formula is obtained,

R = m 2 Ω 3 12 π 3 3 e 2 N 2 0 d E p d f E d E 0 2 p d q q 3 V q 2 . E180

The expression for the resistance depends on the special form of the potential V q . For a pure Coulomb potential e 2 / Ω ϵ 0 q 2 , the integral diverges logarithmically as typical for Coulomb integrals. The divergency at very small values of q is removed if screening due to the plasma is taken into account. Within a many-particle approach, in static approximation the Coulomb potential is replaced by the Debye potential (146). The evaluation yields

σ dc = 3 4 2 π k B 3 / 2 4 π ϵ 0 2 m 1 / 2 e 2 1 Λ p therm E181

where the Coulomb logarithm is approximated by the value of the average p , with 2 p therm 2 / 2 m = 3 k B T / 2 . In the low-density limit, the asymptotic behavior of the Coulomb logarithm Λ is given by 1 / 2 ln n . However, this result for σ d c is not correct and can only be considered as an approximation, as discussed below considering the virial expansion of the resistivity.

#### 3.3.6 Different sets of relevant observables

After fully linearizing the statistical operator (161) with (155), we have for the electrical current density

j el = e m Ω P = m Ω n P B n P B ̇ n i ϵ F n + P P i ϵ e m E = σ dc E . E182

After deriving the Ziman formula from the force-force correlation function in the previous section, we investigate the question to select an appropriate set of relevant observables B n .

#### 3.3.7 Kubo formula

The most simple choice of relevant observables is the empty set. There are no response parameters to be eliminated. According Eq. (182), the Kubo formula

σ dc Kubo = e 2 β m 2 Ω P P irred E183

follows [18, 19]. The index “irred” denotes the irreducible part of the correlation function, because the conductivity is not describing the relation between the current and the external field, but the internal field. We will not discuss this in the present work. A similar expression can also be given for the dynamical, wave-number vector-dependent conductivity σ q ω which is related to other quantities such as the response function, the dielectric function, or the polarization function, see [5, 11, 16, 17]. Equation (183) is a fluctuation-dissipation theorem; equilibrium fluctuations of the current density are related to a dissipative property, the electrical conductivity.

The idea to relate the conductivity with the current-current auto-correlation function in thermal equilibrium looks very appealing because the statistical operator is known. The numerical evaluation by simulations can be performed for any densities and degeneracy. However, the Kubo formula (183) is not appropriate for perturbation theory. In the lowest order of interaction, we have the result σ dc Kubo , 0 = ne 2 / m ϵ (conservation of total momentum) which diverges in the limit ϵ 0 .

#### 3.3.8 Force-force correlation function

The electrical current can be considered as a relevant variable to characterize the nonequilibrium state, when a charged particle system is affected by an electrical field. We can select the total momentum as the relevant observable, B n P . Now, the character of Eq. (182) is changed. According the response equation (162), we have

P P ̇ i ϵ F + P P i ϵ e m E = 0 E184

so that these contributions compensate each other. As a relevant variable, the averaged current density is determined by the response parameter F which follows from the solution of the response equation (184). We obtain the inverse conductivity, the resistance, as a force-force auto-correlation, see Eq. (175). Now, perturbation theory can be applied, and in Born approximation a standard result of transport theory is obtained, the Ziman formula (180). We conclude that the use of relevant observables gives a better starting point for perturbation theory. In contrast to the Kubo formula that starts from thermal equilibrium as initial state, the correct current is already reproduced in the initial state and must not be created by the dynamical evolution.

However, despite the excellent results using the Ziman formula in solid and liquid metals where the electrons are strongly degenerated, we cannot conclude that the result (181) for the conductivity is already correct for low-density plasmas (nondegenerate limit if T remains constant) in the lowest order of perturbation theory considered here. The prefactor 3 / 4 2 π is wrong. If we go to the next order of interaction, divergent contributions arise. These divergences can be avoided performing a partial summation, that will also change the coefficients in Eq. (181) which are obtained in the lowest order of the perturbation expansion. The divergent contributions can also be avoided extending the set of relevant observables B n , see below.

#### 3.3.9 Higher moments of the single-particle distribution function

Besides the electrical current, also other deviations from thermal equilibrium can occur in the stationary nonequilibrium state, such as a thermal current. In general, for homogeneous systems, we can consider a finite set of moments of the single-particle distribution function

P n = p p x β E p n / 2 a p a p E185

as set of relevant observables B n . It can be shown that with increasing number of moments the result

σ dc = s k B 3 / 2 4 π ϵ 0 ) 2 m 1 / 2 e 2 1 Λ p therm E186

is improved, as can be shown with the Kohler variational principle, see [11, 15]. The value s = 3 / 4 2 π obtained from the single moment approach is increasing to the limiting value s = 2 5 / 2 / π 3 / 2 . For details see [5, 15, 16], where also other thermoelectric effects in plasmas are considered.

#### 3.3.10 Single-particle distribution function and the general form of the linearized Boltzmann equation

Kinetic equations are obtained if the occupation numbers n ν of single-(quasi-) particle states ν is taken as the set of relevant observables B n . The single-particle state ν is described by a complete set of quantum numbers, e.g., the momentum, the spin and the species in the case of a homogeneous multi-component plasma. In thermal equilibrium, the averaged occupation numbers of the quasiparticle states are given by the Fermi or Bose distribution function, n ν eq = f ν 0 = Tr ρ eq n ν . These equilibrium occupation numbers are changed under the influence of the external field. We consider the deviation Δ n ν = n ν f ν 0 as relevant observables. They describe the fluctuations of the occupation numbers. The response equations, which eliminate the corresponding response parameters F ν , have the structure of a linear system of coupled Boltzmann equations for the quasiparticles, see [11]

e m E P n ν + P n ̇ ν = ν F ν P ν ν , E187

with P ν ν = n ̇ ν Δ n ν + n ̇ ν n ̇ ν . The response parameters F ν are related to the averaged occupation numbers as

f ν t = Tr ρ t n ν = f ν 0 + β ν F ν Δ n ν Δ n ν . E188

The general form of the linear Boltzmann equation (187) can be compared with the expression obtained from kinetic theory. The left-hand side can be interpreted as the drift term, where self-energy effects are included in the correlation function P n ̇ ν . Because the operators n ν are commuting, from the Kubo identity follows n ̇ ν n ν = 1 / β n ν n ν = 0 . In the general form, the collision operator is expressed in terms of equilibrium correlation functions of fluctuations that can be evaluated by different many-body techniques. In particular, for the Lorentz model the result (186) with s = 2 5 / 2 / π 3 / 2 is obtained [5, 15, 16].

#### 3.3.11 Two-particle distribution function, bound states

Even more information is included if we also consider the nonequilibrium two-particle distributions. As an example, we mention the Debye-Onsager relaxation effect, see [5, 14]. Another important case is the formation of bound states. It seems naturally to consider the bound states as new species and to include the occupation numbers (more precisely, the density matrix) of the bound particle states in the set of relevant observables [20, 21]. It needs a long memory time to produce bound states from free states dynamically in a low-density system, because bound states cannot be formed in binary collisions, a third particle is needed to fulfill the conservation laws.

The inclusion of initial correlation to improve the kinetic theory, in particular to fulfill the conservation of total energy, is an important step worked out during the last decades, see [22] where further references are given. Other approaches to include correlations in the kinetic theory are given, e.g., in [23, 24].

#### 3.3.12 Conclusions

Transport coefficients are expressed in terms of correlation functions in equilibrium. The evaluation can be performed numerically (molecular-dynamic simulations), or using quantum statistical methods such as perturbation theory and the technique of Green functions. The generalized linear response theory has solved problems owing to the evaluation of correlation functions. Perturbation expansions are improved if higher orders are considered. The treatment of singular terms that appear in perturbation expansions is quite complex. Alternatively, the set of relevant observables can be extended. Examples are the virial expansion of the conductivity [14] or the hopping conductivity [5, 12].

It is not clear whether the rigorous evaluation of the correlation functions (i.e., the limit ϵ 0 only after full summation of the perturbation expansion) will give nontrivial results for the conductivity. For instance, arguments can be given that the exact evaluation of the force-force correlation function to calculate the resistance leads to a vanishing result, and the correlation function of stochastic forces must be considered, in analogy to the corresponding term in the Langevin equation [6, 25]. A related projection operator technique was used by Mori [26] for the memory-function approach.

There are close relation to other approaches, such as kinetic theory or quantum master equations, where the response function of the bath is considered. Irreversibility is not inherent in the equilibrium correlation functions, but in the assumption that a nonequilibrium state is considered as a fluctuation in equilibrium with a prescribed value of the relevant quantity. Other degrees of freedom are forced to adopt the distribution of thermal equilibrium.

## 4. Concluding remarks

### 4.1 Information theory

The method of nonequilibrium statistical operator (NSO) to describe irreversible processes is based on a very general concept of entropy, the Shannon information entropy (10). This concept is not restricted to dynamical properties like energy, particle numbers, momentum, etc., occurring in physics, but may be applied also to other properties occurring, e.g., in economics, financial market, and game theory. The generalized Gibbs distributions (13) and (19) are obtained if the averages of a given set of observables are known. Other statistical ensembles may be constructed, where the values of some observables have a given distribution. For instance, the canonical ensemble follows if the particle numbers are fixed, and the microcanonical ensemble has in addition a fixed energy in the interval Δ E around E , see [1, 2]. There exist alternative concepts of entropy to valuate a probability distribution which are not discussed here.

In physics, we have a dynamical evolution that forms the equilibrium distribution for ergodic systems, and any initial distribution that is compatible with the values of the conserved quantities can be used to produce the correct equilibrium distribution. The main problem is the microscopic approach to evaluate the dynamical averages, which can be done using quantum statistical methods such as Green function theory or path integral calculations, or, alternatively, numerical simulations of the microscopic equations of motion such as molecular dynamics. In more general, complex systems, we do not know the exact dynamics of the time evolution. However, we can observe time-dependent correlation functions which reflect the time evolution, and properties such as the fluctuation-dissipation theorem are not related to a specific dynamical model for the complex system. The most interesting issue of the NSO method is the selection of the set of relevant observables to describe a nonequilibrium process. The better the choice of the set of relevant observables is, for which a dynamical model for the time evolution can be found, the less influence is produced by the irrelevant observables which may be described by time-dependent correlation functions.

### 4.2 Hydrodynamics

An important application is the description of hydrodynamic processes and its relation to kinetic theory. The NSO method allows to treat this problem, selecting the single-particle distribution as well as the hydrodynamic variables as set of relevant observables. This approach has been worked out in [23]. A more general presentation is found in [4], and transport processes in multi-component fluids and superfluid systems are investigated. Until now, a rigorous theory of turbulence is not available, but hydrodynamic fluctuations and turbulent flow have been considered using the NSO method [4].

### 4.3 The limit ϵ → 0

It is the source term of the extended von Neumann equation (27) that introduces irreversible behavior. Different choices for the set of relevant observables are elected for different applications, in particular quantum master equations, kinetic theory, and linear response theory. It is claimed that this choice of the set of relevant observables is only a technical issue and has no influence on the result, only if the limit ε 0 is correctly performed in the final result.

However, calculations are not performed this way. For instance, the limit ε 0 is performed already in a finite order of perturbation theory. The self-consistency conditions (18) guarantee that a finite source term will not influence the Hamiltonian dynamics of the relevant observables. A closer investigation of a finite source term and its influence on the nonequilibrium evolution would be of interest.

### 4.4 Heat production and entropy

A serious problem is that irreversibility is connected with the production of entropy [6]. For instance, in the case of electrical conductivity, heat is produced. In principle, we have to consider an open system coupled to a bath that absorbs the produced heat. In the Zubarev NSO method considered here, it is the right-hand of the extended von Neumann equation (27) that contains the source term. We impose the stationary conditions so that ρ rel , in particular T , are not explicitly depending on time. Then, the source term acts like an additional process describing the coupling to a bath without specifying the microscopic process. The parameter ϵ now has the meaning of a relaxation time, and is no longer arbitrarily small but is of the order E 2 .

From a systematic microscopic point of view, one can introduce a process into the system Hamiltonian which describes the cooling of the system via the coupling to a bath, as known from the quantum master equations for open systems. Phonons related to the motion of ions can be absorbed by the bath, but one can calculate the electrical conductivity also for (infinitely) heavy ions so that the scattering of the electrons, accelerated by the field, is elastic. Collisions of electrons with the bath may help, but an interesting process to reduce the energy is radiation. Electrons which are accelerated during the collisions emit bremsstrahlung. This heat transfers the gain of energy of electrons, which are moving in the external field, to the surroundings.

### 4.5 Open systems: coupling to the radiation field

A general approach to scattering theory was given by Gell-Mann and Goldberger [27] (see also [1, 2]) to incorporate the boundary condition into the Schrödinger equation. The equation of motion in the potential V r reads

t ψ ϵ r t + i H ψ ϵ r t = ϵ ψ ϵ r t ψ rel t ̂ r t . E189

With H = H 0 + V , the relevant state is an eigenstate p of H 0 which changes its value at the scattering time t ̂ where the asymptotic state p is formed. As known from the Langevin equation, one can consider ψ ϵ r t = ϱ 1 / 2 exp i S / as a stochastic process [5] related to a stochastic potential V r t ; Eq. (189) appears as an average. The relaxation term is related to the fluctuations of V r t . The average Hamiltonian dynamics is realized by the self-consistency conditions for ψ rel t ̂ r t , see Eq. (28).

An interesting example is the electrical conductivity. In the stationary case which is homogeneous in time, the system remains near thermodynamic equilibrium as long as the electrical field is weak so that the produced heat can be exported. We have to consider an open system. If the conductor is embedded in vacuum, heat export is given by radiation. Bremsstrahlung is emitted during the collision of charged particles. Emission of photons can be considered as a measuring process to localize the charged particle during the collision process. The time evolution of the system is considered as a stochastic process, see also [6].

## References

1. 1. Zubarev DN. Nonequilibrium Statistical Thermodynamics. New York: Plenum Press; 1974
2. 2. Zubarev DN. The statistical operator for nonequilibrium systems. Doklady Akademii Nauk SSSR. 1961;140:92
3. 3. Zubarev D, Morozov V, Röpke G. Statistical Mechanics of Nonequilibrium Processes. Vol. 1. Berlin: Akademie-Verlag; 1996
4. 4. Zubarev D, Morozov V, Röpke G. Statistical Mechanics of Nonequilibrium Processes. Vol. 2. Berlin: Akademie-Verlag; 1997
5. 5. Röpke G. Nonequilibrium Statistical Physics. Weinheim: Wiley-VCH; 2013
6. 6. Röpke G. Electrical conductivity of charged particle systems and zubarevs nonequilibrium statistical operator method. Theoretical and Mathematical Physics. 2018;194:74
7. 7. Boltzmann L. Vorlesungen über Gastheorie. Vol. 2. Leipzig: Verlag; 1898
8. 8. Bogoliubov NN. Problems of Dynamic Theory in Statistical Physics (in Russian). Moscow-Leningrad: Gostekhizdat; 1946
9. 9. Gocke C, Röpke G. Master equation of the reduced statistical operator of an atom in a plasma. Theoretical and Mathematical Physics. 2008;154:26
10. 10. Lin C, Gocke C, Röpke G, Reinholz H. Transition rates for a Rydberg atom surrounded by a plasma. Physical Review A. 2016;93:042711
11. 11. Reinholz H, Röpke G. Dielectric function beyond the random-phase approximation: Kinetic theory versus linear response theory. Physical Review E. 2012;85:036401
12. 12. Christoph V, Röpke G. Theory of inverse linear response coefficients. Physica Status Solidi B. 1985;131:11
13. 13. Reinholz H, Redmer R, Röpke G, Wierling A. Long-wavelength limit of the dynamical local-field factor and dynamical conductivity of a two-component plasma. Physical Review E. 2000;62:5648
14. 14. Röpke G. Quantum-statistical approach to the electrical conductivity of dense, high-temperature plasmas. Physical Review A. 1988;38:3001
15. 15. Redmer R. Physical properties of dense, low-temperature plasmas. Physics Reports. 1997;282:35
16. 16. Reinholz H. Dielectric and optical properties of dense plasmas. Annales de Physique (Paris). 2005;30:1
17. 17. Röpke G. Dielectric function and electrical dc conductivity of nonideal plasmas. Physical Review E. 1998;57:4673
18. 18. Kubo R. Statistical-Mechanical theory of irreversible processes. Journal of the Physical Society of Japan. 1957;12:570
19. 19. Kubo R. The fluctuation-dissipation theorem. Reports on Progress in Physics. 1966;29:255
20. 20. Röpke G. Electrical conductivity of a system of localized and delocalized electrons. Theoretical and Mathematical Physics. 1981;46:184
21. 21. Adams J et al. Coulomb contribution to the direct current electrical conductivity of dense partially ionized plasmas. Physics of Plasmas. 2007;14:062303
22. 22. Morozov VD et al. The “Mixed” green’s function approach to quantum kinetics with initial correlations. Annals of Physics (New York). 1999;278:127
23. 23. Zubarev DN, Morozov VG, Omelyan IP, Tokarchuk MV. Unification of the kinetic and hydrodynamic approaches in the theory of dense gases and liquids. Theoretical and Mathematical Physics. 1993;96:997
24. 24. Morozov VG, Röpke G. Non-Markovian quantum kinetics and conservation laws. Journal of Statistical Physics. 2001;102:285
25. 25. Kalashnikov VP. Linear relaxation equations in the nonequilibrium statistical operator method. Theoretical and Mathematical Physics. 1978;34:263
26. 26. Mori H. A continued-fraction representation of the time-correlation functions. Progress in Theoretical Physics. 1965;34:399
27. 27. Gell-Mann M, Goldberger ML. The formal theory of scattering. Physics Review. 1953;91:398

Written By

Gerd Röpke

Submitted: 27 May 2018 Reviewed: 24 January 2019 Published: 10 May 2019