The density matrix formalism describes the statistical distribution of quantum states in a system [59]. This method allows to treat an ensemble of particles statistically.
Pure States:
The quantum subsystem is described in a Hilbert space of basis functions {Φ_{1},Φ_{2},...} and a wave function Ψ represents a generic particle from the ensemble according to
 (4.1) 
In the linear combination, the coefficents c_{i} are time dependent and describe the time propagation of the quantum subsystem [60], and
 (4.2) 
The density operator describes the probability distribution in a system. It takes the form of a projection operator and is defined by
 (4.3) 
The probability of finding a system in state Φ_{i}⟩ is described by the diagonal elements ρ_{ii}. The degree of coherence is described by the polarization between states Φ_{i}⟩ and Φ_{j}⟩, which is included in the offdiagonal elements ρ_{ij}. The total population density is conserved, which implies
 (4.4) 
The expression for the expectation value of an operator can be written as
Mixed States:
A mixed state is an incoherent mixture of pure states Ψ^{(j)}⟩, (j = 1,2,...N) with statistical weights p_{j} which obey [61]

The average value of Ō in the mixed state is given by
 (4.6) 
The density matrix depends on the pure states involved in the mixed state and their statistical weights according to
 (4.7) 
The time evolution of the density matrix is obtained by taking the time derivative of the density operator
 (4.8) 
with Ĥ = Ĥ _{0 } + Ĥ^{′}, and Ĥ^{′} represents a perturbation. This equation of motion is known as the quantum Liouville equation. Due to the linear transformation of the density operator on the r.h.s in equation (4.8), it is possible to define a linear operator which yields
 (4.9) 
where denotes the so called Liouville operator

If N is the number of states in the system, the superoperator is a N^{2} × N^{2} matrix and ρ is a N^{2} dimensional vector. Thus, systems with many states may give rise to insuperable computational challenges.
We consider a system which is open and in permanent contact with its environment. Under certain conditions the system which is initially in a nonequilibrium state will go over into an equilibrium state after some time. This relaxation process is irreversible.
Let S denote an open quantum system which is coupled to another quantum system R, the environment. The density matrix characterizing the total system is denoted by ρ(t) and the total Hamiltonian by Ĥ = Ĥ_{S} + Ĥ_{R} + Ĥ^{′}. In the Schrödinger picture, the equation of motion for the density matrix is given by the Liouvillevon Neumann equation
 (4.10) 
where Ĥ _{0 } = Ĥ _{S } + Ĥ_{R}. Within the interaction picture, the time evolution of the density matrix is controlled only by the interaction Hamiltonian according to
 (4.11) 
Inserting the formal integration

into the r.h.s of equation (4.11), the equation of motion for ρ_{I} can be expressed as
 (4.12) 
The reduced density matrix ρ_{I}^{S} which describes the system S is obtained by taking the trace over all variables of system R according to [62]
 (4.13) 
It is assumed that the interaction starts at t = 0 and that the two systems are uncorrelated prior to this time. Thus
 (4.14) 
In order to have an irreversible process, it is assumed that R has so many degrees of freedom that it is described by a thermal equlibrium distribution, independent on the amount of energy transferred to it from the system S [63]. Hence
 (4.15) 
This is the fundamental condition for irreversibility. Using these approximations, the equation of motion for the reduced density matrix can be written as
 (4.16) 
According to this equation, the behavior of the system depends on the past events in the time interval [0, t], since the integral contains ρ_{I}^{S}(t^{′}). However, the system S is coupled to the reservoir which causes a damping that destroys the ”knowledge” of the past behavior. These considerations lead to the assumption that the system loses all memory of the past. Thus, the substitution
 (4.17) 
is necessary, which is known as the Markov approximation.
In general, the interaction operator can be expressed as
 (4.18) 
where ^{S } and ^{R} denote the operators of the dynamic system and the reservoir, respectively. In the interaction picture we have
 (4.19) 
and
 (4.20) 
Inserting the interaction operator of this form and the substitution (4.17) into the equation (4.16) and taking into account that the operators _{i}^{S} and _{i}^{R} commute, yields
Let r⟩ denote the eigenstates of Ĥ_{R}. Then, the expectation values appearing in the first term in equation (4.21) can be written as

describe the average correlation between interactions occuring at times t and t^{′}. The reservoir dissipates quickly the effects of its interaction with the system S. Thus

where τ_{R } denotes the correlation time of the reservoir. For t  t^{′} > τ_{R}, interactions become less correlated and the correlation functions satisfy

On the basis of the last considerations we will discuss the Markov approximation. It follows that the integral in equation (4.16) is nonzero only for t^{′} in the time interval [t  τ_{R},t]. Outside this interval, the values of ρ_{I}^{S}(t^{′}) have negligible influence on the time derivative of ρ_{I}^{S}(t). The Markov approximation holds if τ_{R} is much smaller than a characteristic time that is required for ρ_{I}^{S}(t) to change considerably on a macroscopic scale [65].
Since the correlation functions are stationary and depend only on the time difference t^{′′} = t  t^{′}

the equation of motion of the reduced density matrix can be written as

one obtains after some algebra
 (4.23) 
where the following abbreviations are used
 (4.24) 
where
 (4.25) 
and
 (4.26) 
Equation (4.24) is called the generalized Master equation and describes the irreversible behavior of a system. In order to get an interpretation of the occuring parameters, we take a look at the rate of change of the diagonal elements of the density matrix. For the diagonal elements, the Schrödinger picture is equivalent to the interaction picture. From the results determined above, it is straightforward to obtain
 (4.27) 
This equation is called the Pauli Master equation which plays a fundamental role in modern statistics and can be used in many problems in physics and chemical kinetics [67]. The diagonal elements of the density matrix give the probability of finding the corresponding level occupied at time t. Since the probability increases due to transitions from γ⟩ to β⟩, and decreases as a result of transitions from β⟩ to other states γ⟩, the parameters S_{βγ} describe the probability per unit time for the transitions γ⟩ → β⟩.
Now we take a more detailed look at the transitions rates. Making use of equation (4.20), the transition rate can be written as
This is the ”Golden Rule” for a transition from level β⟩ to level γ⟩. In order to ensure the energy conservation

the reservoir undergoes simultaneously a transition from a state r^{′}⟩ to r⟩. The diagonal elements ⟨r^{′}ρ^{R}(0)r^{′}⟩ describe the probability of finding the reservoir in a state with energy E_{r}. Reservoirs are in thermal equilibrium. Hence, such a system is represented by an incoherent sum of energy eigenstates with statistical weights^{1 } . Equation (4.27) holds for closed systems. For open systems, a term describing the interaction with the electrical borders has to be added.
 (4.29) 
The derivation of the Pauli Master equation can be simplified, if the interaction term in equation 4.10 is approximated by a relaxation term
 (4.30) 
The last term on the r.h.s expresses the weak interaction with the surroundings. The limit τ_{s} →∞ will be taken in the end. Considering the matrix element of the above equation over the eigenstates α⟩ of the unperturbed Hamiltonian Ĥ_{0}, one gets
 (4.31) 
(∂ρ∕∂t)_{res } is assumed to be of order O(α_{c}^{0}), where α_{c} denotes the strength of the interaction Hamiltonian. Equation (4.31) implies that the off diagonal elements of the density matrix are of order O(α_{c }^{1 }). To the lowest order, O(α_{c}^{1}), one can obtain

which implies that ρ_{ββ} is of the order O(α_{c}^{2}). Rewriting the last equation as

and inserting into eqaution (4.31), yields
 (4.32) 
Making use of

one gets for τ_{s } →∞
 (4.33) 
which is formally identical to equation 4.29.