A short introduction to the Lindblad master equation

Manzano, Daniel

doi:10.1063/1.5115323

The theory of open quantum systems is one of the most essential tools for the development of quantum technologies. Furthermore, the Lindblad (or Gorini-Kossakowski-Sudarshan-Lindblad) master equation plays a key role as it is the most general generator of Markovian dynamics in quantum systems. In this paper, we present this equation together with its derivation and methods of resolution. The presentation tries to be as self-contained and straightforward as possible to be useful to readers with no previous knowledge of this field.

I. INTRODUCTION

Open quantum system techniques are vital for many studies in quantum mechanics.^1–3 This happens because closed quantum systems are just an idealization of real systems (the same happens with closed classical systems), as in Nature nothing can be isolated. In practical problems, the interaction of the system of interest with the environment cannot be avoided and we require an approach in which the environment can be effectively removed from the equations of motion.

The general problem addressed by open quantum theory is sketched in Fig. 1. In the most general picture, we have a total system that conforms a closed quantum system by itself. We are mostly interested in a subsystem of the total one (we call it just “system” instead of “total system”). Therefore, the whole system is divided into our system of interest and an environment. The goal of open quantum theory is to infer the equations of motions of the reduced system from the equation of motion of the total system. For practical purposes, the reduced equations of motion should be easier to solve than the full dynamics of the system. Because of his requirement, several approximations are usually made in the derivation of the reduced dynamics.

FIG. 1.

View large Download slide

A total system divided into the system of interest, “system” and the environment.

One particular, and interesting, case of study is the dynamics of a system connected to several baths modeled by a Markovian interaction. In this case, the most general quantum dynamics is generated by the Lindblad equation (also called Gorini-Kossakowski-Sudarshan-Lindblad equation).^4,5 It is difficult to overemphasize the importance of this master equation. It plays an important role in fields as quantum optics,^1,6 condensed matter,^7–10 atomic physics,^11,12 quantum information,^13,14 decoherence,^15,16 and quantum biology.^17–19

The purpose of this paper is to provide basic knowledge about the Lindblad master equation. In Sec. II, the mathematical requirements are introduced, while in Sec. III, there is a brief review of quantum mechanical concepts that are required to understand the paper. Section IV includes a description of a mathematical framework, the Fock-Liouville space (FLS) that is especially useful to work in this problem. In Sec. V, we define the concept of completely positive and trace-preserving maps (CPT-maps), we derive the Lindblad master equation from two different approaches, and we discus several properties of the equation. Finally, Sec. VI is devoted to the resolution of the master equation using different methods. To deepen the techniques of solving the Lindblad equation, an example consisting of a two-level system with decay is analyzed, illustrating the content of every section. The problems proposed are solved by the use of Mathematica notebooks that can be found in Ref. 20.

II. MATHEMATICAL BASIS

The primary mathematical tool in quantum mechanics is the theory of Hilbert spaces. This mathematical framework allows extending many results from finite linear vector spaces to infinite ones. In any case, this tutorial deals only with finite systems and, therefore, the expressions “Hilbert space” and “linear space” are equivalent. We assume that the reader is skilled in operating in Hilbert spaces. To deepen the field of Hilbert spaces, we recommend the book by Debnath and Mikusińki.²¹ If the reader needs a brief review of the main concepts required for understanding this paper, we may recommend Nielsen and Chuang’s Quantum Computing book.²² It is also required some basic knowledge about infinitesimal calculus, such as integration, derivation, and the resolution of simple differential equations. To help the readers, we have made a glossary of the mostly used mathematical terms. It can also be used as a checklist of terms that the reader should be familiar with.

Glossary:

$H$ represents a Hilbert space, usually the space of pure states of a system.
$| ψ (⟩ \in H$ represents a vector of the Hilbert space $H$ (a column vector).
$⟨) ψ | \in H$ represents a vector of the dual Hilbert space of $H$ (a row vector).
$⟨ ψ | ϕ ⟩ \in C$ is the scalar product of vectors |ψ⟩ and |ϕ⟩.
|||ψ⟩|| is the norm of vector |ψ⟩. $| | | ψ (⟩ | | \equiv \sqrt{⟨ ψ | ψ ⟩}$ ⁠.
$B (H)$ represents the space of bounded operators acting on the Hilbert space $B : H \to H$ ⁠.
$1_{H} \in B (H)$ is the identity operator of the Hilbert space $H$ s.t. $1_{H} | ψ (⟩ = | ψ (⟩, \forall | ψ (⟩ \in H$ ⁠.
$| ψ (⟩ ⟨) ϕ | \in B (H)$ is the operator such that $(| ψ (⟩ ⟨) ϕ |) | φ (⟩ = ⟨ ϕ | φ ⟩ | ψ (⟩, \forall | φ (⟩ \in H$ ⁠.
$O^{†} \in B (H)$ is the Hermitian conjugate of the operator $O \in B (H)$ ⁠.
$U \in B (H)$ is a unitary operator iff $U U^{†} = U^{†} U = 1$ ⁠.
$H \in B (H)$ is a Hermitian operator iff H = H^†.
$A \in B (H)$ is a positive operator $(A > 0)$ iff $⟨) ϕ | A | ϕ (⟩ \geq 0, \forall | ϕ (⟩ \in H$ ⁠.
$P \in B (H)$ is a projector iff PP = P.
$Tr [B]$ represents the trace of operator B.
$ρ (L)$ represents the space of density matrices, meaning the space of bounded operators acting on $H$ with trace 1 and positive.
|ρ⟩⟩ is a vector in the Fock-Liouville space.
$⟨ ⟨ A | B ⟩ ⟩ = Tr [A^{†} B]$ is the scalar product of operators $A, B \in B (H)$ in the Fock-Liouville space.
$\tilde{L}$ is the matrix representation of a superoperator in the Fock-Liouville space.

III. (VERY SHORT) INTRODUCTION TO QUANTUM MECHANICS

The purpose of this chapter is to refresh the main concepts of quantum mechanics necessary to understand the Lindblad master equation. Of course, this is NOT a full quantum mechanics course. If a reader has no background in this field, just reading this chapter would be insufficient to understand the remaining of this tutorial. Therefore, if the reader is unsure of his/her capacities, we recommend to go first through a quantum mechanics course or to read an introductory book carefully. There are many great quantum mechanics books in the market. For beginners, we recommend Sakurai’s book²³ or Nielsen and Chuang’s Quantum Computing book.²² For more advanced students, looking for a solid mathematical description of quantum mechanics methods, we recommend Galindo and Pascual.²⁴ Finally, for a more philosophical discussion, you should go to Peres’ book.²⁵

We start stating the quantum mechanics postulates that we need to understand the derivation and application of the Lindblad master equation. The first postulate is related to the concept of a quantum state.

Postulate 1.

Associated with any isolated physical system, there is a complex Hilbert space $H$ , known as the state space of the system. The state of the system is entirely described by a state vector, which is a unit vector of the Hilbert space $(| ψ (⟩ \in H)$ .

As quantum mechanics is a general theory (or a set of theories), it does not tell us which is the proper Hilbert space for each system. This is usually done system by system. A natural question to ask is if there is a one-to-one correspondence between unit vectors and physical states, meaning that if every unit vector corresponds to a physical system. This is resolved by the following corollary that is a primary ingredient for quantum computation theory (see Ref. 22, Chapter 7).

Corollary 1.

All unit vectors of a finite Hilbert space correspond to possible physical states of a system.

Unit vectors are also called pure states. If we know the pure state of a system, we have all physical information about it and we can calculate the probabilistic outcomes of any potential measurement (see the next postulate). This is a very improbable situation as experimental settings are not perfect, and in most cases, we have only imperfect information about the state. Most generally, we may know that a quantum system can be in one state of a set

\{| ψ_{i} (⟩\}

with probabilities p_i. Therefore, our knowledge of the system is given by an ensemble of pure states described by the set

\{| ψ_{i} (⟩, p_{i}\}

⁠. If more than one p_i is different from zero, the state is not pure anymore and it is called a mixed state. The mathematical tool that describes our knowledge of the system, in this case, is the density operator (or density matrix),

ρ \equiv \sum_{i} p_{i} | ψ_{i} (⟩ ⟨) ψ_{i} | .

(1)

Density matrices are bounded operators that fulfill two mathematical conditions,

A density matrix ρ has unit trace $(Tr [ρ] = 1)$ ⁠.
A density matrix is a positive matrix ρ > 0.

Any operator fulfilling these two properties is considered a density operator. It can be proved trivially that density matrices are also Hermitian.

If we are given a density matrix, it is easy to verify if it belongs to a pure or a mixed state. For pure states, and only for them, Tr[ρ²] = Tr[ρ] = 1. Therefore, if Tr[ρ²] < 1, the system is mixed. The quantity Tr[ρ²] is called the purity of the states, and it fulfills the bounds $\frac{1}{d} \leq Tr [ρ^{2}] \leq 1$ ⁠, being d the dimension of the Hilbert space.

If we fix an arbitrary basis

{\{| i (⟩\}}_{i = 1}^{N}

of the Hilbert space, the density matrix in this basis is written as

ρ = \sum_{i, j = 1}^{N} ρ_{i, j} | i (⟩ ⟨) j |

or

ρ = (\begin{matrix} ρ_{00} & ρ_{01} & \dots & ρ_{0 N} \\ ρ_{10} & ρ_{11} & \dots & ρ_{1 N} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ ρ_{N 0} & ρ_{N 1} & \dots & ρ_{N N} \end{matrix}),

(2)

where the diagonal elements are called populations

() ρ_{i i} \in R_{0}^{+}

and

\sum_{i} ρ_{i, i} = 1))

⁠, while the off-diagonal elements are called coherences

(ρ_{i, j} \in C and ρ_{i, j} = ρ_{j, i}^{*})

⁠. Note that this notation is base-dependent.Example 1: State of a two-leve system (qubit).

The Hilbert space of a two-level system is just the two-dimensional linear space $H_{2}$ ⁠. Examples of this kind of system are $\frac{1}{2}$ -spins and two-level atoms. We can define a basis of it by the orthonormal vectors, $\{| 0 (⟩, | 1 (⟩\}$ ⁠. A pure state of the system would be any unit vector of $H_{2}$ ⁠. It can always be expressed as |ψ⟩ = a|0⟩ + b|1⟩ with $a, b \in C$ s.t. |a|² + |b|² = 1.

A mixed state is therefore represented by a positive unit trace operator

ρ \in O (H_{2})

⁠,

ρ = (\begin{matrix} ρ_{00} & ρ_{01} \\ ρ_{10} & ρ_{11} \end{matrix}) = ρ_{00} | 0 (⟩ ⟨) 0 | + ρ_{01} | 0 (⟩ ⟨) 1 | + ρ_{10} | 1 (⟩ ⟨) 0 | + ρ_{11} | 1 (⟩ ⟨) 1 |,

(3)

and it should fulfill ρ₀₀ + ρ₁₁ = 1 and

ρ_{01} = ρ_{10}^{*}

⁠.End of example.

Once we know the state of a system, it is natural to ask about the possible outcomes of experiments (see Ref. 23, Sec. 1.4).

Postulate 2.

All possible measurements in a quantum system are described by a Hermitian operator or observable. Due to the spectral theorem, we know that any observable O has a spectral decomposition in the form (for simplicity, we assume a nondegenerated spectrum)

O = \sum_{i} a_{i} | a_{i} (⟩ ⟨) a_{i} |,

(4)

being

a_{i} \in R

the eigenvalues of the observable and |a_i⟩ their corresponding eigenvectors. The probability of obtaining the result a_i when measuring the property described by observable O in a state |ψ⟩ is given by

P (a_{i}) = {|⟨ ψ | a_{i} ⟩|}^{2} .

(5)

After the measurement, we obtain the state |a_i⟩ if the outcome a_i was measured. This is called the postmeasurement state.

This postulate allows us to calculate the possible outputs of a system, the probability of these outcomes, as well as the after-measurement state. A measurement usually changes the state, as it can only remain unchanged if it was already in an eigenstate of the observable.

It is possible to calculate the expectation value of the outcome of a measurement defined by operator O in a state |ψ⟩ by just applying the simple formula,

⟨ O ⟩ = ⟨) ψ | O | ψ (⟩ .

(6)

With a little algebra, we can translate this postulate to mixed states. In this case, the probability of obtaining an output a_i that corresponds to an eigenvector |a_i⟩ is

P (a_{i}) = Tr [| a_{i} (⟩ ⟨) a_{i} | ρ],

(7)

and the expectation value of operator O is

⟨ O ⟩ = Tr [O ρ] .

(8)

Example 2: Measurement in a two-level system.

A possible test to perform in our minimal model is to measure the energetic state of a system, assuming that both states have a different energy. The observable corresponding to this measurement would be

H = E_{0} | 0 (⟩ ⟨) 0 | + E_{1} | 1 (⟩ ⟨) 1 | .

(9)

This operator has two eigenvalues

\{E_{0}, E_{1}\}

with two corresponding eigenvectors

\{| 0 (⟩, | 1 (⟩\}

⁠.

If we have a pure state ψ = a|0⟩ + b|1⟩, the probability of measuring the energy E₀ would be P(E₀) = |⟨0|ψ⟩|² = |a|². The probability of finding E₁ would be P(E₁) = |⟨1|ψ⟩|² = |b|². The expected value of the measurement is ⟨H⟩ = E₀|a|² + E₁|b|².

In the more general case of having a mixed state ρ = ρ₀₀|0⟩⟨0| + ρ₀₁|0⟩⟨1| + ρ₁₀|1⟩⟨0| + ρ₁₁|1⟩⟨1|, the probability of finding the ground state energy is $P (0) = Tr [| 0 (⟩ ⟨) 0 | ρ] = ρ_{00}$ and the expected value of the energy would be $⟨ H ⟩ = Tr [H ρ] = E_{0} ρ_{00} + E_{1} ρ_{11}$ ⁠.End of example.

Another natural question to ask is how quantum systems evolve. The time-evolution of a pure state of a closed quantum system is given by the Schrödinger equation (see Ref. 24, Sec. 2.9).

Postulate 3.

Time evolution of a pure state of a closed quantum system is given by the Schrödinger equation,

\frac{d}{d t} | ψ (t) (⟩ = - i ℏ H | ψ (t) (⟩,

(10)

where H is the Hamiltonian of the system and it is a Hermitian operator of the Hilbert space of the system state (from now on we avoid including Planck’s constant by selecting the units such that ℏ = 1).

The Hamiltonian of a system is the operator corresponding to its energy, and it can be nontrivial to realize.

The Schrödinger equation can be formally solved in the following way. If at t = 0, the state of a system is given by |ψ(0)⟩; at time t, it will be

| ψ (t) (⟩ = e^{- i H t} | ψ (0) (⟩ .

(11)

As H is a Hermitian operator, the operator U = e^−iHt is unitary. This gives us another way of phrasing postulate 3.

Postulate 3′.

The evolution of a closed system is given by a unitary operator of the Hilbert space of the system,

| ψ (t) (⟩ = U | ψ (0) (⟩,

(12)

with

U \in B (H)

s.t.

U U^{†} = U^{†} U = 1

⁠.

It is easy to prove that unitary operators preserve the norm of vectors and, therefore, transform pure states into pure states. As we did with the state of a system, it is reasonable to wonder if any unitary operator corresponds to the evolution of a real physical system. The answer is yes.

Lemma 1.

All unitary evolutions of a state belonging to a finite Hilbert space can be constructed in several physical realizations such as photons and cold atoms.

The proof of this lemma can be found in Ref. 22.

The time evolution of a mixed state can be calculated just by combining Eqs. and , giving the von Neumann equation

\dot{ρ} = - i [H, ρ] \equiv L ρ,

(13)

where we have used the commutator

[A, B] = A B - B A

and

L

is the so-called Liouvillian superoperator.

It is easy to prove that the Hamiltonian dynamics does not change the purity of a system,

\frac{d}{d t} Tr [ρ^{2}] = Tr [\frac{d ρ^{2}}{d t}] = Tr [2 ρ \dot{ρ}] = - 2 i Tr [ρ (H ρ - ρ H)] = 0,

(14)

where we have used the cyclic property of the trace. This result illustrates that the mixing rate of a state does not change due to the quantum evolution.Example 3: Time evolution of a two-level system.

The evolution of our isolated two-level system is described by its Hamiltonian,

H_{f r e e} = E_{0} | 0 (⟩ ⟨) 0 | + E_{1} | 1 (⟩ ⟨) 1 | .

(15)

As the states |0⟩ and |1⟩ are Hamiltonian eigenstates, if at t = 0 the atom is at the excited state |ψ(0)⟩ = |1⟩ after a time t, the state would be

| ψ (t) (⟩ = e^{- i H t} | 1 (⟩ = e^{- i E_{1} t} | 1 (⟩

⁠.

As the system was already in an eigenvector of the Hamiltonian, its time-evolution consists only in adding a phase to the state, without changing its physical properties. (If an excited state does not change, why do atoms decay?) Without losing any generality, we can fix the energy of the ground state as zero, obtaining

H_{f r e e} = E | 1 (⟩ ⟨) 1 |,

(16)

with E ≡ E₁.

To make the model more interesting, we can include a driving that coherently switches between both states. The total Hamiltonian would be then

H = E | 1 (⟩ ⟨) 1 | + Ω (| 0 (⟩ ⟨) 1 | + | 1 (⟩ ⟨) 0 |),

(17)

where Ω is the frequency of driving.

By using the von Neumann equation (13), we can calculate the populations $(ρ_{00}, ρ_{11})$ as a function of time. The system is then driven between the states, and the populations present Rabi oscillations, as it is shown in Fig. 2.End of example.

Finally, as we are interested in composite quantum systems, we need to postulate how to work with them.

FIG. 2.

View large Download slide

Population dynamics under a quantum dynamics (parameters are Ω = 1 and E = 1). The blue line represents ρ₁₁ and the orange one ρ₀₀.

Postulate 4.

The state-space of a composite physical system, composed by N subsystems, is the tensor product of the state space of each component $H = H_{1} \otimes H_{2} \otimes \dots \otimes H_{N}$ . The state of the composite physical system is given by a unit vector of $H$ . Moreover, if each subsystem belonging to $H_{i}$ is prepared in the state |ψ_i⟩, the total state is given by |ψ⟩ = |ψ₁⟩ ⊗ |ψ₂⟩ ⊗⋯⊗|ψ_N⟩.

The symbol ⊗ represents the tensor product of Hilbert spaces, vectors, and operators. If we have a composited mixed state where each component is prepared in the state ρ_i, the total state is given by ρ = ρ₁ ⊗ ρ₂ ⊗⋯⊗ ρ_N.

States that can be expressed in the simple form |ψ⟩ = |ψ₁⟩ ⊗ |ψ₂⟩, in any specific basis, are very particular and they are called separable states.³⁸ In general, any arbitrary state should be described as |ψ⟩ = ∑_i,j|ψ_i⟩ ⊗ |ψ_j⟩ (or ρ = ∑_i,jρ_i ⊗ ρ_j for mixed states). Nonseparable states are called entangled states.

Now, that we know how to compose systems, but we can be interested in going the other way around. If we have a system belonging to a bipartite Hilbert space in the form

H = H_{a} \otimes H_{b}

⁠, we can be interested in studying some properties of the subsystem corresponding to one of the subspaces. To do so, we define the reduced density matrix. If the state of our system is described by a density matrix ρ, the reduced density operator of the subsystem a is defined by the operator,

ρ_{a} \equiv {Tr}_{b} [ρ],

(18)

where Tr_b is the partial trace over subspace b and it is defined as²²

{Tr}_{b} [\sum_{i, j, k, l} | a_{i} (⟩ ⟨) a_{j} | \otimes | b_{k} (⟩ ⟨) b_{l} |] \equiv \sum_{i, j} | a_{i} (⟩ ⟨) a_{j} | Tr [\sum_{k, l} | b_{k} (⟩ ⟨) b_{l} |] .

(19)

The concepts of reduced density matrix and partial trace are essential in the study of open quantum systems. If we want to calculate the equation of motions of a system affected by an environment, we should trace out this environment and deal only with the reduced density matrix of the system. This is the main idea of the theory of open quantum systems.Example 4: Two two-level atoms.

If we have two two-level systems, the total Hilbert space is given by

H = H_{2} \otimes H_{2}

⁠. A basis of this Hilbert space would be given by the set

\{) | 00 (⟩ \equiv | 0 {(⟩}_{1} \otimes | 0 {(⟩}_{2}, | 01 (⟩ \equiv | 0 {(⟩}_{1} \otimes | 1 {(⟩}_{2}, | 10 (⟩ \equiv | 1 {(⟩}_{1} \otimes | 0 {(⟩}_{2}, | 11 (⟩ \equiv | 1 {(⟩}_{1} \otimes | 1 {(⟩}_{2} \})

⁠. If both systems are in their ground state, we can describe the total state by the separable vector,

| ψ {(⟩}_{G} = | 00 (⟩ .

(20)

A more complex, but still separable, state can be formed if both systems are in superposition,

\begin{matrix} | ψ {(⟩}_{S} & = & \frac{1}{\sqrt{2}} (| 0 {(⟩}_{1} + | 1 {(⟩}_{1}) \otimes \frac{1}{\sqrt{2}} (| 0 {(⟩}_{2} + | 1 {(⟩}_{2}) \\ = & \frac{1}{2} (| 00 (⟩ + | 10 (⟩ + | 01 (⟩ + | 11 (⟩) . \end{matrix}

(21)

An entangled state would be

| ψ {(⟩}_{E} = \frac{1}{\sqrt{2}} (| 00 (⟩ + | 11 (⟩) .

(22)

This state cannot be separated into a direct product of each subsystem. If we want to obtain a reduced description of subsystem 1 (or 2), we have to use the partial trace. To do so, we need first to calculate the density matrix corresponding to the pure state |ψ⟩_E,

ρ_{E} = | ψ (⟩ ⟨) ψ |_{E} = \frac{1}{2} (| 00 (⟩ ⟨) 00 | + | 00 (⟩ ⟨) 11 | + | 11 (⟩ ⟨) 00 | + | 11 (⟩ ⟨) 11 |) .

(23)

We can now calculate the reduced density matrix of the subsystem 1 by using the partial trace,

ρ_{E}^{(1)} = ⟨) 0 |_{2} ρ_{E} | 0 {(⟩}_{2} + ⟨) 1 |_{2} ρ_{E} | 1 {(⟩}_{2} = \frac{1}{2} (| 00 (⟩ ⟨) 00 |_{1} + | 11 (⟩ ⟨) 11 |_{2}) .

(24)

From this reduced density matrix, we can calculate all the measurement statistics of subsystem 1.End of example.

IV. THE FOCK-LIOUVILLE HILBERT SPACE. THE LIOUVILLE SUPEROPERATOR

In this section, we revise a useful framework for both analytical and numerical calculations. It is clear that some linear combinations of density matrices are valid density matrices (as long as they preserve positivity and trace 1). Because of that, we can create a Hilbert space of density matrices just by defining a scalar product.³⁹ This allows us to define a linear space of matrices, converting the matrices effectively into vectors (ρ → |ρ⟩⟩). This is called Fock-Liouville space (FLS). The usual definition of the scalar product of matrices ϕ and ρ is defined as $⟨ ⟨ ϕ | ρ ⟩ ⟩ \equiv Tr [ϕ^{†} ρ]$ ⁠. The Liouville superoperator from Eq. (13) is now an operator acting on the Hilbert space of density matrices. The main utility of the FLS is to allow the matrix representation of the evolution operator.Example 5: Time evolution of a two-level atom.

The density matrix of our system (25) can be expressed in the FLS as

| ρ (⟩ (⟩ = (\begin{matrix} ρ_{00} \\ ρ_{01} \\ ρ_{10} \\ ρ_{11} \end{matrix}) .

(25)

The time evolution of a mixed state is given by the von Neumann equation (13). The Liouvillian superoperator can now be expressed as a matrix,

\tilde{L} = (\begin{matrix} 0 & i Ω & - i Ω & 0 \\ i Ω & i E & 0 & - i Ω \\ - i Ω & 0 & - i E & i Ω \\ 0 & - i Ω & i Ω & 0 \end{matrix}),

(26)

where each row is calculated just by observing the output of the operation $- i [H, ρ]$ in the computational basis of the density matrices space. The time evolution of the system now corresponds to the matrix equation $\frac{d | ρ (⟩ (⟩}{d t} = \tilde{L} | ρ (⟩ (⟩$ ⁠, which in matrix notation would be

(\begin{matrix} {\dot{ρ}}_{00} \\ {\dot{ρ}}_{01} \\ {\dot{ρ}}_{10} \\ {\dot{ρ}}_{11} \end{matrix}) = (\begin{matrix} 0 & i Ω & - i Ω & 0 \\ i Ω & i E & 0 & - i Ω \\ - i Ω & 0 & - i E & i Ω \\ 0 & - i Ω & i Ω & 0 \end{matrix}) (\begin{matrix} ρ_{00} \\ ρ_{01} \\ ρ_{10} \\ ρ_{11} \end{matrix}) .

(27)

End of example.

V. CPT-MAPS AND THE LINDBLAD MASTER EQUATION

A. Completely positive maps

The problem we want to study is to find the most general Markovian transformation set between density matrices. Until now, we have seen that quantum systems can evolve in two ways by a coherent evolution (postulate 3) and by collapsing after a measurement (postulate 2). Many efforts have been made to unify these two ways of evolving,¹⁶ without giving a definite answer so far. It is reasonable to ask what is the most general transformation that can be performed in a quantum system and what is the dynamical equation that describes this transformation.

We are looking for maps that transform density matrices into density matrices. We define $ρ (H)$ as the space of all density matrices in the Hilbert space $H$ ⁠. Therefore, we are looking for a map of this space onto itself, $V : ρ (H) \to ρ (H)$ ⁠. To ensure that the output of the map is a density matrix, this should fulfill the following properties:

Trace preserving. $Tr [V A] = Tr [A]$ ⁠, $\forall A \in O (H)$ ⁠.
Completely positive (see below).

Any map that fulfills these two properties is called a completely positive and trace-preserving map (CPT-maps). The first property is quite apparent, and it does not require more thinking. The second one is a little more complicated, and it requires an intermediate definition.

Definition 1.

A map $V$ is positive iff $\forall A \in B (H)$ s.t. $A \geq 0 \Rightarrow V A \geq 0$ .

This definition is based on the idea that, as density matrices are positive, any physical map should transform positive matrices into positive matrices. One could naively think that this condition must be sufficient to guarantee the physical validity of a map; it is not. As we know, there exist composite systems, and our density matrix could be the partial trace of a more complicated state. Because of that, we need to impose a more general condition.

Definition 2.

A map $V$ is completely positive iff $\forall n \in N$ , $V \otimes 1_{n}$ is positive.

To prove that not all positive maps are completely positive, we need a counterexample. A canonical example of an operation that is positive but fails to be completely positive is the matrix transposition. If we have a Bell state in the form

| ψ_{B} (⟩ = \frac{1}{\sqrt{2}} (| 01 (⟩ + | 10 (⟩)

⁠, its density matrix can be expressed as

\begin{matrix} ρ_{B} = \frac{1}{2} (| 0 (⟩ ⟨) 0 | \otimes | 1 (⟩ ⟨) 1 | + | 1 (⟩ ⟨) 1 | \otimes | 0 (⟩ ⟨) 0 | + | 0 (⟩ ⟨) 1 | \otimes | 1 (⟩ ⟨) 0 | + | 1 (⟩ ⟨) 0 | \otimes | 0 (⟩ ⟨) 1 |), \end{matrix}

(28)

with a matrix representation

\begin{matrix} ρ_{B} & = & \frac{1}{2} \{(\begin{matrix} 1 & 0 \\ 0 & 0 \end{matrix}) \otimes (\begin{matrix} 0 & 0 \\ 0 & 1 \end{matrix}) + (\begin{matrix} 0 & 0 \\ 0 & 1 \end{matrix}) \otimes (\begin{matrix} 1 & 0 \\ 0 & 0 \end{matrix})) \\ (\otimes (\begin{matrix} 0 & 0 \\ 1 & 0 \end{matrix}) \otimes (\begin{matrix} 0 & 1 \\ 0 & 0 \end{matrix}) + (\begin{matrix} 0 & 1 \\ 0 & 0 \end{matrix}) \otimes (\begin{matrix} 0 & 0 \\ 1 & 0 \end{matrix})\} . \end{matrix}

(29)

A little algebra shows that the full form of this matrix is

ρ_{B} = (\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}),

(30)

and it is positive.

It is easy to check that the transformation

1 \otimes T_{2}

⁠, meaning that we transpose the matrix of the second subsystem leads to a nonpositive matrix

\begin{matrix} (1 \otimes T_{2}) ρ_{B} & = & \frac{1}{2} \{(\begin{matrix} 1 & 0 \\ 0 & 0 \end{matrix}) \otimes (\begin{matrix} 0 & 1 \\ 0 & 0 \end{matrix}) + (\begin{matrix} 0 & 0 \\ 0 & 1 \end{matrix}) \otimes (\begin{matrix} 0 & 0 \\ 1 & 0 \end{matrix})) \\ (\otimes (\begin{matrix} 0 & 0 \\ 1 & 0 \end{matrix}) \otimes (\begin{matrix} 0 & 0 \\ 1 & 0 \end{matrix}) + (\begin{matrix} 0 & 0 \\ 0 & 1 \end{matrix}) \otimes (\begin{matrix} 0 & 1 \\ 0 & 0 \end{matrix})\} . \end{matrix}

(31)

The total matrix is

(1 \otimes T_{2}) ρ_{B} = (\begin{matrix} 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 \end{matrix}),

(32)

with −1 as an eigenvalue. This example illustrates how the nonseparability of quantum mechanics restricts the operations we can perform in a subsystem. By imposing these two conditions, we can derive a unique master equation as the generator of any possible Markovian CPT-map.

B. Derivation of the Lindblad equation from microscopic dynamics

The most common derivation of the Lindblad master equation is based on open quantum theory. The Lindblad equation is then an effective motion equation for a subsystem that belongs to a more complicated system. This derivation can be found in several textbooks such as Breuer and Petruccione’s² as well as Gardiner and Zoller’s.¹ Here, we follow the derivation presented in Ref. 26. Our initial point is displayed in Fig. 3. A total system belonging to a Hilbert space $H_{T}$ is divided into our system of interest, belonging to a Hilbert space $H$ ⁠, and the environment living in $H_{E}$ ⁠.

FIG. 3.

View large Download slide

A total system (belonging to a Hilbert space $H_{T}$ ⁠, with states described by density matrices ρ_T, and with dynamics determined by a Hamiltonian H_T) divided into the system of interest, “system,” and the environment.

The evolution of the total system is given by the von Neumann equation (13),

\dot{ρ_{T}} (t) = - i [H_{T}, ρ_{T} (t)] .

(33)

As we are interested in the dynamics of the system, without the environment, we trace over the environment degrees of freedom to obtain the reduced density matrix of the system ρ(t) = Tr_E[ρ_T]. To separate the effect of the total Hamiltonian in the system and the environment, we divide it in the form $H_{T} = H_{S} \otimes 1_{E} + 1_{S} \otimes H_{E} + α H_{I}$ ⁠, with $H \in H$ ⁠, $H_{E} \in H_{E}$ ⁠, and $H_{I} \in H_{T}$ ⁠, and being α a measure of the strength of the system-environment interaction. Therefore, we have a part acting on the system, a part acting on the environment, and the interaction term. Without losing any generality, the interaction term can be decomposed in the following way:

H_{I} = \sum_{i} S_{i} \otimes E_{i},

(34)

with $S_{i} \in B (H)$ and $E_{i} \in B (H_{E})$ ⁠.⁴⁰

To better describe the dynamics of the system, it is useful to work in the interaction picture (see Ref. 24 for a detailed explanation about Schrödinger, Heisenberg, and interaction pictures). In the interaction picture, density matrices evolve with time due to the interaction Hamiltonian, while operators evolve with the system and environment Hamiltonian. An arbitrary operator $O \in B (H_{T})$ is represented in this picture by the time-dependent operator Ô(t), and its time evolution is

Ô (t) = e^{i (H + H_{E}) t} O e^{- i (H + H_{E}) t} .

(35)

The time evolution of the total density matrix is given in this picture by

\frac{d {\hat{ρ}}_{T} (t)}{d t} = - i α [Ĥ_{I} (t), {\hat{ρ}}_{T} (t)] .

(36)

This equation can be easily integrated to give

{\hat{ρ}}_{T} (t) = {\hat{ρ}}_{T} (0) - i α \int_{0}^{t} d s [Ĥ_{I} (s), {\hat{ρ}}_{T} (s)] .

(37)

By this formula, we can obtain the exact solution, but it still has the complication of calculating an integral in the total Hilbert space. It is also a troublesome fact that the state $\tilde{ρ} (t)$ depends on the integration of the density matrix in all previous time. To avoid that, we can introduce Eq. (37) into Eq. (36) giving

\begin{matrix} \frac{d {\hat{ρ}}_{T} (t)}{d t} = - i α [Ĥ_{I} (t), {\hat{ρ}}_{T} (0)] - α^{2} \int_{0}^{t} d s [Ĥ_{I} (t), [Ĥ_{I} (s), {\hat{ρ}}_{T} (s)]] . \end{matrix}

(38)

By applying this method one more time, we obtain

\begin{matrix} \frac{d {\hat{ρ}}_{T} (t)}{d t} & = & - i α [Ĥ_{I} (t), {\hat{ρ}}_{T} (0)] \\ - α^{2} \int_{0}^{t} d s [Ĥ_{I} (t), [Ĥ_{I} (s), {\hat{ρ}}_{T} (t)]] + O (α^{3}) . \end{matrix}

(39)

After this substitution, the integration of the previous states of the system is included only in the terms that are O(α³) or higher. At this moment, we perform our first approximation by considering that the strength of the interaction between the system and the environment is small. Therefore, we can avoid high-orders in Eq. (39). Under this approximation, we have

\begin{matrix} \frac{d {\hat{ρ}}_{T} (t)}{d t} = - i α [Ĥ_{I} (t), {\hat{ρ}}_{T} (0)] - α^{2} \int_{0}^{t} d s [Ĥ_{I} (t), [Ĥ_{I} (s), {\hat{ρ}}_{T} (t)]] . \end{matrix}

(40)

We are interested in finding an equation of motion for ρ, so we trace over the environment degrees of freedom,

\begin{matrix} \frac{d \hat{ρ} (t)}{d t} & = & {Tr}_{E} [\frac{d {\hat{ρ}}_{T} (t)}{d t}] = - i α {Tr}_{E} [Ĥ_{I} (t), {\hat{ρ}}_{T} (0)] \\ - α^{2} \int_{0}^{t} d s {Tr}_{E} [Ĥ_{I} (t), [Ĥ_{I} (s), {\hat{ρ}}_{T} (t)]] . \end{matrix}

(41)

This is not a closed time-evolution equation for $\hat{ρ} (t)$ because the time derivative still depends on the full density matrix ${\hat{ρ}}_{T} (t)$ ⁠. To proceed, we need to make two more assumptions. First, we assume that at t = 0 the system and the environment have a separable state in the form ρ_T(0) = ρ(0) ⊗ ρ_E(0). This means that there are not correlations between the system and the environment. This may be the case if the system and the environment have not interacted at previous times or if the correlations between them are short-lived. Second, we assume that the initial state of the environment is thermal, meaning that it is described by a density matrix in the form $ρ_{E} (0) = \exp (- H_{E} / T) / Tr [\exp (- H_{E} / T)]$ ⁠, being T the temperature and taking the Boltzmann constant as k_B = 1. By using these assumptions, and the expansion of H_I(34), we can calculate an expression for the first element of the rhs of Eq. (41),

\begin{matrix} {Tr}_{E} [Ĥ_{I} (t), {\hat{ρ}}_{T} (0)] & = & \sum_{i} () Ŝ_{i} (t) \hat{ρ} (0) {Tr}_{E} [Ê_{i} (t) {\hat{ρ}}_{E} (0)] \\ - \hat{ρ} (0) Ŝ_{i} (t) {Tr}_{E} [{\hat{ρ}}_{E} (0) Ê_{i} (t)])) . \end{matrix}

(42)

To calculate the explicit value of this term, we may use that $⟨E_{i}⟩ = Tr [E_{i} ρ_{E} (0)] = 0$ for all values of i. This looks like a strong assumption, but it is not. If our total Hamiltonian does not fulfill it, we can always rewrite it as $H_{T} = (H + α \sum_{i} ⟨E_{i}⟩ S_{i}) + H_{E} + α H_{i}^{'}$ ⁠, with $H_{i}^{'} = \sum_{i} S_{i} \otimes (E_{i} - ⟨E_{i}⟩)$ ⁠. It is clear that now $⟨E_{i}^{'}⟩ = 0$ ⁠, with $E_{i}^{'} = E_{i} - ⟨E_{i}⟩$ ⁠, and the system Hamiltonian is changed just by the addition of an energy shift that does not affect the system dynamics. Because of that, we can assume that $⟨E_{i}⟩ = 0$ for all i. Using the cyclic property of the trace, it is easy to prove that the term of Eq. (42) is equal to zero and the equation of motion (41) reduces to

\dot{\hat{ρ}} (t) = - α^{2} \int_{0}^{t} d s {Tr}_{E} [Ĥ_{I} (t), [Ĥ_{I} (s), {\hat{ρ}}_{T} (t)]] .

(43)

This equation still includes the entire state of the system and environment. To unravel the system from the environment, we have to make a more restrictive assumption. As we are working in the weak coupling regime, we may suppose that the system and the environment are noncorrelated during all the time evolution. Of course, this is only an approximation. Due to the interaction Hamiltonian, some correlations between system and environment are expected to appear. On the other hand, we may assume that the time scales of correlation (τ_corr) and relaxation of the environment (τ_rel) are much smaller than the typical system time scale (τ_sys), as the coupling strength is very small (α≪). Therefore, under this strong assumption, we can assume that the environment state is always thermal and is decoupled from the system state, ${\hat{ρ}}_{T} (t) = \hat{ρ} (t) \otimes {\hat{ρ}}_{E} (0)$ ⁠. Equation (43) then transforms into

\dot{\hat{ρ}} (t) = - α^{2} \int_{0}^{t} d s {Tr}_{E} [Ĥ_{I} (t), [Ĥ_{I} (s), \hat{ρ} (t) \otimes {\hat{ρ}}_{E} (0)]] .

(44)

The equation of motion is now independent for the system and local in time. It is still non-Markovian, as it depends on the initial state preparation of the system. We can obtain a Markovian equation by realizing that the kernel in the integration decays fast enough and that we can extend the upper limit of the integration to infinity with no real change in the outcome. By doing so, and by changing the integral variable to s → t − s, we obtain the famous Redfield equation⁴¹

\dot{\hat{ρ}} (t) = - α^{2} \int_{0}^{\infty} d s {Tr}_{E} [Ĥ_{I} (t), [Ĥ_{I} (s - t), \hat{ρ} (t) \otimes {\hat{ρ}}_{E} (0)]] .

(45)

It is known that this equation does not warrant the positivity of the map, and it sometimes gives rise to density matrices that are nonpositive. To ensure complete positivity, we need to perform one further approximation, the rotating wave approximation. To do so, we need to use the spectrum of the superoperator $\tilde{H} A \equiv [H, A]$ ⁠, $\forall A \in B (H)$ ⁠. The eigenvectors of this superoperator form a complete basis of space $B (H)$ ⁠, and, therefore, we can expand the system-environment operators from Eq. (34) in this basis,

S_{i} = \sum_{ω} S_{i} (ω),

(46)

where the operators S_i(ω) fulfill

[H, S_{i} (ω)] = - ω S_{i} (ω),

(47)

being ω the eigenvalues of $\tilde{H}$ ⁠. It is easy to take also the Hermitian conjugated

[H, S_{i}^{†} (ω)] = ω S_{i}^{†} (ω) .

(48)

To apply this decomposition, we need to change back to the Schrödinger picture for the term of the interaction Hamiltonian acting on the system’s Hilbert space. This is done by the expression $Ŝ$ _k = e^itHS_ke^−itH. By using the eigenexpansion (46), we arrive to

{\tilde{H}}_{i} (t) = \sum_{k, ω} e^{- i ω t} S_{k} (ω) \otimes {\tilde{E}}_{k} (t) = \sum_{k, ω} e^{i ω t} S_{k}^{†} (ω) \otimes {\tilde{E}}_{k}^{†} (t) .

(49)

To combine this decomposition with Redfield equation (45), we first may expand the commutators,

\begin{matrix} \dot{\hat{ρ}} (t) & = & - α^{2} Tr [) \int_{0}^{\infty} d s Ĥ_{I} (t) Ĥ_{I} (t - s) \hat{ρ} (t) \otimes {\hat{ρ}}_{E} (0) \\ - \int_{0}^{\infty} d s Ĥ_{I} (t) \hat{ρ} (t) \otimes {\hat{ρ}}_{E} (0) Ĥ_{I} (t - s) \\ - \int_{0}^{\infty} d s Ĥ_{I} (t - s) \hat{ρ} (t) \otimes {\hat{ρ}}_{E} (0) Ĥ_{I} (t) \\ + \int_{0}^{\infty} d s \hat{ρ} (t) \otimes {\hat{ρ}}_{E} (0) Ĥ_{I} (t - s) Ĥ_{I} (t)]) . \end{matrix}

(50)

We now apply the eigenvalue decomposition in terms of S_k(ω) for Ĥ_I(t − s) and in terms of $S_{k}^{†} (ω^{'})$ for Ĥ_I(t). By using the permutation property of the trace and the fact that $[H_{E}, ρ_{E} (0)] = 0$ ⁠, and after some nontrivial algebra, we obtain

\begin{matrix} \dot{\hat{ρ}} (t) & = & \sum_{\begin{array}{c} ω, ω^{'} \\ k, l \end{array}} () e^{i (ω^{'} - ω) t} Γ_{k l} (ω) [S_{l} (ω) \hat{ρ} (t), S_{k}^{†} (ω^{'})] \\ + e^{i (ω - ω^{'}) t} Γ_{l k}^{*} (ω^{'}) [S_{l} (ω), \hat{ρ} (t) S_{k}^{†} (ω^{'})])), \end{matrix}

(51)

where the effect of the environment has been absorbed into the factors,

Γ_{k l} (ω) \equiv \int_{0}^{\infty} d s e^{i ω s} {Tr}_{E} [{\tilde{E}}_{k}^{†} (t) {\tilde{E}}_{l} (t - s) ρ_{E} (0)],

(52)

where we are writing the environment operators of the interaction Hamiltonian in the interaction picture (⁠ $Ê_{l} (t) = e^{i H_{E} t} E_{l} e^{- i H_{E} t}$ ⁠). At this point, we can already perform the rotating wave approximation. By considering the time-dependency on Eq. (51), we conclude that the terms with $|ω - ω^{'}| ≫ α^{2}$ will oscillate much faster than the typical time scale of the system evolution. Therefore, they do not contribute to the evolution of the system. In the low-coupling regime (α → 0), we can consider that only the resonant terms, ω = ω′, contribute to the dynamics and remove all the others. By applying this approximation to Eq. (51) reduces to

\begin{matrix} \dot{\hat{ρ}} (t) = \sum_{\begin{array}{c} ω \\ k, l \end{array}} (Γ_{k l} (ω) [S_{l} (ω) \hat{ρ} (t), S_{k}^{†} (ω)] + Γ_{l k}^{*} (ω) [S_{l} (ω), \hat{ρ} (t) S_{k}^{†} (ω)]) . \end{matrix}

(53)

To divide the dynamics into Hamiltonian and non-Hamiltonian, we now decompose the operators Γ_kl into Hermitian and non-Hermitian parts, $Γ_{k l} (ω) = \frac{1}{2} γ_{k l} (ω) + i π_{k l}$ ⁠, with

\begin{aligned} π_{k l} (ω) & \equiv \frac{- i}{2} (Γ_{k l} (ω) - Γ_{k l}^{*} (ω)), \\ γ_{k l} (ω) & \equiv Γ_{k l} (ω) + Γ_{k l}^{*} (ω) = \int_{- \infty}^{\infty} d s e^{i ω s} Tr [Ê_{k}^{†} (s) E_{l} {\hat{ρ}}_{E} (0)] . \end{aligned}

(54)

By these definitions, we can separate the Hermitian and non-Hermitian parts of the dynamics and we can transform back to the Schrödinger picture,

\begin{matrix} \dot{ρ} (t) & = & - i [H + H_{L s}, ρ (t)] + \sum_{\begin{array}{c} ω \\ k, l \end{array}} γ_{k l} (ω) \\ \times (S_{l} (ω) ρ (t) S_{k}^{†} (ω) - \frac{1}{2} \{S_{k}^{†} S_{l} (ω), ρ (t)\}) . \end{matrix}

(55)

The Hamiltonian dynamics now is influenced by a term $H_{L s} = \sum_{ω, k, l} π_{k l} (ω) S_{k}^{†} (ω) S_{l} (ω)$ ⁠. This is usually called a Lamb shift Hamiltonian, and its role is to renormalize the system energy levels due to the interaction with the environment. Equation (55) is the first version of the Markovian master equation, but it is not in the Lindblad form yet.

It can be easily proved that the matrix formed by the coefficients γ_kl(ω) is positive as they are the Fourier’s transform of a positive function $(Tr [Ê_{k}^{†} (s) E_{l} {\hat{ρ}}_{E} (0)])$ ⁠. Therefore, this matrix can be diagonalized. This means that we can find a unitary operator, O, s.t.

O γ (ω) O^{†} = (\begin{matrix} d_{1} (ω) & 0 & \dots & 0 \\ 0 & d_{2} (ω) & \dots & 0 \\ ⋮ & ⋮ & ⋱ & 0 \\ 0 & 0 & \dots & d_{N} (ω) \end{matrix}) .

(56)

We can now write the master equation in a diagonal form,

\begin{matrix} \dot{ρ} (t) & = & - i [H + H_{L s}, ρ (t)] \\ + \sum_{i, ω} (L_{i} (ω) ρ (t) L_{i}^{†} (ω) - \frac{1}{2} \{L_{i}^{†} L_{i} (ω), ρ (t)\}) \equiv L ρ (t) . \end{matrix}

(57)

This is the celebrated Lindblad (or Lindblad-Gorini-Kossakowski-Sudarshan) master equation. In the simplest case, there will be only one relevant frequency ω and the equation can be further simplified to

\begin{matrix} \dot{ρ} (t) & = & - i [H + H_{L s}, ρ (t)] + \sum_{i} (L_{i} ρ (t) L_{i}^{†} - \frac{1}{2} \{L_{i}^{†} L_{i}, ρ (t)\}) \equiv L ρ (t) . \end{matrix}

(58)

The operators L_i are usually referred to as jump operators.

C. Derivation of the Lindblad equation as a CPT generator

The second way of deriving the Lindblad equation comes from the following question: What is the most general (Markovian) way of mapping density matrix onto density matrices? This is usually the approach from quantum information researchers that look for general transformations of quantum systems. We analyze this problem following mainly Ref. 28.

To start, we need to know what is the form of a general CPT-map.

Lemma 2.

Any map $V : B (H) \to B (H)$ that can be written in the form $V ρ = V^{†} ρ V$ with $V \in B (H)$ is positive.

The proof of the lemma requires a little algebra and a known property of normal matrices.

Proof.

If $ρ \geq 0 \Rightarrow ρ = A^{†} A$ ⁠, with $A \in B (H)$ ⁠. Therefore, $V ρ = V^{†} ρ V \Rightarrow ⟨) ψ | V^{†} ρ V | ψ (⟩ = ⟨) ψ | V^{†} A^{†} A V | ψ (⟩ = | | A V | ψ (⟩ | | \geq 0$ ⁠. Therefore, if $ρ$ is positive, the output of the map is also positive.End of the proof.

This is a sufficient condition for the positivity of a map, but it is not necessary. It could happen that there are maps that cannot be written in this form, but they are still positive. To go further, we need a more general condition and this comes in the form of the following theorem:

Theorem 1.

Choi’s theorem.

A linear map

V : B (H) \to B (H)

is completely positive iff it can be expressed as

V ρ = \sum_{i} V_{i}^{†} ρ V_{i},

(59)

with

V_{i} \in B (H)

⁠.

The proof of this theorem requires some algebra.

Proof.

The “if” implication is a trivial consequence of Lemma 2. To prove the converse, we need to extend the dimension of our system by the use of an auxiliary system. If d is the dimension of the Hilbert space of pure states, $H$ ⁠, we define a new Hilbert space of the same dimension $H_{A}$ ⁠.

We define a maximally entangled pure state in the bipartition

H_{A} \otimes H

in the way,

| Γ (⟩ \equiv \sum_{i = 0}^{d - 1} | i {(⟩}_{A} \otimes | i (⟩,

(60)

being

\{| i (⟩\}

and

\{| i {(⟩}_{A}\}

arbitrary orthonormal bases for

H

and

H_{A}

⁠, respectively.

We can extend the action of our original map

V

that acts on

B (H)

to our extended Hilbert space by defining the map

V_{2} : B (H_{A}) \otimes B (H) \to B (H_{A}) \otimes B (H)

as

V_{2} \equiv 1_{B (H_{A})} \otimes V .

(61)

Note that the idea behind this map is to leave the auxiliary subsystem invariant while applying the original map to the original system. This map is positive because

V

is completely positive. This may appear trivial, but as it has been explained before complete positivity is a more restrictive property than positivity, and we are looking for a condition to ensure complete positivity.

We can now apply the extended map to the density matrix corresponding to the maximally entangled state , obtaining

V_{2} | Γ (⟩ ⟨) Γ | = \sum_{i, j = 0}^{d - 1} | i (⟩ ⟨) j | \otimes V | i (⟩ ⟨) j | .

(62)

Now, we can use the maximal entanglement of the state |Γ⟩ to relate the original map

V

and the action

V_{2} | Γ (⟩ ⟨) Γ |

by taking the matrix elements with respect to

H_{A}

⁠,

V | i (⟩ ⟨) j | = ⟨) i |_{A} (V_{2} | Γ (⟩ ⟨) Γ |) | j {(⟩}_{A} .

(63)

To relate this operation to the action of the map to an arbitrary vector

| ψ (⟩ \in H_{A} \otimes H

⁠, we can expand it in this basis as

| ψ (⟩ = \sum_{i = 0}^{d - 1} \sum_{j = 0}^{d - 1} α_{i j} | i {(⟩}_{A} \otimes | j (⟩ .

(64)

We can also define an operator

V_{| ψ (⟩} \in B (H)

s.t. it transforms |Γ⟩ into |ψ⟩. Its explicit action would be written as

\begin{matrix} (1_{A} \otimes V_{| ψ (⟩}) | Γ (⟩ & = & \sum_{i, j = 0}^{d - 1} α_{i j} (1_{A} \otimes | j (⟩ ⟨) i |) (\sum_{k = 0}^{d - 1} | k (⟩ \otimes | k (⟩) \\ = & \sum_{i, j, k = 0}^{d - 1} α_{i j} (| k (⟩ \otimes | j (⟩) ⟨ i | k ⟩ = \sum_{i, j, k = 0}^{d - q} α_{i j} (| k (⟩ \otimes | j (⟩) δ_{i, k} \\ = & \sum_{i, j = 0}^{d - 1} α_{i j} | i (⟩ \otimes | j (⟩ = | ψ (⟩ . \end{matrix}

(65)

At this point, we have related the vectors in the extended space

H_{A} \otimes H

to operators acting on

H

⁠. This can only be done because the vector |Γ⟩ is maximally entangled. We go now back to our extended map

V_{2}

⁠. Its action on |Γ⟩⟨Γ| is given by Eq. , and as it is a positive map, it can be expanded as

V_{2} (| Γ (⟩ ⟨) Γ |) = \sum_{l = 0}^{d^{2} - 1} | v_{l} (⟩ ⟨) v_{l} |,

(66)

with

| v_{l} (⟩ \in H_{A} \otimes H

⁠. The vectors |

v

_l⟩ can be related to operators in

H

as in Eq. ,

| v_{l} (⟩ = (1_{A} \otimes V_{l}) | Γ (⟩ .

(67)

Based on this result, we can calculate the product of an arbitrary vector

| i {(⟩}_{A} \in H_{A}

⁠, with |

v

_l⟩,

⟨) i |_{A} | v_{l} (⟩ = ⟨) i |_{A} (1_{A} \otimes V_{l}) | Γ (⟩ = V_{l} \sum_{k = 0}^{d - 1} {⟨ i | k ⟩}_{A} \otimes | k (⟩ .

(68)

This is the last ingredient we need for the proof.

We come back to the original question, we want to characterize the map

V

⁠. We do so by applying it to an arbitrary basis element |i⟩⟨j| of

B (H)

⁠,

\begin{matrix} V (| i (⟩ ⟨) j |) & = & (⟨) i |_{A} \otimes 1_{A}) V_{2} (| Γ (⟩ ⟨) Γ |) (| j {(⟩}_{A} \otimes 1_{A}) \\ = & (⟨) i |_{A} \otimes 1_{A}) [\sum_{l = 0}^{d^{2} - 1} | v_{l} (⟩ ⟨) v_{l} |] (| j {(⟩}_{A} \otimes 1_{A}) \\ = & \sum_{l = 0}^{d^{2} - 1} [(⟨) i |_{A} \otimes 1_{A}) | v_{l} (⟩] [⟨) v_{l} | (| j {(⟩}_{A} \otimes 1_{A})] \\ = & \sum_{l = 0}^{d^{2} - 1} V_{l}^{†} | i (⟩ ⟨) j | V_{l} . \end{matrix}

(69)

As |i⟩⟨j| is an arbitrary element of a basis, any operator can be expanded in this basis. Therefore, it is straightforward to prove that

V ρ = \sum_{l}^{d^{2} - l} V_{l} ρ V_{l}^{†} .

End of the proof.

Thanks to Choi’s theorem, we know the general form of CP-maps, but there is still an issue to address. As density matrices should have trace one, we need to require any physical maps to be also trace-preserving. This requirement gives as a new constraint that completely defines all CPT-maps. This requirement comes from the following theorem:

Theorem 2.

Choi-Kraus’ theorem.

A linear map

V : B (H) \to B (H)

is completely positive and trace-preserving iff it can be expressed as

V ρ = \sum_{l} V_{l}^{†} ρ V_{l},

(70)

with

V_{l} \in B (H)

fulfilling

\sum_{l} V_{l} V_{l}^{†} = 1_{H} .

(71)

Proof.

We have already proved that this is a completely positive map and we only need to prove that it is also trace-preserving and that all trace preserving-maps fulfill Eq. . The “if” proof is quite simple by applying the cyclic permutations and linearity properties of the trace operator,

Tr [V ρ] = Tr [\sum_{l = 1}^{d^{2} - 1} V_{l}^{†} ρ V_{l}] = Tr [(\sum_{l = 1}^{d^{2} - 1} V_{l} V_{l}^{†}) ρ] = Tr [ρ] .

(72)

We have to prove also that any map in the form is trace-preserving only if the operators V_l fulfill . We start by stating that if the map is trace-preserving by applying it to an any arbitrary element of a basis of

B (H)

⁠, we should obtain

Tr [V (| i (⟩ ⟨) j |)] = Tr [| i (⟩ ⟨) j |] = δ_{i, j} .

(73)

As the map has a form given by , we can calculate this same trace in an alternative way,

\begin{matrix} Tr [V (| i (⟩ ⟨) j |)] & = & Tr [\sum_{l = 1}^{d^{2} - 1} V_{l}^{†} | i (⟩ ⟨) j | V_{l}] = Tr [\sum_{l = 1}^{d^{2} - 1} V_{l} V_{l}^{†} | i (⟩ ⟨) j |] \\ = & \sum_{k} ⟨) k | (\sum_{l = 1}^{d^{2} - 1} V_{l} V_{l}^{†} | i (⟩ ⟨) j |) | k (⟩ = ⟨) j | (\sum_{l = 1}^{d^{2} - 1} V_{l} V_{l}^{†}) | i (⟩, \end{matrix}

(74)

where

\{| k (⟩\}

is an arbitrary basis of

H

⁠. As both equalities should be right, we obtain

⟨) j | (\sum_{l = 1}^{d^{2} - 1} V_{l} V_{l}^{†}) | i (⟩ = δ_{i, j},

(75)

and therefore, condition should be fulfilled.End of the proof.

Operators V_i of a map fulfilling condition (71) are called Krauss operators. Because of that, sometimes CPT-maps are also called Krauss maps, especially when they are presented as a collection of Krauss operators. Both concepts are ubiquitous in quantum information science. Krauss operators can also be time-dependent as long as they fulfill relation (71) for all times.

At this point, we already know the form of CPT-maps, but we do not have a master equation that is a continuous set of differential equations. This means that we know how to perform an arbitrary operation in a system, but we do not have an equation to describe its time evolution. To do so, we need to find a time-independent generator

L

such that

\frac{d}{d t} ρ (t) = L ρ (t),

(76)

and therefore our CPT-map could be expressed as

V (t) = e^{L t}

⁠. The following calculation is about founding the explicit expression of

L

⁠. We start by choosing an orthonormal basis of the bounded space of operators

B (H)

⁠,

{\{F_{i}\}}_{i = 1}^{d^{2}}

⁠. To be orthonormal, it should satisfy the following condition:

⟨ ⟨ F_{i} | F_{j} ⟩ ⟩ \equiv Tr [F_{i}^{†} F_{j}] = δ_{i, j} .

(77)

Without any loss of generality, we select one of the elements of the basis to be proportional to the identity,

F_{d^{2}} = \frac{1}{\sqrt{d}} 1_{H}

⁠. It is trivial to prove that the norm of this element is one, and it is easy to see from Eq. that all the other elements of the basis should have trace zero,

Tr [F_{i}] = 0 \forall i = 1, \dots, d^{2} - 1 .

(78)

The closure relation of this basis is

1_{B (H)} = \sum_{i} | F_{i} (⟩ (⟩ ⟨) ⟨) F_{i} |

⁠. Therefore, the Krauss operators can be expanded in this basis by using the Fock-Liouville notation,

V_{l} (t) = \sum_{i = 1}^{d^{2}} ⟨ ⟨ F_{i} | V_{l} (t) ⟩ ⟩ | F_{i} (⟩ (⟩ .

(79)

As the map

V (t)

is in the form , we can apply to obtain (For simplicity, in this discussion, we omit the explicit time-dependency of the density matrix.)

\begin{matrix} V (t) ρ = \sum_{l} [\sum_{i = 1}^{d^{2}} ⟨ ⟨ F_{i} | V_{l} (t) ⟩ ⟩ F_{i} ρ \sum_{j = 1}^{d^{2}} F_{j}^{†} ⟨ ⟨ V_{l} (t) | F_{j} ⟩ ⟩] = \sum_{i, j = 1}^{d^{2}} c_{i, j} (t) F_{i} ρ F_{j}^{†}, \end{matrix}

(80)

where we have absorbed the summation over the Krauss operators in the terms c_i,j(t) = ∑_l⟨⟨F_i|V_l⟩⟩⟨⟨V_l|F_j⟩⟩. We go back now to the original problem by applying this expansion into the time-derivative of Eq. ,

\begin{matrix} \frac{d ρ}{d t} & = & lim_{Δ t \to 0} \frac{1}{Δ t} (V (Δ t) ρ - ρ) = lim_{Δ t \to 0} (\sum_{i, j = 1}^{d^{2}} c_{i, j} (Δ t) F_{i} ρ F_{j}^{†} - ρ) \\ = & lim_{Δ t \to 0} (\sum_{i, j = 0}^{d^{2} - 1} c_{i, j} (Δ t) F_{i} ρ F_{j}^{†} + \sum_{i = 1}^{d^{2} - 1} c_{i, d^{2}} F_{i} ρ F_{d^{2}}^{†}) \\ (+ \sum_{j = 1}^{d^{2} - 1} c_{d^{2}, j} (Δ t) F_{d^{2}} ρ F_{j}^{†} + c_{d^{2}, d^{2}} (Δ t) F_{d^{2}} ρ F_{d^{2}}^{†} - ρ), \end{matrix}

(81)

where we have separated the summations to take into account that

F_{d^{2}} = \frac{1}{\sqrt{d}} 1_{H}

⁠. By using this property, this equation simplifies to

\begin{matrix} \frac{d ρ}{d t} & = & lim_{Δ t \to 0} \frac{1}{Δ t} () \sum_{i, j = 1}^{d^{2} - 1} c_{i, j} (Δ t) F_{i} ρ F_{j}^{†} + \frac{1}{\sqrt{d}} \sum_{i = 1}^{d^{2} - 1} c_{i, d^{2}} (Δ t) F_{i} ρ \\ + \frac{1}{\sqrt{d}} \sum_{j = 1}^{d^{2} - 1} c_{d^{2}, j} (Δ t) ρ F_{j}^{†} + \frac{1}{d} c_{d^{2}, d^{2}} (Δ t) ρ - ρ)) . \end{matrix}

(82)

The next step is to eliminate the explicit dependence with time. To do so, we define new constants to absorb all the time intervals,

\begin{aligned} g_{i, j} & \equiv lim_{Δ t \to 0} \frac{c_{i, j} (Δ t)}{Δ t} (i, j < d^{2}), \\ g_{i, d^{2}} & \equiv lim_{Δ t \to 0} \frac{c_{i, d^{2}} (Δ t)}{Δ t} (i < d^{2}), \\ g_{d^{2}, j} & \equiv lim_{Δ t \to 0} \frac{c_{d^{2}, j} (Δ t)}{Δ t} (j < d^{2}), \\ g_{d^{2}, d^{2}} & \equiv lim_{Δ t \to 0} \frac{c_{d^{2}, d^{2}} (Δ t) - d}{Δ t} . \end{aligned}

(83)

Introducing these coefficients in Eq. , we obtain an equation with no explicit dependence in time,

\begin{matrix} \frac{d ρ}{d t} = \sum_{i, j = 1}^{d^{2} - 1} g_{i, j} F_{i} ρ F_{j}^{†} + \frac{1}{\sqrt{d}} \sum_{i = 1}^{d^{2} - 1} g_{i, d^{2}} F_{i} ρ + \frac{1}{\sqrt{d}} \sum_{j = 1}^{d^{2} - 1} g_{d^{2}, j} ρ F_{j}^{†} + \frac{g_{d^{2}, d^{2}}}{d} ρ . \end{matrix}

(84)

As we are already summing up over all the Krauss operators, it is useful to define a new operator,

F \equiv \frac{1}{\sqrt{d}} \sum_{i = 1}^{d^{2} - 1} g_{i, d^{2}} F_{i} .

(85)

Applying it to Eq. ,

\frac{d ρ}{d t} = \sum_{i, j = 1}^{d^{2} - 1} g_{i, j} F_{i} ρ F_{j}^{†} + F ρ + ρ F^{†} + \frac{g_{d^{2}, d^{2}}}{d} ρ .

(86)

At this point, we want to separate the dynamics of the density matrix into a Hermitian (equivalent to the von Neumann equation) and an incoherent part. We split the operator F into two to obtain a Hermitian and anti-Hermitian part,

F = \frac{F + F^{†}}{2} + i \frac{F - F^{†}}{2 i} \equiv G - i H,

(87)

where we have used the notation H for the Hermitian part for obvious reasons. If we take this definition to Eq. , we obtain

\frac{d ρ}{d t} = \sum_{i, j = 1}^{d^{2} - 1} g_{i, j} F_{i} ρ F_{j}^{†} + \{G, ρ\} - i [H, ρ] + \frac{g_{d^{2}, d^{2}}}{d} ρ .

(88)

We now define the last operator for this proof,

G_{2} \equiv G + \frac{g_{d^{2}, d^{2}}}{2 d}

⁠, and the expression of the time derivative leads to

\frac{d ρ}{d t} = \sum_{i, j = 1}^{d^{2} - 1} g_{i, j} F_{i} ρ F_{j}^{†} + \{G_{2}, ρ\} - i [H, ρ] .

(89)

Until now, we have imposed the complete positivity of the map, as we have required it to be written in terms of Krauss maps, but we have not used the trace-preserving property. We impose now this property, and by using the cyclic property of the trace, we obtain a new condition,

Tr [\frac{d ρ}{d t}] = Tr [\sum_{i, j = 1}^{d^{2} - 1} F_{j}^{†} F_{i} ρ + 2 G_{2} ρ] = 0 .

(90)

Therefore, G₂ should fulfill

G_{2} = \frac{1}{2} \sum_{i, j = 1}^{d^{2} - 1} g_{i, j} F_{j}^{†} F_{i} ρ .

(91)

By applying this condition, we arrive at the Lindblad master equation,

\frac{d ρ}{d t} = - i [H, ρ] + \sum_{i, j = 1}^{d^{2} - 1} g_{i, j} (F_{i} ρ F_{j}^{†} - \frac{1}{2} \{F_{j}^{†} F_{i}, ρ\}) .

(92)

Finally, by definition of the coefficients, g_i,j can be arranged to form a Hermitian, and therefore diagonalizable, matrix. By diagonalizing it, we obtain the diagonal form of the Lindblad master equation,

\frac{d}{d t} ρ = - i [H, ρ] + \sum_{k} Γ_{k} (L_{k} ρ L_{k}^{†} - \frac{1}{2} \{L_{k} L_{k}^{†}, ρ\}) \equiv L ρ .

(93)

D. Properties of the Lindblad master equation

Some interesting properties of the Lindblad equation are the following:

Under a Lindblad dynamics, if all the jump operators are Hermitian, the purity of a system fulfills $\frac{d}{d t} (Tr [ρ^{2}]) \leq 0$ ⁠. The proof is given in A.
The Lindblad master equation is invariant under unitary transformations of the jump operators,

\sqrt{Γ_{i}} L_{i} \to \sqrt{Γ_{i}^{'}} L_{i}^{'} = \sum_{j} v_{i j} \sqrt{Γ_{j}} L_{j},

(94)

with $v$ representing a unitary matrix. It is also invariant under inhomogeneous transformations in the form

\begin{aligned} L_{i} & \to L_{i}^{'} = L_{i} + a_{i}, \\ H & \to H^{'} = H + \frac{1}{2 i} \sum_{j} Γ_{j} (a_{j}^{*} A_{j} - a_{j} A_{j}^{†}) + b, \end{aligned}

(95)

where $a_{i} \in C$ and $b \in R$ ⁠. The proof of this can be found in Ref. 2 (Sec. III).

Thanks to the previous properties, it is possible to find traceless jump operators without loss of generality.

Example 6: A master equation for a two-level system with decay.

Continuing our example of a two-level atom, we can make it more realistic by including the possibility of atom decay by the emission of a photon. This emission happens due to the interaction of the atom with the surrounding vacuum state. (This is why atoms decay.) The complete quantum system would be in this case the “atom + vacuum” system, and its time evolution should be given by the von Neumann equation (13), where H represents the total “atom + vacuum” Hamiltonian. This system belongs to an infinite-dimension Hilbert space, as the radiation field has infinite modes. If we are interested only in the time dependence of the state of the atom, we can derive a Markovian master equation for the reduced density matrix of the atom (see for instance, Refs. 1 and 2). The master equation we will study is

\frac{d}{d t} ρ (t) = - i [H, ρ] + Γ (σ^{-} ρ σ^{+} - \frac{1}{2} \{σ^{+} σ^{-}, ρ\}),

(96)

where Γ is the coupling between the atom and the vacuum.

In the Fock-Liouvillian space [following the same ordering as in Eq. (25)], the Liouvillian corresponding to evolution (96) is

L = (\begin{matrix} 0 & i Ω & - i Ω & Γ \\ i Ω & - i E - \frac{Γ}{2} & 0 & - i Ω \\ - i Ω & 0 & - i E - \frac{Γ}{2} & i Ω \\ 0 & - i Ω & i Ω & - Γ \end{matrix}) .

(97)

Expressing explicitly the set of differential equations, we obtain

\begin{aligned} {\dot{ρ}}_{00} & = i Ω ρ_{01} - i Ω ρ_{10} + Γ ρ_{11}, \\ {\dot{ρ}}_{01} & = i Ω ρ_{00} - (i E - \frac{Γ}{2}) ρ_{01} - i Ω ρ_{11}, \\ {\dot{ρ}}_{10} & = - i Ω ρ_{00} (- i E - \frac{Γ}{2}) ρ_{10} + i Ω ρ_{11}, \\ {\dot{ρ}}_{10} & = - i Ω ρ_{01} + i Ω ρ_{10} - Γ ρ_{11} . \end{aligned}

(98)

End of example.

VI. RESOLUTION OF THE LINDBLAD MASTER EQUATION

A. Integration

To calculate the time evolution of a system determined by a master equation in the form (96), we need to solve a set of equations with as many equations as the dimension of the density matrix. In our example, this means to solve a 4 variable set of equations, but the dimension of the problem increases exponentially with the system size. Because of this, for bigger systems, techniques for dimension reduction are required.

To solve systems of partial differential equations, there are several canonical algorithms. This can be done analytically only for a few simple systems and by using sophisticated techniques as damping bases.²⁹ In most cases, we have to rely on numerical approximated methods. One of the most popular approaches is the 4th-order Runge-Kutta algorithm (see, for instance, Ref. 30, for an explanation of the algorithm). By integrating the equations of motion, we can calculate the density matrix at any time t.

The steady-state of a system can be obtained by evolving it for a long time $(t \to \infty)$ ⁠. Unfortunately, this method presents two difficulties. First, if the dimension of the system is big, the computing time would be huge. This means that for systems beyond a few qubits, it will take too long to reach the steady-state. Even worse is the problem of stability of the algorithms for integrating differential equations. Due to small errors in the calculation of derivatives by the use of finite differences, the trace of the density matrix may not be constantly equal to one. This error accumulates during the propagation of the state, giving nonphysical results after a finite time. One solution to this problem is the use of algorithms specifically designed to preserve the trace, as the Crank-Nicholson algorithm.³¹ The problem with this kind of algorithm is that they consume more computational power than Runge-Kutta, and therefore they are not useful to calculate the long-time behavior of big systems. An analysis of different methods and their advantages and disadvantages can be found in Ref. 32.Example 7: Time dependency of the two-level system with decay.

In this section, we show some results of solving Eq. (96) and calculating the density matrix as a function of time. A Mathematica notebook solving this problem can be found in Ref. 20. To illustrate the time behavior of this system, we calculate the evolution for different state parameters. In all cases, we start with an initial state that represents the state being excited ρ₁₁ = 1, with no coherence between different states, meaning ρ₀₁ = ρ₁₀ = 0. If the decay parameter Γ is equal to zero, the problem reduces to solve the von Neumann equation, and the result is displayed in Fig. 2. The other extreme case would be a system with no coherent dynamics (Ω = 0) but with decay. In this case, as shown in Fig. 4, we observe an exponential decay of the population of the excited state.

FIG. 4.

View large Download slide

Population dynamics under a pure incoherent dynamics (Γ = 0.1, n = 1, Ω = 0, and E = 1). The blue line represents ρ₁₁ and the orange one ρ₀₀.

Finally, we can calculate the dynamics of a system with both coherent driving and decay. In this case, both behaviors coexist and there are oscillations and decay. This is displayed in Fig. 5.End of example.

FIG. 5.

View large Download slide

Population dynamics under a pure incoherent dynamics (Γ = 0.2, n = 1, Ω = 0, and E = 1). The blue line represents ρ₁₁ and the orange one ρ₀₀.

B. Diagonalization

As we have discussed before, in the Fock-Liouville space, the Liouvillian corresponds to a matrix (in general complex, non-Hermitian, and nonsymmetric). By diagonalizing it, we can calculate both the time-dependent and the steady-state of the density matrices. For most purposes, in the short time regime, integrating the differential equations may be more efficient than diagonalizing. This is due to the high dimensionality of the Liouvillian that makes the diagonalization process very costly in computing power. On the other hand, in order to calculate the steady-state, the diagonalization is the mostly used method due to the problems of integrating the equation of motions discussed in Sec. VI A.

Let us see first how we use diagonalization to calculate the time evolution of a system. As the Liouvillian matrix is non-Hermitian, we cannot apply the spectral theorem to it and it may have different left and right eigenvectors. For a specific eigenvalue Λ_i, we can obtain the eigenvectors $| Λ_{i}^{R} (⟩ (⟩$ and $| Λ_{i}^{L} (⟩ (⟩$ s.t.

\begin{aligned} \tilde{L} | Λ_{i}^{R} (⟩ (⟩ = Λ_{i} | Λ_{i}^{R} (⟩ (⟩, \\ ⟨) ⟨) Λ_{i}^{L} | \tilde{L} = Λ_{i} ⟨) ⟨) Λ_{i}^{L} | . \end{aligned}

(99)

An arbitrary system can be expanded in the eigenbasis of $\tilde{L}$ as^1,33

| ρ (0) (⟩ (⟩ = \sum_{i} | Λ_{i}^{R} (⟩ (⟩ ⟨ ⟨ Λ_{i}^{L} | ρ (0) ⟩ ⟩ .

(100)

Therefore, the state of the system at a time t can be calculated in the form

| ρ (t) (⟩ (⟩ = \sum_{i} e^{Λ_{i} t} | Λ_{i}^{R} (⟩ (⟩ ⟨ ⟨ Λ_{i}^{L} | ρ (0) ⟩ ⟩ .

(101)

Note that in this case to calculate the state a time t we do not need to integrate into the interval $[0, t]$ ⁠, as we have to do if we use a numerical solution of the differential set of equations. This is an advantage when we want to calculate long-time behavior. Furthermore, to calculate the steady-state of a system, we can look to the eigenvector that has zero eigenvalue, as this is the only one that survives when t → ∞.

For any finite system, Evans’ theorem ensures the existence of at least one zero eigenvalue of the Liouvillian matrix.^34,35 The eigenvector corresponding to this zero eigenvalue would be the steady-state of the system. In exceptional cases, a Liouvillian can present more than one zero eigenvalues due to the presence of symmetry in the system.^26,27,36 This is a nongeneric case, and for most purposes, we can assume the existence of a unique fixed point in the dynamics of the system. Therefore, diagonalizing can be used to calculate the steady-state without calculating the full evolution of the system. This can be done even analytically for small systems, and when numerical approaches are required, this technique gives better precision than integrating the equations of motion. The spectrum of Liouvillian superoperators has been analyzed in several recent papers.^33,37Example 8: Spectrum-analysis of the Liouvillian for the two-level system with decay.

Here, we diagonalize (97) and obtain its steady state. A Mathematica notebook solving this problem can be downloaded from Ref. 20. This specific case is straightforward to diagonalize as the dimension of the system is very low. We obtain 4 different eigenvalues, two of them are real while the other two form a conjugated pair. Figure 6 displays the spectrum of the superoperator $L$ given in (97).

FIG. 6.

View large Download slide

Spectrum of the Liouvillian matrix given by (97) for the general case of both coherent and incoherent dynamics (Γ = 0.2, n = 1, Ω = 0, E = 1).

As there is only one zero eigenvalue, we can conclude that there is only one steady-state, and any initial density matrix will evolve to it after an infinite-time evolution. By selecting the right eigenvector corresponding to the zero-eigenvalue and normalizing it we obtain the density matrix. This can be done even analytically. The result is the matrix,

\begin{matrix} ρ_{S S} = (\begin{matrix} \frac{(1 + n) (4 E^{2} + {(Γ + 2 n Γ)}^{2}) + 4 (1 + 2 n) Ω^{2}}{(1 + 2 n) (4 E^{2} + {(Γ + 2 n Γ)}^{2} + 8 Ω^{2})} & \frac{2 (- 2 E - i (Γ + 2 n Γ)) Ω}{(1 + 2 n) (4 E^{2} + {(Γ + 2 n Γ)}^{2} + 8 Ω^{2})} \\ \frac{2 (- 2 E + i (Γ + 2 n Γ)) Ω}{(1 + 2 n) (4 E^{2} + {(Γ + 2 n Γ)}^{2} + 8 Ω^{2})} & \frac{n (4 E^{2} + {(Γ + 2 n Γ)}^{2}) + 4 (1 + 2 n) Ω^{2}}{(1 + 2 n) (4 E^{2} + {(Γ + 2 n Γ)}^{2} + 8 Ω^{2})} \end{matrix}) . \end{matrix}

(102)

End of example.

ACKNOWLEDGMENTS

The author wants to acknowledge the Spanish Ministry and the Agencia Española de Investigación (AEI) for financial support under Grant No. FIS2017-84256-P (FEDER funds).

APPENDIX: PROOF OF $\frac{d}{d t} T r [ρ^{2}] \leq 0$

In this Appendix, we proof that under the Lindblad dynamics given by Eq. (93), the purity of a density matrix fulfills that $\frac{d}{d t} Tr [ρ^{2}] \leq 0$ if all the jump operators of the Lindblad dynamics are Hermitian.

We start just by interchanging the trace and the derivative. As the trace is a linear operation, it commutes with the derivation, and we have

\frac{d}{d t} (Tr [ρ^{2}]) = Tr [\frac{d ρ^{2}}{d t}] = Tr [2 ρ \dot{ρ}],

(A1)

where we have used the cyclic property of the trace operator. (This property is used along all the demonstration without explicitly mentioning it.) By inserting the Lindblad Eq. (93) into the rhs of (A1), we obtain

\begin{matrix} \frac{d}{d t} (Tr [ρ^{2}]) & = & - \frac{i}{ℏ} Tr [(2 ρ (H ρ - ρ H))] + 2 \sum_{k} Γ_{k} Tr [ρ L_{k} ρ L_{k}^{†}] \\ - 2 \sum_{k} Γ_{k} Tr [ρ^{2} L_{k}^{†} L_{k}] . \end{matrix}

(A2)

The first term is zero. Therefore, the inequality we want to prove becomes equivalent to

\sum_{k} Γ_{k} Tr [ρ L_{k} ρ L_{k}^{†}] \leq \sum_{k} Γ_{k} Tr [ρ^{2} L_{k}^{†} L_{k}] .

(A3)

As the density matrix is Hermitian, we can diagonalize it to obtain its eigenvalues (⁠ $Λ_{i} \in R$ ⁠) and its corresponding eigenvectors (|Λ_i⟩). The density matrix is diagonal in its own eigenbasis and can be expressed as (this eigenbasis changes with time, of course, but the proof is valid as the inequality should be fulfilled at any time)

ρ \to \tilde{ρ} = \sum_{i} Λ_{i} | Λ_{i} (⟩ ⟨) Λ_{i} |,

(A4)

where we assume an ordering of the eigenvalues in the form Λ₀ ≥ Λ₁ ≥ ⋯ ≥Λ_d.

We rename the jump operators in this basis as ${\tilde{L}}_{i}$ a. Expanding each term of the inequality (A3) in this basis, we obtain

\begin{matrix} \sum_{k} Γ_{k} Tr [ρ L_{k} ρ L_{k}^{†}] & = & \sum_{k} Γ_{k} Tr [(\sum_{i} Λ_{i} | Λ_{i} (⟩ ⟨) Λ_{i} |) {\tilde{L}}_{k} (\sum_{j} Λ_{j} | Λ_{j} (⟩ ⟨) Λ_{j} |) {\tilde{L}}_{k}] \\ = & \sum_{k} Γ_{k} \sum_{i, j} Λ_{i} Λ_{j} Tr [{\tilde{L}}_{k}^{†} | Λ_{i} (⟩ ⟨) Λ_{i} | {\tilde{L}}_{k} | Λ_{j} (⟩ ⟨) Λ_{j} |] \\ = & \sum_{k} Γ_{k} \sum_{i, j} Λ_{i} Λ_{j} Tr [{|⟨) Λ_{i} | {\tilde{L}}_{k} | Λ_{j} (⟩|}^{2}] \\ = & \sum_{k} Γ_{k} \sum_{i, j} Λ_{i} Λ_{j} x_{i j}^{(k)}, \end{matrix}

(A5)

where we have introduced the coefficients $x_{i j}^{(k)} \equiv {|⟨) Λ_{i} | {\tilde{L}}_{k} | Λ_{j} (⟩|}^{2}$ ⁠. As the operators L_k are Hermitian, these coefficients fulfill $x_{i j}^{(k)} = x_{j i}^{(k)}$ ⁠.

The second term is expanded as

\begin{matrix} \sum_{k} Γ_{k} Tr [ρ^{2} L_{k}^{†} L_{k}] & = & \sum_{k} Γ_{k} Tr [(\sum_{i} Λ_{i} | Λ_{i} (⟩ ⟨) Λ_{i} |) (\sum_{j} Λ_{j} | Λ_{j} (⟩ ⟨) Λ_{j} |) {\tilde{L}}_{k}^{†} {\tilde{L}}_{k}] \\ = & \sum_{k} Γ_{k} \sum_{i j} Λ_{i} Λ_{j} Tr [{\tilde{L}}_{k} | Λ_{i} (⟩ ⟨) Λ_{j} | {\tilde{L}}_{k}^{†} ⟨ Λ_{i} | Λ_{j} ⟩] \\ = & \sum_{k} Γ_{k} \sum_{i} Λ_{i}^{2} Tr [{\tilde{L}}_{k} | Λ_{i} (⟩ ⟨) Λ_{i} | {\tilde{L}}_{k}^{†}] \\ = & \sum_{k} Γ_{k} \sum_{i} Λ_{i}^{2} Tr [{\tilde{L}}_{k} | Λ_{i} (⟩ ⟨) Λ_{i} | {\tilde{L}}_{k}^{†} (\sum_{j} | Λ_{j} (⟩ ⟨) Λ_{j} |)] \\ = & \sum_{k} Γ_{k} \sum_{i j} Λ_{i}^{2} Tr [⟨) Λ_{j} | {\tilde{L}}_{k} | Λ_{i} (⟩ + ⟨) Λ_{i} | {\tilde{L}}_{k} | Λ_{j} (⟩] \\ = & \sum_{k} Γ_{k} \sum_{i j} Λ_{i}^{2} x_{i j}, \end{matrix}

(A6)

where we have used the closure relation in the density matrix eigenbasis, $1_{H} = \sum_{j} | Λ_{j} (⟩ ⟨) Λ_{j} |$ ⁠. The inequality can now be written as

\sum_{k} Γ_{k} \sum_{i j} Λ_{i} Λ_{j} x_{i j} \leq \sum_{k} Γ_{k} \sum_{i j} Λ_{i}^{2} x_{i j} .

(A7)

As x_ij = x_ji, we can reorder the ij sum in the following way:

\begin{matrix} \sum_{k} Γ_{k} \sum_{i} (\sum_{j \leq i} 2 Λ_{i} Λ_{j} x_{i j}^{(k)} + Λ_{i}^{2} x_{i i}^{(k)}) \\ \leq \sum_{k} Γ_{k} \sum_{i} (\sum_{j < i} (Λ_{i}^{2} + Λ_{j}^{2}) x_{i j}^{(k)} + Λ_{i}^{2} x_{i i}^{(k)}) . \end{matrix}

(A8)

Therefore, we can reduce the proof of this inequality to the proof of a set of inequalities,

2 Λ_{i} Λ_{j} x_{i j}^{(k)} \leq (Λ_{i}^{2} + Λ_{j}^{2}) x_{i j}^{(k)}, \forall (k, i, j) .

(A9)

It is obvious that (A9) ⇒ (A8) (but not the other way around). The inequalities (A9) are easily proved just by taking into account that $x_{i j}^{(k)} \geq 0$ and applying the triangular inequality.

REFERENCES

1.

C. W.

Gardiner

and

P.

Zoller

,

Quantum Noise

(

Springer

,

Berlin

,

2000

).

Google Scholar

Crossref

2.

H. P.

Breuer

and

F.

Petruccione

,

The Theory of Open Quantum Systems

(

Oxford University Press

,

2002

).

Google Scholar

3.

A.

Rivas

and

S.

Huelga

,

Open Quantum Systems. An Introduction

(

Springer

,

New York

,

2012

).

Google Scholar

Crossref

4.

G.

Lindblad

, “

On the generators of quantum dynamical semigroups

,”

Commun. Math. Phys.

119

,

48

(

1976

).

https://doi.org/10.1007/bf01608499

Google Scholar

Crossref

5.

V.

Gorini

,

A.

Kossakowski

, and

E. C.

Sudarsahan

, “

Completely positive semigroups of n-level systems

,”

J. Math. Phys.

17

,

821

(

1976

).

https://doi.org/10.1063/1.522979

Google Scholar

Crossref

6.

D.

Manzano

and

E.

Kyoseva

, “

An atomic symmetry-controlled thermal switch

,”

Sci. Rep.

6

,

31161

(

2016

).

https://doi.org/10.1038/srep31161

Google Scholar

Crossref

PubMed

7.

T.

Prosen

, “

Open XXZ spin chain: Nonequilibrium steady state and a strict bound on ballistic transport

,”

Phys. Rev. Lett.

106

,

217206

(

2011

).

https://doi.org/10.1103/physrevlett.106.217206

Google Scholar

Crossref

PubMed

8.

D.

Manzano

,

M.

Tiersch

,

A.

Asadian

, and

H. J.

Briegel

, “

Quantum transport efficiency and Fourier’s law

,”

Phys. Rev. E

86

,

061118

(

2012

).

https://doi.org/10.1103/physreve.86.061118

Google Scholar

Crossref

9.

D.

Manzano

,

C.

Chuang

, and

J.

Cao

, “

Quantum transport in d-dimensional lattices

,”

New J. Phys.

18

,

043044

(

2015

).

https://doi.org/10.1088/1367-2630/18/4/043044

Google Scholar

Crossref

10.

B.

Olmos

,

I.

Lesanovsky

, and

J. P.

Garrahan

, “

Facilitated spin models of dissipative quantum glasses

,”

Phys. Rev. Lett.

109

,

020403

(

2012

).

https://doi.org/10.1103/physrevlett.109.020403

Google Scholar

Crossref

PubMed

11.

J.

Metz

,

M.

Trupke

, and

A.

Beige

, “

Robust entanglement through macroscopic quantum jumps

,”

Phys. Rev. Lett.

97

,

040503

(

2006

).

https://doi.org/10.1103/physrevlett.97.040503

Google Scholar

Crossref

PubMed

12.

R.

Jones

,

J. A.

Needham

,

I.

Lesanovsky

, and

F.

Intravaia

, “

Beatriz Olmos modified dipole-dipole interaction and dissipation in an atomic ensemble near surfaces

,”

Phys. Rev. A

97

,

053841

(

2018

).

https://doi.org/10.1103/physreva.97.053841

Google Scholar

Crossref

13.

D. A.

Lidar

,

I. L.

Chuang

, and

K. B.

Whaley

, “

Decoherence-free subspaces for quantum computation

,”

Phys. Rev. Lett.

81

(

12

),

2594

(

1998

).

https://doi.org/10.1103/physrevlett.81.2594

Google Scholar

Crossref

14.

B.

Kraus

,

H. P.

Büchler

,

S.

Diehl

,

A.

Kantian

,

A.

Micheli

, and

P.

Zoller

, “

Preparation of entangled states by quantum Markov processes

,”

Phys. Rev. A

78

,

042307

(

2008

).

https://doi.org/10.1103/physreva.78.042307

Google Scholar

Crossref

15.

T. A.

Brun

, “

Continuous measurements, quantum trajectories, and decoherent histories

,”

Phys. Rev. A

61

,

042107

(

2000

).

https://doi.org/10.1103/physreva.61.042107

Google Scholar

Crossref

16.

M.

Schlosshauer

,

Decoherence and the Quantum-to-Classical Transition

(

Springer

,

New York

,

2007

).

Google Scholar

17.

M.

Plenio

and

S.

Huelga

, “

Dephasing-assisted transport: Quantum networks and biomolecules

,”

New J. Phys.

10

,

113019

(

2008

).

https://doi.org/10.1088/1367-2630/10/11/113019

Google Scholar

Crossref

18.

M.

Mohseni

,

P.

Rebentrost

,

S.

Lloyd

, and

A.

Aspuru-Guzik

, “

Enviroment-assisted quantum walks in photosynthetic energy transfer

,”

J. Chem. Phys.

129

,

174106

(

2008

).

https://doi.org/10.1063/1.3002335

Google Scholar

Crossref

PubMed

19.

D.

Manzano

, “

Quantum transport in quantum networks and photosynthetic complexes at the steady state

,”

PLoS One

8

(

2

),

e57041

(

2013

).

https://doi.org/10.1371/annotation/b73c06e9-1e54-497c-bafb-f56de94f2f18

Google Scholar

Crossref

PubMed

20.

See https://ic1.ugr.es/manzano/Descargas/Lindblad/Lindblad_Manzano.zip for downloading the notebooks used to generate all the data of this paper.

21.

L.

Debnath

and

P.

Mikusińki

,

Introduction to Hilbert Spaces with Applications

(

Elsevier Academic Press

,

2005

).

Google Scholar

22.

M. A.

Nielsen

and

I. L.

Chuang

,

Quantum Computation and Quantum Information

(

Cambridge University Press

,

Cambridge

,

2000

).

Google Scholar

23.

J. J.

Sakurai

,

Modern Quantum Mechanics

(

Addison-Wesley Publishing Co.

,

1994

).

Google Scholar

24.

A.

Galindo

and

P.

Pascual

,

Quantum Mechanics I

(

Springer

,

Berlin

,

1990

).

Google Scholar

Crossref

25.

A.

Peres

,

Quantum Theory: Concepts and Methods

(

Kluwer Academic Publishers

,

1995

).

Google Scholar

26.

D.

Manzano

and

P. I.

Hurtado

, “

Harnessing symmetry to control quantum transport

,”

Adv. Phys.

67

(

1

),

1

(

2018

).

https://doi.org/10.1080/00018732.2018.1519981

Google Scholar

Crossref

27.

D.

Manzano

and

P. I.

Hurtado

, “

Symmetry and the thermodynamics of currents in open quantum systems

,”

Phys. Rev. B

90

,

125138

(

2014

).

https://doi.org/10.1103/physrevb.90.125138

Google Scholar

Crossref

28.

M. M.

Wilde

,

Quantum Information Theory

(

Cambridge University Press

,

Cambridge

,

2017

).

Google Scholar

29.

H. J.

Briegel

and

B. G.

Englert

, “

Quantum optical master equation: The use of damping bases

,”

Phys. Rev. A

47

,

3311

(

1993

).

https://doi.org/10.1103/physreva.47.3311

Google Scholar

Crossref

PubMed

30.

W. H.

Press

,

S. A.

Teukolsky

,

W. T.

Vetterling

, and

B. P.

Flannery

,

Numerical Recipes

(

Cambridge University Press

,

Cambridge

,

2007

).

Google Scholar

31.

A.

Goldberg

,

H.

Schey

, and

J. L.

Schwartz

, “

Computer-generated motion pictures of one-dimensional quantum-mechanical transmission and reflection phenomena

,”

Am. J. Phys.

35

,

177

(

1967

).

https://doi.org/10.1119/1.1973991

Google Scholar

Crossref

32.

M.

Riesch

and

C.

Jirauschek

, “

Analyzing the positivity preservation of numerical methods for the Liouville-von Neumann equation

,”

J. Comput. Phys.

390

,

290

(

2019

).

https://doi.org/10.1016/j.jcp.2019.04.006

Google Scholar

Crossref

33.

J.

Thingna

,

D.

Manzano

, and

J.

Cao

, “

Dynamical signatures of molecular symmetries in nonequilibrium quantum transport

,”

Sci. Rep.

6

,

28027

(

2016

).

https://doi.org/10.1038/srep28027

Google Scholar

Crossref

PubMed

34.

D. E.

Evans

, “

Irreducible quantum dynamical semigroups

,”

Commun. Math. Phys.

54

,

293

(

1977

).

https://doi.org/10.1007/bf01614091

Google Scholar

Crossref

35.

D. E.

Evans

and

H.

Hance-Olsen

, “

The generators of positive semigroups

,”

J. Funct. Anal.

32

,

207

(

1979

).

https://doi.org/10.1016/0022-1236(79)90054-5

Google Scholar

Crossref

36.

B.

Buča

and

T.

Prosen

, “

A note on symmetry reductions of the Lindblad equation: Transport in constrained open spin chains

,”

New J. Phys.

14

,

073007

(

2012

).

https://doi.org/10.1088/1367-2630/14/7/073007

Google Scholar

Crossref

37.

V. V.

Albert

and

L.

Jiang

, “

Symmetries and conserved quantities in Lindblad master equations

,”

Phys. Rev. A

89

,

022118

(

2014

).

https://doi.org/10.1103/physreva.89.022118

Google Scholar

Crossref

38.

For this discussion, we use a bipartite system as an example. The extension to a general multipartite system is straightforward.

39.

This is clear for finite systems because in this case scalar space and Hilbert space are the same things. It also happens to be true for infinite spaces.

40.

From now on, we will not write the identity operators of the Hamiltonian parts explicitly when they can be inferred from the context.

41.

A. G.

Redfield

,

Advances in Magnetic and Optical Resonance

(

Elsevier

,

1965

),

1

, Vol. 1, p.

1

.

Google Scholar

Crossref

2020

Author(s)

A short introduction to the Lindblad master equation

I. INTRODUCTION

II. MATHEMATICAL BASIS

Glossary:

III. (VERY SHORT) INTRODUCTION TO QUANTUM MECHANICS

IV. THE FOCK-LIOUVILLE HILBERT SPACE. THE LIOUVILLE SUPEROPERATOR

V. CPT-MAPS AND THE LINDBLAD MASTER EQUATION

A. Completely positive maps

B. Derivation of the Lindblad equation from microscopic dynamics

C. Derivation of the Lindblad equation as a CPT generator

D. Properties of the Lindblad master equation

VI. RESOLUTION OF THE LINDBLAD MASTER EQUATION

A. Integration

B. Diagonalization

ACKNOWLEDGMENTS

APPENDIX: PROOF OF $\frac{d}{d t} T r [ρ^{2}] \leq 0$

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

A short introduction to the Lindblad master equation

I. INTRODUCTION

II. MATHEMATICAL BASIS

Glossary:

III. (VERY SHORT) INTRODUCTION TO QUANTUM MECHANICS

IV. THE FOCK-LIOUVILLE HILBERT SPACE. THE LIOUVILLE SUPEROPERATOR

V. CPT-MAPS AND THE LINDBLAD MASTER EQUATION

A. Completely positive maps

B. Derivation of the Lindblad equation from microscopic dynamics

C. Derivation of the Lindblad equation as a CPT generator

D. Properties of the Lindblad master equation

VI. RESOLUTION OF THE LINDBLAD MASTER EQUATION

A. Integration

B. Diagonalization

ACKNOWLEDGMENTS

APPENDIX: PROOF OF ddtTrρ2≤0

REFERENCES

Citing articles via

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

APPENDIX: PROOF OF $\frac{d}{d t} T r [ρ^{2}] \leq 0$