A short introduction to the Lindblad Master Equation

TL;DR

This paper provides a clear, self-contained introduction to the Lindblad Master Equation, essential for understanding Markovian dynamics in open quantum systems, with derivation and solution methods explained for newcomers.

Contribution

It offers a straightforward, accessible presentation of the Lindblad Master Equation, including its derivation and solution techniques, aimed at readers new to the field.

Findings

Provides a self-contained derivation of the Lindblad Equation

Explains methods for solving the Lindblad Master Equation

Clarifies the role of the Lindblad Equation in quantum dynamics

Abstract

The theory of open quantum system 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.

Equations

$$\rho \equiv \sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|.$$

$$\rho =\left(\begin{array}{cccc}{\rho }_{00}& {\rho }_{01}& \cdots & {\rho }_{0N}\\ {\rho }_{10}& {\rho }_{11}& \cdots & {\rho }_{1N}\\ ⋮& ⋮& \ddots & ⋮\\ {\rho }_{N0}& {\rho }_{N1}& \cdots & {\rho }_{NN}\end{array}\right),$$

$$\rho =\left(\begin{array}{cc}{\rho }_{00}& {\rho }_{01}\\ {\rho }_{10}& {\rho }_{11}\end{array}\right)={\rho }_{00}|0⟩⟨0|+{\rho }_{01}|0⟩⟨1|+{\rho }_{10}|1⟩⟨0|+{\rho }_{11}|1⟩⟨1|,$$

$$O=\sum _i{a}_i|a_i⟩⟨a_i|,$$

$$P(a_i)={\left|⟨\psi |a_i⟩\right|}^2.$$

$$⟨O⟩=⟨\psi |O|\psi ⟩.$$

$$P(a_i)=\text{Tr}\left[|a_i⟩⟨a_i|\rho \right],$$

$$⟨O⟩=\text{Tr}\left[O\rho \right].$$

$$H=E_0|0⟩⟨0|+E_1|1⟩⟨1|.$$

$$\frac{d}{dt}|\psi (t)⟩=-i\hslash H|\psi (t)⟩,$$

$$|\psi (t)⟩=e^{-iHt}|\psi (0)⟩.$$

$$|\psi (t)⟩=U|\psi (0)⟩,$$

$$\dot{\rho }=-i\left[H,\rho \right]\equiv \mathcal{L}\rho ,$$

$$\frac{d}{dt}\text{Tr}\left[{\rho }^2\right]=\text{Tr}\left[\frac{d{\rho }^2}{dt}\right]=\text{Tr}\left[2\rho \dot{\rho }\right]=-2i\text{Tr}\left[\rho \left(H\rho -\rho H\right)\right]=0,$$

$$H_{\text{free}}=E_0|0⟩⟨0|+E_1|1⟩⟨1|,$$

$$H_{\text{free}}=E|1⟩⟨1|,$$

$$H=E|1⟩⟨1|+\Omega \left(|0⟩⟨1|+|1⟩⟨0|\right),$$

$${\rho }_a\equiv {\text{Tr}}_b\left[\rho \right],$$

$${\text{Tr}}_b\left[\sum _{i,j,k,l}|a_i⟩⟨a_j|\otimes |b_k⟩⟨b_l|\right]\equiv \sum _{i,j}|a_i⟩⟨a_j|\text{Tr}\left[\sum _{k,l}|b_k⟩⟨b_l|\right].$$

$$|\psi {⟩}_G=|00⟩.$$

$$|\psi {⟩}_S$$

$$|\psi {⟩}_E=\frac{1}{\sqrt{2}}\left(|00⟩+|11⟩\right).$$

$${\rho }_E=|\psi ⟩⟨\psi {|}_E=\frac{1}{2}\left(|00⟩⟨00|+|00⟩⟨11|+|11⟩⟨00|+|11⟩⟨11|\right).$$

$${\rho }_E^{(1)}=⟨0{|}_2{\rho }_E|0{⟩}_2+⟨1{|}_2{\rho }_E|1{⟩}_2=\frac{1}{2}\left(|00⟩⟨00{|}_1+|11⟩⟨11{|}_2\right).$$

$$|\rho ⟩⟩=\left(\begin{array}{c}{\rho }_{00}\\ {\rho }_{01}\\ {\rho }_{10}\\ {\rho }_{11}\end{array}\right).$$

$$\tilde{{\cal L}}=\left(\begin{array}[]{cccc}0&i\Omega&-i\Omega&0\\ i\Omega&iE&0&-i\Omega\\ -i\Omega&0&-iE&i\Omega\\ 0&-i\Omega&i\Omega&0\end{array}\right),$$

$$\begin{pmatrix}\dot{\rho}_{00}\\ \dot{\rho}_{01}\\ \dot{\rho}_{10}\\ \dot{\rho}_{11}\end{pmatrix}=\left(\begin{array}[]{cccc}0&i\Omega&-i\Omega&0\\ i\Omega&iE&0&-i\Omega\\ -i\Omega&0&-iE&i\Omega\\ 0&-i\Omega&i\Omega&0\end{array}\right)\begin{pmatrix}\rho_{00}\\ \rho_{01}\\ \rho_{10}\\ \rho_{11}\end{pmatrix}$$

$${\rho }_B=\frac{1}{2}\left(|0⟩⟨0|\otimes |1⟩⟨1|+|1⟩⟨1|\otimes |0⟩⟨0|+|0⟩⟨1|\otimes |1⟩⟨0|+|1⟩⟨0|\otimes |0⟩⟨1|\right),$$

$$\displaystyle\rho_{B}=\frac{1}{2}\left\{\left(\begin{array}[]{cc}1&0\\ 0&0\\ \end{array}\right)\otimes\left(\begin{array}[]{cc}0&0\\ 0&1\\ \end{array}\right)+\left(\begin{array}[]{cc}0&0\\ 0&1\\ \end{array}\right)\otimes\left(\begin{array}[]{cc}1&0\\ 0&0\\ \end{array}\right)\right.$$

$$\displaystyle\left.\otimes\left(\begin{array}[]{cc}0&0\\ 1&0\\ \end{array}\right)\otimes\left(\begin{array}[]{cc}0&1\\ 0&0\\ \end{array}\right)+\left(\begin{array}[]{cc}0&1\\ 0&0\\ \end{array}\right)\otimes\left(\begin{array}[]{cc}0&0\\ 1&0\\ \end{array}\right)\right\}.$$

$$\rho_{B}=\left(\begin{array}[]{cccc}0&0&0&0\\ 0&1&1&0\\ 0&1&1&0\\ 0&0&0&0\end{array}\right),$$

Full text

I Introduction

Open quantum system techniques are vital for many studies in quantum mechanics gardiner_00 ; breuer_02 ; rivas_12 . This happens because closed quantum systems are just an idealisation of real systems111The 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 Figure 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 “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 systems 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.

One particular, and interesting, case of study is the dynamics of a system connected to several baths modelled 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) lindblad:cmp76 ; gorini:jmp76 . It is difficult to overemphasize the importance of this Master Equation. It plays an important role in fields as quantum optics gardiner_00 ; manzano:sr16 , condensed matter prosen:prl11 ; manzano:pre12 ; manzano:njp16 ; olmos:prl12 , atomic physics metz:prl06 ; jones:pra18 , quantum information lidar:prl98 ; kraus:08 , decoherence brun:pra00 ; schlosshauer_07 , and quantum biology plenio:njp08 ; mohseni:jcp08 ; manzano:po13 .

The purpose of this paper is to provide basic knowledge about the Lindblad Master Equation. In Section II, the mathematical requirements are introduced while in Section 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, that is especially useful to work in this problem. In Section V, we define the concept of CPT-Maps, derive the Lindblad Master Equation from two different approaches, and we discus several properties of the equation. Finally, Section VI is devoted to the resolution of the master equation using different methods. To deepen in the techniques of solving the Lindblad equation, an example consisting of a two-level system with decay is analysed, illustrating the content of every section. The problems proposed are solved by the use of Mathematica notebooks that can be found at notebook .

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 in the field of Hilbert spaces we recommend the book by Debnath and Mikusińki debnath_05 . 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 nielsen_00 . It is also required some basic knowledge about infinitesimal calculus, like integration, derivation, and the resolution of simple differential equations, To help the readers, we have made a glossary of the most used mathematical terms. It can be used also as a checklist of terms the reader should be familiar with.

Glossary:

• •

${\cal H}$ represents a Hilbert space, usually the space of pure states of a system.

• •

$|{\psi}\rangle\in{\cal H}$ represents a vector of the Hilbert space ${\cal H}$ (a column vector).

• •

$\langle{\psi}|\in{\cal H}$ represents a vector of the dual Hilbert space of ${\cal H}$ (a row vector).

• •

$\langle{\psi}|{\phi}\rangle\in\mathbb{C}$ is the scalar product of vectors $|{\psi}\rangle$ and $|{\phi}\rangle$.

• •

$\left|\left||{\psi}\rangle\right|\right|$ is the norm of vector $|{\psi}\rangle$. $\left|\left||{\psi}\rangle\right|\right|\equiv\sqrt{\langle{\psi}|{\psi}\rangle}$.

• •

$B({\cal H})$ represents the space of bounded operators acting on the Hilbert space $B:{\cal H}\to{\cal H}$.

• •

$\mathbb{1}_{{\cal H}}\in B({\cal H})$ is the Identity Operator of the Hilbert space ${\cal H}$ s.t. $\mathbb{1}_{{\cal H}}|{\psi}\rangle=|{\psi}\rangle,\;\;\forall|{\psi}\rangle\in{\cal H}$.

• •

$|{\psi}\rangle\!\langle{\phi}|\in B({\cal H})$ is the operator such that $\left(|{\psi}\rangle\!\langle{\phi}|\right)|{\varphi}\rangle=\langle{\phi}|{\varphi}\rangle|{\psi}\rangle,\;\;\forall|{\varphi}\rangle\in{\cal H}$.

• •

$O^{\dagger}\in B({\cal H})$ is the Hermitian conjugate of the operator $O\in B({\cal H})$.

• •

$U\in B({\cal H})$ is a unitary operator iff $UU^{\dagger}=U^{\dagger}U=\mathbb{1}$.

• •

$H\in B({\cal H})$ is a Hermitian operator iff $H=H^{\dagger}$.

• •

$A\in B({\cal H})$ is a positive operator $\left(A>0\right)$ iff $\langle{\phi}|A|{\phi}\rangle\geq 0,\;\;\forall|{\phi}\rangle\in{\cal H}$

• •

$P\in B({\cal H})$ is a proyector iff $PP=P$.

• •

$\textrm{Tr}\left[B\right]$ represents the trace of operator $B$.

• •

$\rho\left({\cal L}\right)$ represents the space of density matrices, meaning the space of bounded operators acting on ${\cal H}$ with trace $1$ and positive.

• •

$|{\rho}\rangle\rangle$ is a vector in the Fock-Liouville space.

• •

$\langle\langle{A}|{B}\rangle\rangle=\textrm{Tr}\left[A^{\dagger}B\right]$ is the scalar product of operators $A,B\in B({\cal H})$ in the Fock-Liouville space.

• •

$\tilde{{\cal 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 sakurai_94 or Nielsen and Chuang’s Quantum Computing book nielsen_00 . For more advanced students, looking for a solid mathematical description of quantum mechanics methods, we recommend Galindo and Pascual galindo_pascual_90 . Finally, for a more philosophical discussion, you should go to Peres’ book peres_95 .

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 to any isolated physical system, there is a complex Hilbert space ${\cal 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 $(|{\psi}\rangle\in{\cal 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. nielsen_00 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 $\left\{|{\psi_{i}}\rangle\right\}$ with probabilities $p_{i}$. Therefore, our knowledge of the system is given by an ensemble of pure states described by the set $\left\{|{\psi_{i}}\rangle,\;p_{i}\right\}$. 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).

[TABLE]

Density matrices are bounded operators that fulfil two mathematical conditions

A density matrix $\rho$ has unit trace $\left(\textrm{Tr}[\rho]=1\right)$. 2. 2.

A density matrix is a positive matrix $\rho>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, $\textrm{Tr}[\rho^{2}]=\textrm{Tr}[\rho]=1$. Therefore, if $\textrm{Tr}[\rho^{2}]<1$ the system is mixed. The quantity $\textrm{Tr}[\rho^{2}]$ is called the purity of the states, and it fulfils the bounds $\frac{1}{d}\leq\textrm{Tr}[\rho^{2}]\leq 1$, being $d$ the dimension of the Hilbert space.

If we fix an arbitrary basis $\left\{|{i}\rangle\right\}_{i=1}^{N}$ of the Hilbert space the density matrix in this basis is written as $\rho=\sum_{i,j=1}^{N}\rho_{i,j}|{i}\rangle\!\langle{j}|$, or

[TABLE]

where the diagonal elements are called populations $\left(\rho_{ii}\in\mathbb{R}_{0}^{+}\text{ and }\sum_{i}\rho_{i,i}=1\right)$, while the off-diagonal elements are called coherences $\left(\rho^{\phantom{*}}_{i,j}\in\mathbb{C}\text{ and }\rho^{\phantom{*}}_{i,j}=\rho_{j,i}^{*}\right)$. Note that this notation is base-dependent.

Box 1. State of a two-level system (qubit)

The Hilbert space of a two-level system is just the two-dimension lineal space ${\cal 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: $\left\{|{0}\rangle,\;|{1}\rangle\right\}$. A pure state of the system would be any unit vector of ${\cal H}_{2}$. It can always be expressed as a $|{\psi}\rangle=a|{0}\rangle+b|{1}\rangle$ with $a,b\in\mathbb{C}$ s. t. $\left|a\right|^{2}+\left|b\right|^{2}=1$.

A mixed state is therefore represented by a positive unit trace operator $\rho\in O({\cal H}_{2})$.

$\rho=\begin{pmatrix}\rho_{00}&\rho_{01}\\ \rho_{10}&\rho_{11}\end{pmatrix}=\rho_{00}|{0}\rangle\!\langle{0}|+\rho_{01}|{0}\rangle\!\langle{1}|+\rho_{10}|{1}\rangle\!\langle{0}|+\rho_{11}|{1}\rangle\!\langle{1}|,$

(3)

ant it should fulfil $\rho_{00}+\rho_{11}=1$ and $\rho_{01}^{\phantom{*}}=\rho_{10}^{*}$.

Once we know the state of a system, it is natural to ask about the possible outcomes of experiments (see Ref. sakurai_94 , Section 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 form222For simplicity, we assume a non-degenerated spectrum.

[TABLE]

being $a_{i}\in\mathbb{R}$ the eigenvalues of the observable and $|{a_{i}}\rangle$ their corresponding eigenvectors. The probability of obtaining the result $a_{i}$ when measuring the property described by observable $O$ in a state $|{\psi}\rangle$ is given by

[TABLE]

After the measurement we obtain the state $|{a_{i}}\rangle$ if the outcome $a_{i}$ was measured. This is called the post-measurement state.

This postulate allow 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 $|{\psi}\rangle$ by just applying the simple formula

[TABLE]

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}}\rangle$ is

[TABLE]

and the expectation value of operator $O$ is

[TABLE]

Box 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}\rangle\!\langle{0}|+E_{1}|{1}\rangle\!\langle{1}|.$

(9)

This operator has two eigenvalues $\left\{E_{0},\;E_{1}\right\}$ with two corresponding eigenvectors $\left\{|{0}\rangle,\;|{1}\rangle\right\}$.

If we have a pure state $\psi=a|{0}\rangle+b|{1}\rangle$ the probability of measuring the energy $E_{0}$ would be $P(E_{0})=\left|\langle{0}|{\psi}\rangle\right|^{2}=\left|a\right|^{2}$. The probability of finding $E_{1}$ would be $P(E_{1})=\left|\langle{1}|{\psi}\rangle\right|^{2}=\left|b\right|^{2}$. The expected value of the measurement is $\langle H\rangle=E_{0}\left|a\right|^{2}+E_{1}\left|b\right|^{2}$.

In the more general case of having a mixed state $\rho=\rho_{00}|{0}\rangle\!\langle{0}|+\rho_{01}|{0}\rangle\!\langle{1}|+\rho_{10}|{1}\rangle\!\langle{0}|+\rho_{11}|{1}\rangle\!\langle{1}|$ the probability of finding the ground state energy is $P(0)=\textrm{Tr}\left[|{0}\rangle\!\langle{0}|\rho\right]=\rho_{00}$, and the expected value of the energy would be $\langle H\rangle=\textrm{Tr}\left[H\rho\right]=E_{0}\rho_{00}+E_{1}\rho_{11}$.

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 galindo_pascual_90 , Section 2.9).

Postulate 3

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

[TABLE]

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 $\hbar=1)$.

The Hamiltonian of a system is the operator corresponding to its energy, and it can be non-trivial to realise.

Schrödinger equation can be formally solved in the following way. If at $t=0$ the state of a system is given by $|{\psi(0)}\rangle$ at time $t$ it will be

[TABLE]

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 *

[TABLE]

with $U\in{\cal B}\left({\cal H}\right)$ s.t. $UU^{\dagger}=U^{\dagger}U=\mathbb{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 realisations like photons and cold atoms.

The proof of this lemma can be found at nielsen_00 .

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

[TABLE]

where we have used the commutator $\left[A,B\right]=AB-BA$, and ${\cal L}$ is the so-called Liouvillian superoperator.

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

[TABLE]

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.

Box 3. Time evolution of a two-level system.

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

$H_{\text{free}}=E_{0}|{0}\rangle\!\langle{0}|+E_{1}|{1}\rangle\!\langle{1}|,$

(15)

As the states $|{0}\rangle$ and $|{1}\rangle$ are Hamiltonian eigenstates if at $t=0$ the atom is at the excited state $|{\psi(0)}\rangle=|{1}\rangle$ after a time $t$ the state would be $|{\psi(t)}\rangle=e^{-iHt}|{1}\rangle=e^{-iE_{1}t}|{1}\rangle$.

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_{\text{free}}=E|{1}\rangle\!\langle{1}|,$

(16)

with $E\equiv E_{1}$. 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}\rangle\!\langle{1}|+\Omega\left(|{0}\rangle\!\langle{1}|+|{1}\rangle\!\langle{0}|\right),$

(17)

where $\Omega$ is the frequency of driving. By using the von-Neumann equation (13) we can calculate the populations $\left(\rho_{00},\rho_{11}\right)$ 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.

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

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 ${\cal H}={\cal H}_{1}\otimes{\cal H}_{2}\otimes\cdots\otimes{\cal H}_{N}$. The state of the composite physical system is given by a unit vector of ${\cal H}$. Moreover, if each subsystem belonging to ${\cal H}_{i}$ is prepared in the state $|{\psi_{i}}\rangle$ the total state is given by $|{\psi}\rangle=|{\psi_{1}}\rangle\otimes|{\psi_{2}}\rangle\otimes\cdots\otimes|{\psi_{N}}\rangle$.

The symbol $\otimes$ 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 $\rho_{i}$ the total state is given by $\rho=\rho_{1}\otimes\rho_{2}\otimes\cdots\otimes\rho_{N}$.

States that can be expressed in the simple form $|{\psi}\rangle=|{\psi_{1}}\rangle\otimes|{\psi_{2}}\rangle$, in any specific basis, are very particular and they are called separable states (For this discussion, we use a bipartite system as an example. The extension to a general multipartite system is straightforward.) . In general, any arbitrary state should be described as $|{\psi}\rangle=\sum_{i,j}|{\psi_{i}}\rangle\otimes|{\psi_{j}}\rangle$ (or $\rho=\sum_{i,j}\rho_{i}\otimes\rho_{j}$ for mixed states). Non-separable 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 ${\cal H}={\cal H}_{a}\otimes{\cal 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 $\rho$ the reduced density operator of the subsystem $a$ is defined by the operator

[TABLE]

were $\textrm{Tr}_{b}$ is the partial trace over subspace $b$ and it is defined as nielsen_00

[TABLE]

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.

Box 4. Two two-level atoms

If we have two two-level systems, the total Hilbert space is given by ${\cal H}={\cal H}_{2}\otimes{\cal H}_{2}$. A basis of this Hilbert space would be given by the set $\left\{|{00}\rangle\equiv|{0}\rangle_{1}\otimes|{0}\rangle_{2},\;|{01}\rangle\equiv|{0}\rangle_{1}\otimes|{1}\rangle_{2},\;|{10}\rangle\equiv|{1}\rangle_{1}\otimes|{0}\rangle_{2},\;|{11}\rangle\equiv|{1}\rangle_{1}\otimes|{1}\rangle_{2}\right\}$. If both systems are in their ground state, we can describe the total state by the separable vector

$|{\psi}\rangle_{G}=|{00}\rangle.$

(20)

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

$\displaystyle|{\psi}\rangle_{S}$ $\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2}}\left(|{0}\rangle_{1}+|{1}\rangle_{1}\right)\otimes\frac{1}{\sqrt{2}}\left(|{0}\rangle_{2}+|{1}\rangle_{2}\right)$

(21)

$\displaystyle=$ $\displaystyle\frac{1}{2}\left(|{00}\rangle+|{10}\rangle+|{01}\rangle+|{11}\rangle\right)$

An entangled state would be

$|{\psi}\rangle_{E}=\frac{1}{\sqrt{2}}\left(|{00}\rangle+|{11}\rangle\right).$

(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 $|{\psi}\rangle_{E}$.

$\rho_{E}=|{\psi}\rangle\langle{\psi}|_{E}=\frac{1}{2}\left(|{00}\rangle\!\langle{00}|+|{00}\rangle\!\langle{11}|+|{11}\rangle\!\langle{00}|+|{11}\rangle\!\langle{11}|\right).$

(23)

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

$\rho_{E}^{(1)}=\langle{0}|_{2}\rho_{E}|{0}\rangle_{2}+\langle{1}|_{2}\rho_{E}|{1}\rangle_{2}=\frac{1}{2}\left(|{00}\rangle\!\langle{00}|_{1}+|{11}\rangle\!\langle{11}|_{2}\right).$

(24)

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

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 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. This allows us to define a linear space of matrices, converting the matrices effectively into vectors ($\rho\to|{\rho}\rangle\rangle$). This is called Fock-Liouville space (FLS). The usual definition of the scalar product of matrices $\phi$ and $\rho$ is defined as $\langle\langle{\phi}|{\rho}\rangle\rangle\equiv\textrm{Tr}\left[\phi^{\dagger}\rho\right]$. The Liouville super-operator 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.

Box 5. Time evolution of a two-level system.

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

$|{\rho}\rangle\rangle=\begin{pmatrix}\rho_{00}\\ \rho_{01}\\ \rho_{10}\\ \rho_{11}\end{pmatrix}.$

(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{{\cal L}}=\left(\begin{array}[]{cccc}0&i\Omega&-i\Omega&0\\ i\Omega&iE&0&-i\Omega\\ -i\Omega&0&-iE&i\Omega\\ 0&-i\Omega&i\Omega&0\end{array}\right),

(26)

where each row is calculated just by observing the output of the operation $-i\left[H,\rho\right]$ in the computational basis of the density matrices space. The time evolution of the system now corresponds to the matrix equation $\frac{d|{\rho}\rangle\rangle}{dt}=\tilde{{\cal L}}|{\rho}\rangle\rangle$, that in matrix notation would be

\begin{pmatrix}\dot{\rho}_{00}\\ \dot{\rho}_{01}\\ \dot{\rho}_{10}\\ \dot{\rho}_{11}\end{pmatrix}=\left(\begin{array}[]{cccc}0&i\Omega&-i\Omega&0\\ i\Omega&iE&0&-i\Omega\\ -i\Omega&0&-iE&i\Omega\\ 0&-i\Omega&i\Omega&0\end{array}\right)\begin{pmatrix}\rho_{00}\\ \rho_{01}\\ \rho_{10}\\ \rho_{11}\end{pmatrix}

(27)

V CPT-maps and the Lindblad Master Equation.

V.1 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 way, by a coherent evolution given (Postulate 3) and by collapsing after a measurement (Postulate 2). Many efforts have been made to unify these two ways of evolving schlosshauer_07 , 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 $\rho({\cal H})$ as the space of all density matrices in the Hilbert space ${\cal H}$. Therefore, we are looking for a map of this space onto itself, ${\cal V}:\rho({\cal H})\to\rho({\cal H})$. To ensure that the output of the map is a density matrix this should fulfil the following properties

• •

Trace preserving. $\textrm{Tr}\left[{\cal V}A\right]=\textrm{Tr}\left[A\right],$ $\forall A\in O({\cal H})$.

• •

Completely positive (see below).

Any map that fulfils 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 ${\cal V}$ is positive iff $\forall A\in B({\cal H})$ s.t. $A\geq 0\Rightarrow{\cal V}A\geq 0$.

This definition is based in 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 ${\cal V}$ is completely positive iff $\forall n\in\mathbb{N}$, ${\cal V}\otimes\mathbb{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 $|{\psi_{B}}\rangle=\frac{1}{\sqrt{2}}\left(|{01}\rangle+|{10}\rangle\right)$ its density matrix can be expressed as

[TABLE]

with a matrix representation

[TABLE]

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

[TABLE]

and it is positive.

It is easy to check that the transformation $\mathbb{1}\otimes T_{2}$, meaning that we transpose the matrix of the second subsystem leads to a non-positive matrix

[TABLE]

The total matrix is

[TABLE]

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

V.2 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 like Breuer and Petruccione’s breuer_02 as well as Gardiner and Zoller’s gardiner_00 . Here, we follow the derivation presented in Ref. manzano:av18 . Our initial point is displayed in Figure 3. A total system belonging to a Hilbert space ${\cal H}_{T}$ is divided into our system of interest, belonging to a Hilbert space ${\cal H}$, and the environment living in ${\cal H}_{E}$.

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

[TABLE]

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 $\rho(t)=\textrm{Tr}_{E}[\rho_{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\mathbb{1}_{E}+\mathbb{1}_{S}\otimes H_{E}+\alpha H_{I}$, with $H\in{\cal H}$, $H_{E}\in{\cal H}_{E}$, and $H_{I}\in{\cal H}_{T}$, and being $\alpha$ 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

[TABLE]

with $S_{i}\in B({\cal H}$) and $E_{i}\in B({\cal H}_{E})$333From now on we will not writethe identity operators of the Hamiltonian parts explicitly when they can be inferred from the context..

To better describe the dynamics of the system, it is useful to work in the interaction picture (see Ref. galindo_pascual_90 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{\cal B}({\cal H}_{T})$ is represented in this picture by the time-dependent operator $\hat{O}(t)$, and its time evolution is

[TABLE]

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

[TABLE]

This equation can be easily integrated to give

[TABLE]

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 troublesome the fact that the state $\tilde{\rho}(t)$ depends on the integration of the density matrix in all previous time. To avoid that we can introduce Eq. (71) into Eq. (70) giving

[TABLE]

By applying this method one more time we obtain

[TABLE]

After this substitution, the integration of the previous states of the system is included only in the terms that are $O(\alpha^{3})$ 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. (73). Under this approximation we have

[TABLE]

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

[TABLE]

This is not a closed time-evolution equation for $\hat{\rho}(t)$, because the time derivative still depends on the full density matrix $\hat{\rho}_{T}(t)$. To proceed, we need to make two more assumptions. First, we assume that $t=0$ the system and the environment have a separable state in the form $\rho_{T}(0)=\rho(0)\otimes\rho_{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 $\rho_{E}(0)=\exp\left(-H_{E}/T\right)/\textrm{Tr}[\exp\left(-H_{E}/T\right)]$, being $T$ the temperature and taking the Boltzmann constant as $k_{B}=1$. By using these assumptions, and the expansion of $H_{I}$ (68), we can calculate an expression for the first element of the r.h.s of Eq. (75).

[TABLE]

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

[TABLE]

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 non-correlated 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 timescales of correlation ($\tau_{\text{corr}}$) and relaxation of the environment ($\tau_{\text{rel}}$) are much smaller than the typical system timescale ($\tau_{\text{sys}}$), as the coupling strength is very small ($\alpha<<$). Therefore, under this strong assumption, we can assume that the environment state is always thermal and is decoupled from the system state, $\hat{\rho}_{T}(t)=\hat{\rho}(t)\otimes\hat{\rho}_{E}(0)$. Eq. (77) then transforms into

[TABLE]

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 realising that the kernel in the integration 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\rightarrow t-s$, we obtain the famous Redfield equation redfield:IBM57.

[TABLE]

It is known that this equation does not warrant the positivity of the map, and it sometimes gives rise to density matrices that are non-positive. 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\left[H,A\right]$, $\forall A\in{\cal B}({\cal H})$. The eigenvectors of this superoperator form a complete basis of space ${\cal B}({\cal H})$ and, therefore, we can expand the system-environment operators from Eq. (68) in this basis

[TABLE]

where the operators $S_{i}(\omega)$ fulfils

[TABLE]

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

[TABLE]

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 $\hat{S}_{k}=e^{itH}S_{k}e^{-itH}$. By using the eigen-expansion (80) we arrive to

[TABLE]

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

[TABLE]

We now apply the eigenvalue decomposition in terms of $S_{k}(\omega)$ for $\hat{H}_{I}(t-s)$ and in terms of $S_{k}^{\dagger}(\omega^{\prime})$ for $\hat{H}_{I}(t)$. By using the permutation property of the trace and the fact that $\left[H_{E},\rho_{E}(0)\right]=0$, and after some non-trivial algebra we obtain

[TABLE]

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

[TABLE]

where we are writing the environment operators of the interaction Hamiltonian in the interaction picture ($\hat{E}_{l}(t)=e^{iH_{E}t}E_{l}e^{-iH_{E}t}$). At this point, we can already perform the rotating wave approximation. By considering the time-dependency on Eq. (85), we conclude that the terms with $\left|\omega-\omega^{\prime}\right|>>\alpha^{2}$ will oscillate much faster than the typical timescale of the system evolution. Therefore, they do not contribute to the evolution of the system. In the low-coupling regime $(\alpha\rightarrow 0)$ we can consider that only the resonant terms, $\omega=\omega^{\prime}$, contribute to the dynamics and remove all the others. By applying this approximation to Eq. (85) reduces to

[TABLE]

To divide the dynamics into Hamiltonian and non-Hamiltonian we now decompose the operators $\Gamma_{kl}$ into Hermitian and non-Hermitian parts, $\Gamma_{kl}(\omega)=\frac{1}{2}\gamma_{kl}(\omega)+i\pi_{kl}$, with

[TABLE]

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

[TABLE]

The Hamiltonian dynamics now is influenced by a term $H_{Ls}=\sum_{\omega,k,l}\pi_{kl}(\omega)S_{k}^{\dagger}(\omega)S_{l}(\omega)$. 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. Eq. (89) 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 $\gamma_{kl}(\omega)$ is positive as they are the Fourier’s transform of a positive function $\left(\textrm{Tr}\left[\hat{E}_{k}^{\dagger}(s)E_{l}\hat{\rho}_{E}(0)\right]\right)$. Therefore, this matrix can be diagonalised. This means that we can find a unitary operator, $O$, s.t.

[TABLE]

We can now write the master equation in a diagonal form

[TABLE]

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

[TABLE]

The operators $L_{i}$ are usually referred to as jump operators.

V.3 Derivation of the Lindblad Equation as a CPT generator

The second way of deriving 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 analyse this problem following mainly Ref. wilde_17 .

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

Lemma 2

Any map ${\cal V}:B\left({\cal H}\right)\to B\left({\cal H}\right)$ that can be written in the form ${\cal V}\rho=V^{\dagger}\rho V^{\phantom{\dagger}}$ with $V\in B\left({\cal H}\right)$ is positive.

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

Proof.

If $\rho\geq 0\Rightarrow\rho=A^{\dagger}A^{\phantom{\dagger}}$, with $A\in B({\cal H})$. Therefore, ${\cal V}\rho=V^{\dagger}\rho V^{\phantom{\dagger}}\Rightarrow\langle{\psi}|V^{\dagger}\rho V|{\psi}\rangle=\langle{\psi}|V^{\dagger}A^{\dagger}AV|{\psi}\rangle=\left|\left|AV|{\psi}\rangle\right|\right|\geq 0$. Therefore, if $\rho$ 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 next theorem.

Theorem 1

Choi’s Theorem.

A linear map ${\cal V}:B({\cal H})\to B({\cal H})$ is completely positive iff it can be expressed as

[TABLE]

with $V_{i}\in B({\cal H})$.

The proof of this theorem requires some algebra.

Proof

The ‘if’ implication is a trivial consequence of the previous lemma. 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, ${\cal H}$, we define a new Hilbert space of the same dimension ${\cal H}_{A}$.

We define a maximally entangled pure state in the bipartition ${\cal H}_{A}\otimes{\cal H}$ in the way

[TABLE]

being $\left\{|{i}\rangle\right\}$ and ${\left\{|{i}\rangle_{A}\right\}}$ arbitrary orthonormal bases for ${\cal H}$ and ${\cal H}_{A}$.

We can extend the action of our original map ${\cal V}$, that acts on ${\cal B}({\cal H})$ to our extended Hilbert space by defining the map ${\cal V}_{2}:{\cal B}({\cal H}_{A})\otimes{\cal B}({\cal H})\to{\cal B}({\cal H}_{A})\otimes{\cal B}({\cal H})$ as

[TABLE]

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 ${\cal 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 (94), obtaining

[TABLE]

Now we can use the maximal entanglement of the state $|{\Gamma}\rangle$ to relate the original map ${\cal V}$ and the action ${\cal V}_{2}|{\Gamma}\rangle\!\langle{\Gamma}|$ by taking the matrix elements with respect to ${\cal H}_{A}$.

[TABLE]

To relate this operation to the action of the map to an arbitrary vector $|{\psi}\rangle\in{\cal H}_{A}\otimes{\cal H}$, we can expand it in this basis as

[TABLE]

We can also define an operator $V_{|{\psi}\rangle}\in{\cal B}\left({\cal H}\right)$ s.t. it transforms $|{\Gamma}\rangle$ into $|{\psi}\rangle$. Its explicit action would be written as

[TABLE]

At this point, we have related the vectors in the extended space ${\cal H}_{A}\otimes{\cal H}$ to operators acting on ${\cal H}$. This can only be done because the vector $|{\Gamma}\rangle$ is maximally entangled. We go now back to our extended map ${\cal V}_{2}$. Its action on $|{\Gamma}\rangle\!\langle{\Gamma}|$ is given by Eq. (96) and as it is a positive map it can be expanded as

[TABLE]

with $|{v_{l}}\rangle\in{\cal H}_{A}\otimes{\cal H}$. The vectors $|{v_{l}}\rangle$ can be related to operators in ${\cal H}$ as in Eq. (99).

[TABLE]

Based on this result we can calculate the product of an arbitrary vector $|{i}\rangle_{A}\in{\cal H}_{A}$ with $|{v_{l}}\rangle$.

[TABLE]

This is the last ingredient we need for the proof.

We come back to the original question, we want to characterise the map ${\cal V}$. We do so by applying it to an arbitrary basis element $|{i}\rangle\!\langle{j}|$ of ${\cal B}\left({\cal H}\right)$.

[TABLE]

As $|{i}\rangle\!\langle{j}|$ is an arbitrary element of a basis any operator can be expanded in this basis. Therefore, it is straightforward to prove that

[TABLE]

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 ${\cal V}:B({\cal H})\to B({\cal H})$ is completely positive and trace-preserving iff it can be expressed as

[TABLE]

with $V_{l}\in B({\cal H})$ fulfilling

[TABLE]

Proof.

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

[TABLE]

We have to prove also that any map in the form (104) is trace-preserving only if the operators $V_{l}$ fulfil (105). We start by stating that if the map is trace-preserving by applying it to an any arbitrary element of a basis of ${\cal B}\left({\cal H}\right)$ we should obtain

[TABLE]

As the map has a form given by (104) we can calculate this same trace in an alternative way.

[TABLE]

where $\left\{|{k}\rangle\right\}$ is an arbitrary basis of ${\cal H}$. As both equalities should be right we obtain

[TABLE]

and therefore, the condition (105) should be fulfilled.

End of the proof.

Operators $V_{i}$ of a map fulfilling condition (105) 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 fulfil relation (105) 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 ${\cal L}$ such that

[TABLE]

and therefore our CPT-map could be expressed as ${\cal V}(t)=e^{{\cal L}t}$. The following calculation is about founding the explicit expression of ${\cal L}$. We start by choosing an orthonormal basis of the bounded space of operators ${\cal B}({\cal H})$, $\left\{F_{i}\right\}_{i=1}^{d^{2}}$. To be orthonormal it should satisfy the following condition

[TABLE]

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}}\mathbb{1}_{{\cal H}}$. It is trivial to prove that the norm of this element is one, and it is easy to see from Eq. (111) that all the other elements of the basis should have trace zero.

[TABLE]

The closure relation of this basis is $\mathbb{1}_{{\cal B}({\cal H})}=\sum_{i}|{F_{i}}\rangle\rangle\!\langle\langle{F_{i}}|$. Therefore, the Krauss operators can be expanded in this basis by using the Fock-Liouville notation

[TABLE]

As the map ${\cal V}(t)$ is in the form (93) we can apply (113) to obtain444 For simplicity, in this discussion we omit the explicit time-dependency of the density matrix..

[TABLE]

where we have absorved the sumation over the Krauss operators in the terms $c_{i,j}(t)=\sum_{l}\langle\langle{F_{i}}|{V_{l}}\rangle\rangle\langle\langle{V_{l}}|{F_{j}}\rangle\rangle$. We go back now to the original problem by applying this expansion into the time-derivative of Eq. (110)

[TABLE]

where we have separated the summations to take into account that $F_{d^{2}}=\frac{1}{\sqrt{d}}\mathbb{1}_{{\cal H}}$. By using this property this equation simplifies to

[TABLE]

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

[TABLE]

Introducing these coefficients in Eq (116) we obtain an equation with no explicit dependence in time.

[TABLE]

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

[TABLE]

Applying it to Eq. (116).

[TABLE]

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

[TABLE]

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

[TABLE]

We define now the last operator for this proof, $G_{2}\equiv G+\frac{g_{d^{2},d^{2}}}{2d}$, and the expression of the time derivative leads to

[TABLE]

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

[TABLE]

Therefore, $G_{2}$ should fulfil

[TABLE]

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

[TABLE]

Finally, by definition the coefficients $g_{i,j}$ can be arranged to form a Hermitian, and therefore diagonalisable, matrix. By diagonalising it, we obtain the diagonal form of the Lindblad master equation.

[TABLE]

V.4 Properties of the Lindblad Master Equation

Some interesting properties of the Lindblad equation are:

• •

Under a Lindblad dynamics, if all the jump operators are Hermitian, the purity of a system fulfils $\frac{d}{dt}\left(\textrm{Tr}\left[\rho^{2}\right]\right)\leq 0$. The proof is given in A.

• •

The Lindblad Master Equation is invariant under unitary transformations of the jump operators

[TABLE]

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

[TABLE]

where $a_{i}\in\mathbb{C}$ and $b\in\mathbb{R}$. The proof of this can be found in Ref. breuer_02 (Section 3).

• •

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

Box 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 state555This 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 state of the atom, we can derive a Markovian master equation for the reduced density matrix of the atom (see for instance Refs. breuer_02 ; gardiner_00 ). The master equation we will study is

$\displaystyle\frac{d}{dt}\rho(t)=-i\left[H,\rho\right]+\Gamma\left(\sigma^{-}\rho\sigma^{+}-\frac{1}{2}\left\{\sigma^{+}\sigma^{-},\rho\right\}\right),$

(130)

where $\Gamma$ is the coupling between the atom and the vacuum.

In the Fock-Liouvillian space (following the same ordering as in Eq. (3)) the Liouvillian corresponding to evolution (130) is

{\cal L}=\left(\begin{array}[]{cccc}0&i\Omega&-i\Omega&\Gamma\\ i\Omega&-iE-\frac{\Gamma}{2}&0&-i\Omega\\ -i\Omega&0&-iE-\frac{\Gamma}{2}&i\Omega\\ 0&-i\Omega&i\Omega&-\Gamma\\ \end{array}\right).

(131)

Expressing explicitly the set of differential equations we obtain

$\displaystyle\dot{\rho}_{00}$ $\displaystyle=$ $\displaystyle i\Omega\rho_{01}-i\Omega\rho_{10}+\Gamma\rho_{11}$

$\displaystyle\dot{\rho}_{01}$ $\displaystyle=$ $\displaystyle i\Omega\rho_{00}-\left(iE-\frac{\Gamma}{2}\right)\rho_{01}-i\Omega\rho_{11}$

$\displaystyle\dot{\rho}_{10}$ $\displaystyle=$ $\displaystyle-i\Omega\rho_{00}\left(-iE-\frac{\Gamma}{2}\right)\rho_{10}+i\Omega\rho_{11}$

(132)

$\displaystyle\dot{\rho}_{10}$ $\displaystyle=$ $\displaystyle-i\Omega\rho_{01}+i\Omega\rho_{10}-\Gamma\rho_{11}$

VI Resolution of the Lindblad Master Equation

VI.1 Integration

To calculate the time evolution of a system determined by a Master Equation in the form (130) 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 briegel:pra93 . In most cases, we have to rely on numerical approximated methods. One of the most popular approaches is the $4^{th}$-order Runge-Kutta algorithm (see, for instance, numericalrecipes 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 $\left(t\rightarrow\infty\right)$. 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 non-physical results after a finite time. One solution to this problem is the use of algorithms specifically designed to preserve the trace, as Crank-Nicholson algorithm goldberg:ajp67 . The problem with this kind of algorithms is that they consume more computational power than Runge-Kutta, and therefore they are not useful to calculate the long-time behaviour of big systems. An analysis of different methods and their advantages and disadvantages can be found at Ref. riesch:jcp19 .

Box 7. Time dependency of the two-level system with decay.

In this box we show some results of solving Eq (130) and calculating the density matrix as a function of time. A Mathematica notebook solving this problem can be found at notebook . To illustrate the time behaviour 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 $\rho_{11}=1$, with no coherence between different states, meaning $\rho_{01}=\rho_{10}=0$. If the decay parameter $\Gamma$ is equal to zero, the problem reduces to solve von Neumann equation, and the result is displayed in Figure 2. The other extreme case would be a system with no coherent dynamics ($\Omega=0$) but with decay. In this case, we observe an exponential decay of the population of the excited state. Finally, we can calculate the dynamics of a system with both coherent driving and decay. In this case, both behaviours coexist, and there are oscillations and decay.

VI.2 Diagonalisation

As we have discussed before, in the Fock-Liouville space the Liouvillian corresponds to a complex matrix (in general complex, non-hermitian, and non-symmetric). By diagonalising 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 diagonalising. This is due to the high dimensionality of the Liouvillian that makes the diagonalisation process very costly in computing power. On the other hand, in order to calculate the steady-state, the diagonalisation is the most used method due to the problems of integrating the equation of motions discussed in the previous section.

Let see first how we use diagonalisation 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 $\Lambda_{i}$ we can obtain the eigenvectors $|{\Lambda_{i}^{R}}\rangle\rangle$ and $|{\Lambda_{i}^{L}}\rangle\rangle$ s. t.

[TABLE]

An arbitrary system can be expanded in the eigenbasis of $\tilde{{\cal L}}$ as thingna:sr16 ; gardiner_00

[TABLE]

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

[TABLE]

Note that in this case to calculate the state a time $t$ we do not need to integrate into the interval $\left[0,t\right]$, 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 behaviour. 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\to\infty$.

For any finite system, Evans’ Theorem ensures the existence of at least one zero eigenvalue of the Liouvillian matrix evans:cmp77 ; evans:jfa79 . 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 buca:njp12 ; manzano:prb14 ; manzano:av18 . This is a non-generic case, and for most purposes, we can assume the existence of a unique fixed point in the dynamics of the system. Therefore, diagonalising 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 analysed in several recent papers albert:pra14 ; thingna:sr16 .

Box 8. Spectrum-analysis of the Liouvillian for the two-level system with decay.

Here we diagonalise (131) and obtain its steady state. A Mathematica notebook solving this problem can be downloaded from notebook . 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 5 sisplays the spectrum of the superoperator ${\cal L}$ given in (131).

As there 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:

\rho_{SS}=\left(\begin{array}[]{cc}\frac{(1+n)\left(4\,E^{2}+(\Gamma+2n\,\Gamma)^{2}\right)+4(1+2n)\Omega^{2}}{(1+2n)\left(4\,E^{2}+(\Gamma+2n\,\Gamma)^{2}+8\Omega^{2}\right)}&\frac{2(-2\,E-i(\Gamma+2n\Gamma))\Omega}{(1+2n)\left(4\,E^{2}+(\Gamma+2n\,\Gamma)^{2}+8\,\Omega^{2}\right)}\\ \frac{2(-2\,E+i(\Gamma+2n\,\Gamma))\Omega}{(1+2n)\left(4\,E^{2}+(\Gamma+2n\,\Gamma)^{2}+8\Omega^{2}\right)}&\frac{n\left(4E^{2}+(\Gamma+2n\Gamma)^{2}\right)+4(1+2n)\Omega^{2}}{(1+2n)\left(4\,E^{2}+(\Gamma+2n\Gamma)^{2}+8\,\Omega^{2}\right)}\\ \end{array}\right)

(136)

VII Acknowledgements

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

Appendix A Proof of $\frac{d}{dt}\textrm{Tr}\left[\rho^{2}\right]\leq 0$

In this appendix we proof that under the Lindblad dynamics given by Eq. (127) the purity of a density matrix fulfils that $\frac{d}{dt}\textrm{Tr}\left[\rho^{2}\right]\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

[TABLE]

where we have used the cyclic property of the trace operator666This property is used along all the demonstration without explicitly mentioning it.. By inserting the Lindblad Eq. (127) into the r.h.s of (137) we obtain

[TABLE]

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

[TABLE]

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

[TABLE]

where we assume an ordering of the eigenvalues in the form $\Lambda_{0}\geq\Lambda_{1}\geq\cdots\geq\Lambda_{d}$.

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

[TABLE]

where we have introduced the oefficients $x_{ij}^{(k)}\equiv\left|\langle{\Lambda_{i}}|\tilde{L}_{k}|{\Lambda_{j}}\rangle\right|^{2}$. As the operators $L_{k}$ are Hermitian these coefficients fulfil $x_{ij}^{(k)}=x_{ji}^{(k)}$

The second term is expanded as

[TABLE]

where we have used the closure relation in the density matrix eigenbasis, $\mathbb{1}_{{\cal H}}=\sum_{j}|{\Lambda_{j}}\rangle\!\langle{\Lambda_{j}}|$. The inequality can be written now as

[TABLE]

As $x_{ij}=x_{ji}$ we can re-order the $ij$ sum in the following way

[TABLE]

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

[TABLE]

It is obvious that (145) $\Rightarrow$ (144) (but not the other way around). The inequalities (145) are easily proved just by taking into account that $x_{ij}^{(k)}\geq 0$ and applying the Triangular Inequality.