Density Matrix Formulation

4.1 Density Matrix Formulation

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 {Φ₁,Φ₂,...} and a wave function Ψ represents a generic particle from the ensemble according to

∑ |Ψ⟩ = ci|Φi⟩ i

(4.1)

In the linear combination, the coefficents c_i are time dependent and describe the time propagation of the quantum subsystem [60], and

ci = ⟨Φi |Ψ ⟩

(4.2)

The density operator describes the probability distribution in a system. It takes the form of a projection operator and is defined by

ˆρ(t) = |Ψ(t)⟩⟨Ψ (t)|

(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 off-diagonal elements ρ_ij. The total population density is conserved, which implies

∑ Tr(ˆρ) = ρii = 1 i

(4.4)

The expression for the expectation value of an operator can be written as

∑ ∑ ⋆ ∑ ∑ ⋆ ⟨Ψ |Oˆ|Ψ ⟩ = cicj⟨Φi |Oˆ|Φj ⟩ = cicj ˆOij i j i j = ∑ ∑ ρ⋆ ˆO = ∑ ∑ ρ Oˆ ji ij ij ij i j i j = Tr [ρOˆ] (4.5)

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]

N ∑ p = 1 j j=1

The average value of Ô in the mixed state is given by

∑N ⟨ˆO⟩ = pj⟨Ψ (j)|Oˆ|Ψ (j)⟩ j=1

(4.6)

The density matrix depends on the pure states involved in the mixed state and their statistical weights according to

N ∑ (j) (j) ˆρ = pj|Ψ ⟩⟨Ψ | j=1

(4.7)

The time evolution of the density matrix is obtained by taking the time derivative of the density operator

∂ρˆ ˆ iℏ∂t = [ˆρ,H ]

(4.8)

with Ĥ = Ĥ ₀ + Ĥ^′, 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 ˆ
L which yields

∂ˆρ i ∂t-= - ℏ-ˆLˆρ

(4.9)

where ˆ
L denotes the so called Liouville operator

Lˆ = ˆH δ - ˆH ⋆δ ij,mn im jn jn im

If N is the number of states in the system, the superoperator is a N² × N² matrix and ρ is a N² dimensional vector. Thus, systems with many states may give rise to insuperable computational challenges.

4.1.1 Density Matrix Equations for Dissipative Systems

We consider a system which is open and in permanent contact with its environment. Under certain conditions the system which is initially in a non-equilibrium 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 Liouville-von Neumann equation

∂ρ(t) i i ′ -∂t--= -ℏ-[H ˆ0, ρ(t)]- ℏ[ ˆH ,ρ(t)]

(4.10)

where Ĥ ₀ = Ĥ _S + Ĥ_R. Within the interaction picture, the time evolution of the density matrix is controlled only by the interaction Hamiltonian according to

∂ρI(t) i ′ --∂t--= -ℏ-[H ˆ (t),ρI(t)]

(4.11)

Inserting the formal integration

∫t i- ′ ˆ′ ′ ′ ρI(t) = ρI(0)- ℏ dt [H (t),ρI(t )] 0

into the r.h.s of equation (4.11), the equation of motion for ρ_I can be expressed as

∫t ∂-ρI(t) = - i[ ˆH (t),ρ (0)]+-1- dt′[Hˆ′(t),[H ˆ′(t′),ρ (t′)]] ∂t ℏ I ℏ2 I 0

(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]

S ρI(t) = TrR (ρI(t))

(4.13)

It is assumed that the interaction starts at t = 0 and that the two systems are uncorrelated prior to this time. Thus

ρ (0) = ρS(0)ρR(0) I

(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

ρ (t) = ρS(t)ρR(0) I I

(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

∂ρS(t) i 1 ∫t I = - -TrR [H ˆ′(t),ρS(0)ρR(0)]---2 dt′TrR[ ˆH ′(t),[H ˆ′(t′),ρSI(t′)ρR(0)]] ∂t ℏ ℏ 0

(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

S ′ S ρI(t) → ρI (t)

(4.17)

is necessary, which is known as the Markov approximation.

4.1.1.1 Markov Approximation

In general, the interaction operator can be expressed as

∑ Hˆ′ = ˆQSi ˆQRi i

(4.18)

where ^S and ^R denote the operators of the dynamic system and the reservoir, respectively. In the interaction picture we have

QˆS (t) = eiˆHSt∕ℏQˆS e-iˆHSt∕ℏ i i

(4.19)

and

ˆQRi (t) = eiˆHRt∕ℏQˆRi e- iHˆRt∕ℏ

(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

[ ] ∂ρSI(t) i-∑ ˆS S ˆR R S ˆS ˆR R ∂t=- ℏ Q i (t)ρI(0)TrR (Q i (t)ρ (0))- ρI(0)Q i (t)TrR (Qi (t)ρ (0)) i 1 ∑ ∫t [ - -2- dt′ (ˆQSi (t)QˆSj (t′)ρSI(t)- ˆQSj (t′)ρSI(t)QˆSi (t))TrR (ˆQRi (t)ˆQRj (t′)ρR(0)) ℏ ij 0 ] - (ˆQSi (t)ρSI(t)ˆQSj (t′)- ρSI(t)ˆQSj (t′)QˆSi (t))TrR (ˆQRj (t′)QˆRi (t)ρR(0)) (4.21)

Let |r⟩ denote the eigenstates of Ĥ_R. Then, the expectation values appearing in the first term in equation (4.21) can be written as

⟨ˆQR (t)⟩ = Tr (ˆQR (t)ρR (0)) i ∑ R i = ⟨r|ˆQRi (t)ρR (0 )|r⟩ r

Under the assumption that the interaction operators don’t have diagonal elements in this representation (no average energy shifts), the first term in equation (4.21) vanishes, since ⟨

_i^R(t)⟩ = 0 [64]. The time correlation functions

⟨ˆQR (t)QˆR (t′)⟩ = Tr (ˆQR (t)QˆR (t′)ρR(0)) i j R i j

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

R R ′ ′ ⟨Qˆi (t)ˆQj (t )⟩ ⁄= 0, fort- t ≲ τR

where τ_R denotes the correlation time of the reservoir. For t - t^′ > τ_R, interactions become less correlated and the correlation functions satisfy

ˆ R ˆR ′ ⟨Q i (t)Q j (t)⟩ ≈ 0

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].

4.1.1.2 The Relaxation Equation

Since the correlation functions are stationary and depend only on the time difference t^′′ = t - t^′

⟨ˆQR(t)ˆQR (t′)⟩ = ⟨QˆR (t- t′)QˆR ⟩ i j i j

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

S ∑ ∫∞ ( ∂ρI-(t) = - 1-- dt′′ [QˆSi (t),Q ˆSj (t- t′′)ρSI(t)]⟨ˆQRi (t′′)QˆRj ⟩ ∂t ℏ2 ij 0 ) - [QˆS (t),ρS (t)QˆS (t- t′′)]⟨ˆQR ˆQR (t′′)⟩ (4.22) i I j j i

where the upper integration limit is extended to infinity with negligible error under the Markov approximation. Let {|β⟩,|β^′⟩,...} denote the eigenstates of Ĥ_S, where |β⟩ = |λ,η,ν,k_∥⟩. Here, η is the subband index, ν denotes the valley index, λ is the index of the stage, and k_∥ is the in-plane wave vector. Then

′ˆ S iωβ′β ′ˆS ⟨β |Q i (t)|β⟩ = e ⟨β |Qi |β⟩, ω β′β = (Eβ′ - Eβ)∕ℏ

one obtains after some algebra

S ∑ ′ ⟨β′|∂-ρI(t) |β⟩ = ⟨γ′|ρSI(t)|γ⟩R γβγ′βei(Eβ′- Eβ-Eγ′+Eγ)t∕ℏ ∂t γ′γ

(4.23)

where the following abbreviations are used

( ) β′β ∑ ∑ Rγ′γ= - δγβΓ +γ′μμβ′ + Γ +βγγ′β′ + Γ -βγγ′β′ - δγ′β′Γ -βμμγ ei(ωγ′β′+ ωγβ)t μ μ ∑ ∞∫ Γ+= 1-- ⟨ζ|ˆQS |μ ⟩⟨γ|ˆQS |β ⟩ dt′′e-iωγβt′′⟨ˆQR (t′′)QˆR ⟩ ζμγβ ℏ2 ij i j i j 0∞ 1 ∑ ∫ ′′ Γ-ζμγβ= -2- ⟨ζ|ˆQSj |μ ⟩⟨γ|ˆQSi |β ⟩ dt′′e-iωζμt ⟨ˆQRj ˆQRi (t′′)⟩ ℏ ij 0

Considering the case of energy conservation, the time dependent exponential vanishes. It can be shown that after transformation into the Schrödinger picture the equations of motion of the reduced density matrices can be written as [66]

′∂ρS(t) i- ′ ˆ S ∑ S ′ S ⟨β|∂t |β⟩ = - ℏ⟨β |[HS, ρ (t)]|β ⟩+ δββ′ ⟨γ|ρ (t)|γ⟩Sβγ - χβ′β⟨β |ρ (t)|β⟩ γ

(4.24)

where

+ - Sγβ = Γβγγβ + Γβγγβ

(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

d ∑ ∑ --ρββ(t) = ργγ(t)Sβγ - ρββ(t) Sγβ dt γ⁄=β γ⁄=β

(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

+ - Sγβ=Γβγγβ + Γβγγβ ∑ ∑ ∫∞ =1- ⟨r|ˆQSi |γ ⟩⟨γ|ˆQSj |r⟩⟨r′|QˆRi |r⟩⟨r|ˆQRj |r′⟩⟨r′|ρR|r′⟩ dt′′ei(Er′- Er- ℏωγβ)t′′∕ℏ ℏ2 ij rr′ 0∞ 1 ∑ S S ′ R R ′ ′ R ′∫ ′′ i(E -E - ℏω )t′′∕ℏ +ℏ2- ⟨r|Qˆj |γ⟩⟨γ |Qˆi |r⟩⟨r |Qˆj |r⟩⟨r|Qˆi |r⟩⟨r|ρ (0)|r⟩ dt e r′ r βγ rr′ 0 2π ∑ ′ ′ 2 ′ R ′ =ℏ- |⟨γ,r|Hˆ|β,r ⟩|⟨r|ρ (0)|r⟩δ(Er′ - Er - ℏω γβ) (4.28) rr′

This is the ”Golden Rule” for a transition from level |β⟩ to level |γ⟩. In order to ensure the energy conservation

Er - Er′ = Eβ - E γ

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¹ . Equation (4.27) holds for closed systems. For open systems, a term describing the interaction with the electrical borders has to be added.

d ∑ ∑ (∂ ρββ) --ρββ(t) = ργγ(t)S βγ - ρββ(t) Sγβ + ----- dt γ⁄= β γ⁄=β ∂t res

(4.29)

4.1.1.3 From the Liouville-von Neumann Equation to the PME

The derivation of the Pauli Master equation can be simplified, if the interaction term in equation 4.10 is approximated by a relaxation term

dρ i ( ∂ρ) ρ- ρeq ---= -[ρ,Hˆ]+ --- - ------- dt ℏ ∂t res τs

(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 Ĥ₀, one gets

dρββ′= iρ ′(E - E ′ + iℏ ∕τ) + i ∑ [ρ ′′⟨β ′′|Hˆ |β ⟩- ⟨β′|H ˆ |β′′⟩ρ ′′ ′] dt ℏ ββ β β s ℏ ′′ ′ ββ int int β β β ⁄=β,β ( ) + i⟨β|Hˆ |β ′⟩(ρ - ρ ′′)+ -iρ ′[⟨β| ˆH |β⟩- ⟨β ′| ˆH |β′⟩]+ ∂ρ- ℏ int ββ ββ ℏ ββ int int ∂t res

Introducing renormalized energies E_β → E_β + ⟨β|Ĥ_int|β⟩, the time evolution of the diagonal elements of the density matrix can be written as

d ρββ i ∑ ′ ′ ( ∂ρ ) --dt- = ℏ- [ρββ′⟨β | ˆHint|β ⟩- ⟨β| ˆHint|β ⟩ρβ′β]+ ∂t- β′⁄=β res

(4.31)

(∂ρ∕∂t)_res is assumed to be of order O(α_c⁰), 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