Open Quantum Systems/Dissipative preparation of many-body states

In this chapter, we will see how one can use the interplay between coherent and dissipative dynamics can give rise to the dissipative preparation of quantum many-body states with useful and interesting properties. For simplicity, we focus on Markovian master equations in Lindblad form, given by

$\frac{d}{dt}\rho =-\frac{i}{\hslash }[H,\rho ]+\sum _i{\gamma }_i(c_i\rho c_i^{†}-\frac{1}{2}\{c_i^{†}{c}_i,\rho \}).$

Stationary states

The stationary states of a master equation are found by solving the equation $\text{Re}(d/dt\rho )=0$. In most cases, the stationary states do not have an imaginary component of the dynamics so this is equivalent to solving the equation $d/dt\rho =0$. The detailed nature of the stationary state of course depends on the details of the master equation, but we can make a few general statements.

For example, consider the case where all jump operators are Hermitian, i.e., $c_i^{†}=c_i$. This corresponds to a pure dephasing where only the off-diagonal elements of the density matrix get damped out, but there is no reshuffling of the diagonal elements (probabilities) taking place. In this case, one can show that the maximally mixed state given by

$\rho =\left(\begin{array}{c}1/d\\ & 1/d\\ & & \ddots \end{array}\right),$

with $d$ being the Hilbert space dimension, is a stationary state of the dynamics. To prove this, we first note that the commutator with the Hamiltonian vanishes since the maximally mixed state is proportional to the identity. The dissipative part can be written as

$\frac{d}{dt}\rho =\sum _i\frac{{\gamma }_i}{d}(c_i{c}_i^{†}-c_i^{†}{c}_i),$

which also vanishes in the case of all jump operators being Hermitian.

Dark states

On the other hand, there can be situations where the stationary state is a pure state, i.e., $\rho =|\psi ⟩⟨\psi |$. To see how this comes about, we first introduce an effective non-Hermitian Hamiltonian,

$H^{\prime }=H-\frac{i\hslash }{2}\sum _i{c}_i^{†}{c}_i.$

Then, we can write the Lindblad master equation as

$\frac{d}{dt}\rho =-\frac{i}{\hslash }(H^{\prime }\rho -\rho H{\text{'}}^{†})+\sum _i{\gamma }_i{c}_i\rho c_i^{†}.$

Then, to have a pure stationary state $|\psi ⟩$ two conditions have to be fulfilled :

1. $|\psi ⟩$ is an eigenstate of the effective Hamiltonian, i.e., $H^{\prime }|\psi ⟩={\epsilon }_{\psi }|\psi ⟩$.
2. $|\psi ⟩$ is an eigenstate of all jump operators. One can always re-define the jump operators such that one has $c_i|\psi ⟩=0$ for all $i$ without changing the master equation.

As an example, consider a three level atom interacting with the radiation field, consisting of two electronic ground states (e.g., different hyperfine states) and one electronically excited state that can decay into both ground states. In this case, the Hamiltonian in the rotating wave approximation is given by

$H=\frac{\hslash \mathrm{\Omega }}{2}(|1⟩⟨2|+|3⟩⟨2|+\text{H.c.}),$

where we have assumed that the Rabi frequencies on both transitions are identical. The jump operators take the form

$\begin{array}{cc}{c}_1& =|1⟩⟨2|\\ c_2& =|3⟩⟨2|,\end{array}$

and again we assume that both decay rates have the same value, $\gamma$. Then, we find that the antisymmetric combination $|\psi ⟩=(|1⟩-|3⟩)/\sqrt{2}$ is a stationary state of the quantum optical master equation. It is an eigenstate of the Hamiltonian with eigenvalue $0$ as the two excitation paths to the excited state interfere destructively. It is also annihilated by both jump operators as it does not contain any contribution from the state $|2⟩$. As this state does not couple to the radiation field anymore, it is often called a "dark state".

Importantly, the dark state is also the only stationary state of the master equation. Any density matrix that has a contribution from ${\rho }_2=|2⟩⟨2|$ will couple to the dark state as the quantity

$⟨\psi |ℒ{\rho }_2|\psi ⟩=\gamma ⟨\psi |(c_1{\rho }_2{c}_1^{†}+c_2{\rho }_2{c}_2^{†})|\psi ⟩=\gamma$

is nonzero. Hence, the only remaining possiblity is the symmetric state $|{\psi }_s=(|1⟩+|3⟩)/\sqrt{2}$. However, $|{\psi }_s⟩$ is not an eigenstate of the Hamiltonian as we find

$H|{\psi }_s⟩=\hslash \mathrm{\Omega }|2⟩.$

Consequently, there are no other stationary density matrices besides $\rho =|\psi ⟩⟨\psi |$. This phenomenon is known as "coherent population trapping" and can also be generalized to unequal decay rates .

Dissipative preparation of a Bose condensate

We will now see how the concept of coherent population trapping can also be generalized to many-body systems. Consider a system of $N$ bosons in a lattice. Then, in the formalism of second quantization, we can assign a bosonic creation operator, $a_i^{†}$ to each site $i$. The Hilbert space of a single state is a Fock space encoding the number of bosons located at the site, i.e.,

$|n{⟩}_i=\frac{(a_i^{†}{)}^n}{\sqrt{n!}}|0{⟩}_i$,

where $|0{⟩}_i$ is the vacuum state of zero particles that is destroyed by the annihilation operator, i.e., $a_i|0⟩=0$. The total Hilbert space of the system is a product space of all the Fock spaces, with the commutation relation for bosons being

$[a_i,a_j^{†}]={\delta }_{ij}.$

With the help of a Fourier transform, we can also define creation and annihilation operators in momentum space,

$\begin{array}{cc}{a}_k^{†}& =\sum _j\frac{\exp (ikj)}{\sqrt{V}}{a}_j^{†}\\ a_k& =\sum _j\frac{\exp (ikj)}{\sqrt{V}}{a}_j,\end{array}$

where $V$ is the volume of the system.

In such a lattice system without interactions, the kinetic part of the Hamiltonian is given by the discretization of the Laplace operator,

$H=J\sum _{⟨ij⟩}{a}_i{a}_j^{†}+\text{H.c},$

where $J$ denotes the strength of the tunneling between adjacent lattice sites. Its ground state is a Bose-Einstein condensate (BEC) of all particles put into the mode with zero momentum, i.e.,

$|\psi ⟩=\frac{{\left(a_{k=0}^{†}\right)}^N}{\sqrt{N!}}|0⟩.$

This state characterizes a distinct state of matter exhibiting long-range phase coherence of the form

$\lim {\mathrm{\backslash limits}}_{|i-j|\to \mathrm{\infty }}⟨\psi |a_i^{†}{a}_j|\psi ⟩=\frac{N}{V}.$

The construction of suitable jump operators such that $|\psi ⟩$ is a dark state can in general be split up into a product $AB$. Here, $B$ is the "interrogation part" that checks whether the system is already in the correct state. Hence, we need to satisfy $B|\psi ⟩=0$ and obtain a finite value for all other states. The "pump part" $A$ then reshuffles the populations. The precise choice of $A$ is usually not important, as long as it does not create a subspace in the Hilbert space, from which a decay into the dark state is no longer possible.

It might be tempting to choose $B$ to involve the projection onto the dark state, i.e.,

$B=1-|\psi ⟩⟨\psi |.$

However, this jump operator contains terms that involve a product of operators involving sites separated by an arbitrary distance, so such a jump operator is unphysical. Nevertheless, it is possible to achieve the desired result using jump operators involving only adjacent sites $i$ and $j$ , by making the choice

$c_{ij}=(a_i^{†}+a_j^{†})(a_i-a_j).$

This jump operator can be understood as its interrogation part probing whether the phase is constant between adjacent sites. This makes sense because a BEC is a state of matter with long-range phase coherence. The pump part ensures that a phase fluctuation with nonzero $a_i-a_j$ will be projected onto a phase coherent state. To show that the BEC state is indeed a dark state of the dynamics, it is instructive to look at the Fourier transform of the jump operators ,

$c_k=\frac{1}{\sqrt{V}}\sum _q[1+e^{i(q-k)}][1-e^{-ik}]a_{q-k}^{†}{a}_q.$

From this one can see that the state with $k=0$ (i.e., the BEC) is the only dark state of the dynamics.

Frustration-free Hamiltonians

It is now very interesting, of course, to ask what kind of states we can prepare dissipatively. It should be noted that it is not possible to efficiently prepare arbitrary states because the preparation process might take a time that grows exponentially with the system size (e.g., in glassy systems) . A particularly important class of models whose ground states can be prepared efficiently, are frustration-free Hamiltonians . Such Hamiltonian are of the form

$H=\sum _{\lambda }{H}_{\lambda },$

where the parts $H_{\lambda }$ are projectors, i.e, $H_{\lambda }^2=H_{\lambda }$, and they act only on a local subspace of the system. Additionally, while the parts do not have to commute in general, they have to commute within the ground state manifold, i.e.,

$⟨{\psi }_0|[H_{\lambda },H_{\mu }]|{\psi }_0⟩,$

where $|{\psi }_0⟩$ is the ground state that can be found from the simultaneous minimization of all sub-parts, i.e., $H_{\lambda }|{\psi }_0⟩=0$. If the parts commute with each other, the Hamiltonian belongs to a specific subclass of frustration-free Hamiltonians known as "stabilizer Hamiltonians" . An example for a stabilizer Hamiltonian is Kitaev's toric code Hamiltonian ,

$H=-E_0(\sum _p{A}_p+\sum _v{B}_v),$

where the "plaquette" operators $A_p=\prod _{i\in p}{\sigma }_x^{(i)}$ and the "site" operators $B_s=\prod _{i\in s}{\sigma }_z^{(i)}$ contain Pauli matrices representing four-body spin interactions. The ground state of the toric code is given by the constraints

$\begin{array}{cc}⟨{\psi }_0|A_p|{\psi }_0⟩& =1\\ ⟨{\psi }_0|B_s|{\psi }_0⟩& =1\end{array}$

for all plaquettes and sites, respectively, which allows to solve the system exactly. For periodic boundary conditions on a torus the stabilizers satisfy the relations $\prod _p{A}_p=1$ and $\prod _s{B}_s=1$. for a system of $N$ atoms there are $N-2$ independent stabilizers. Consequently, the ground state manifold of the system will be four-fold degenerate. Excitations of the toric code Hamiltonian are violations of the stabilizer constraints and have energy gap of $2E_0$ as every violation will affect two sites or plaquettes, respectively. Interestingly, these quasi-particles are neither bosonic nor fermionic, but will pick up minus sign when moved around each other .

To show that it is possible to efficiently cool into the ground state, we consider a dynamical map of the form

$𝒱\rho =\sum _{\lambda }{p}_{\lambda }(P_{\lambda }\rho P_{\lambda }+\frac{1_{\lambda }}{\text{Tr}\left\{1_{\lambda }\right\}}\otimes {\text{Tr}}_{\lambda }\left\{H_{\lambda }\rho \right\}),$

where $p_{\lambda }$ are probabilities and $P_{\lambda }=1-H_{\lambda }$ are projectors onto the high-energy subspace with $H_{\lambda }\ne 0$ . This dynamics can be generated using local Hamiltonians and jump operators and any state in the ground state manifold, ${\rho }_0=|{\psi }_0⟩⟨{\psi }_0|$, is a stationary state of the dynamics satisfying $𝒱\rho =\rho$. It can also be shown that this dissipative preparation process requires only a time that is polynomially increasing with the size of the system .

Wave-function Monte-Carlo method

Finally, we want to conclude the course by taking a look at the practical problem of the numerical simulation of such dissipative quantum-many body systems. In general this is a very hard problem as the Hilbert space of a spin-1/2 many-body system grows exponentially like $2^N$ for $N$ spins. The size density matrix grows with the square of that, and including the cost for diagonalization of matrix to calculate the dynamics , the whole problem requires $O(2^{6N})$ steps, which becomes prohibitive for anything but very few spins. A way out of this difficulty is the reduction of the problem to work only on the level of wave functions and sample over many of them.

Consider an statistical ensemble of wave functions,

$\rho =\sum _{i=1}^M\frac{1}{M}|{\psi }_i⟩⟨{\psi }_i|.$

We now choose the size of the ensemble, $M$ to be much smaller than the Hilbert space dimension $2^N$. For each member of the ensemble, $|{\psi }_i⟩$, and at each timestep $\tau$, we perform the following procedure, i :

1. Evolve the state $|{\psi }_i$ under the non-Hermitian Hamiltonian

   $H^{\prime }=H-\frac{i\hslash }{2}\sum _i{\gamma }_i{c}_i^{†}{c}_i.$

   This can be done approximately by

   $|{\psi }_i(t+\tau )⟩=(1-i{H}^{\prime }\tau )|{\psi }_i(t)⟩+O({\tau }^2).$

2. For each jump operator $c_j$, calculate the probability for a quantum jump to occur, which is given by

   $p_j=\tau ⟨{\psi }_i(t)|c_j^{†}{c}_j|{\psi }_i(t)⟩$

3. Choose a uniformly distributed random number $r$ between 0 and 1.
4. Apply the jump operator to the quantum state, $c_j|{\psi }_i(t+\tau )⟩$, with $j$ chosen such that

   $\sum _{k=1}^{j-1}{p}_k<r<\sum _{k=1}^j{p}_k.$

   If $r$ is larger than the sum of all jump probabilities, no quantum jump occurs.

5. Normalize the resulting state.

If we now average observables over all states of the ensemble, we will reproduce the dynamics of the observables under the full master equation according to

$X(t)=\text{Tr}\{X\rho (t)\}=\frac{1}{M}\sum _i⟨{\psi }_i(t)|X|{\psi }_i(t)⟩+O\left(\frac{1}{\sqrt{M}}\right).$

Note that this method is only a first-order approximation to the dynamics and higher-order schemes exist . Alternatively, one can approximate the dynamics in the same way by a stochastic version of the Schrödinger equation .