Open Quantum Systems/Projection Operator Techniques

Decay into two bands

We consider a two-level systems coupled to an environment consisting of two distinct energy bands. The system is described by the Hamiltonian

$H=H_S+H_E+H_I.$

The system Hamiltonian $H_S$ is an ordinary two-level system with an energy difference of $\mathrm{\Delta }E$, according to

$H_S=\mathrm{\Delta }E|e⟩⟨e|=\mathrm{\Delta }E{\sigma }_{+}{\sigma }_{-},$

where we have set the energy of the ground state $|g⟩$ to zero. For the energy bands of the environment, we assume a linear spectrum, i.e,

$H_E=\sum _{n_1=1}^{N_1}\frac{\delta \epsilon }{N_1}{n}_1|n_1⟩⟨n_1|+\sum _{n_2=1}^{N_2}(\mathrm{\Delta }E+\frac{\delta \epsilon }{N_2}{n}_2)|n_2⟩⟨n_2|,$

where $N_1$ and $N_2$ are the total number of states in each band. The interaction Hamiltonian is given by

$H_I=\sum _{n_1=1}^{N_1}\sum _{n_2=1}^{N_2}{\lambda }_{n_1,n_2}({\sigma }_{+}\otimes 1_E)(1_S\otimes |n_1⟩⟨n_2|)+\text{H.c.},$

i.e., whenever the two-level system is de-excited, an excitation is created in the environment, and vice versa. The matrix elements ${\lambda }_{n1,n2}$ are random Gaussianly distributed variables with a variance ${\lambda }^2$. In the following, we assume a separation of energy scales according to $\mathrm{\Delta }E\gg \delta E\gg \lambda$.

This model has some very interesting properties. First of all, the size of the environment is clearly finite, meaning that excitations in the environment may not decay fast enough and the Markov approximation may not be applicable. This is even further strengthened when taking into account that the interaction will lead to very strong correlations between system and environment. At the same time, the weak interaction naturally gives rise to a peturbative expansion.

We proceed by defining a correlated projection operator for the relevant part of the dynamics , according to

$𝒫\rho =\text{Tr}\left\{{\mathrm{\Pi }}_1\rho \right\}\frac{1}{N_1}{\mathrm{\Pi }}_1+\text{Tr}\left\{{\mathrm{\Pi }}_2\rho \right\}\frac{1}{N_2}{\mathrm{\Pi }}_2=P_e\frac{1}{N_1}{\mathrm{\Pi }}_1+P_g\frac{1}{N_2}{\mathrm{\Pi }}_2,$

where we have introducted the projectors

$\begin{array}{cc}{\mathrm{\Pi }}_1& =\sum _{n_1=1}^{N_1}{1}_S\otimes |n_1⟩⟨n_1|\\ {\mathrm{\Pi }}_2& =\sum _{n_2=1}^{N_2}{1}_S\otimes |n_2⟩⟨n_2|\end{array}$

and $P_e$ and $P_g$ denote the probability to find the two-level system in its excited or ground state, respectively. Because of conservation of probability, we have $P_g=1-P_e$, meaning that the relevant part of the dynamics consists of a single variable!

In the following, we will derive its equation of motion. The second order TCL master equation is given by

$𝒫\frac{d}{dt}\rho =\underset{0}{\overset{t}{\int }}ds𝒫ℒ(t)ℒ(s)𝒫\rho ,$

which in our case can be written as

${\dot{P}}_e\frac{1}{N_1}{\mathrm{\Pi }}_1+(1-{\dot{P}}_e)\frac{1}{N_2}{\mathrm{\Pi }}_2=\underset{0}{\overset{t}{\int }}ds\text{Tr}\{{\mathrm{\Pi }}_1ℒ(t)ℒ(s)𝒫\rho \}\frac{1}{N_1}{\mathrm{\Pi }}_1+\text{Tr}\{{\mathrm{\Pi }}_2ℒ(t)ℒ(s)𝒫\rho \}\frac{1}{N_2}{\mathrm{\Pi }}_2.$

since the projection and the time derivative commute. Multiplying with ${\mathrm{\Pi }}_1$ and taking the trace yields

$\begin{array}{cc}{\dot{P}}_e& =\underset{0}{\overset{t}{\int }}ds\text{Tr}\{{\mathrm{\Pi }}_1ℒ(t)ℒ(s)𝒫\rho \}\\ & =-\underset{0}{\overset{t}{\int }}ds\text{Tr}\left\{{\mathrm{\Pi }}_1[H_I(t),[H_I(s),P_e\frac{1}{N_1}{\mathrm{\Pi }}_1+(1-P_e)\frac{1}{N_2}{\mathrm{\Pi }}_2]]\right\}.\end{array}$

Expanding the commutator and making use of the relation ${\mathrm{\Pi }}_i{\mathrm{\Pi }}_j={\delta }_{ij}{\mathrm{\Pi }}_i$ results in

${\dot{P}}_e=-\underset{0}{\overset{t}{\int }}ds\text{Tr}\{P_e\frac{1}{N_1}{\mathrm{\Pi }}_1{H}_I(s)H_I(t)-(1-P_e)\frac{1}{N_2}{\mathrm{\Pi }}_1{H}_I(t)H_I(s)+\text{H.c.}\}.$

We can now introduce two relaxation rates ${\gamma }_{1,2}$ describing the transition rates between the two bands. Specifically, we have

${\gamma }_i=\frac{1}{N_i}\underset{0}{\overset{t}{\int }}ds\sum _{n1}⟨n_1|H_I(s)H_I(t)+\text{H.c}|n_1⟩.$

As the interaction Hamiltonian couples does not contain terms that leave the band index unchanged, we may write the matrix elements as

$⟨n_1|H_I(s)H_I(t)+\text{H.c}|n_1⟩=\sum _{n_2}⟨n_1|H_I(s)|n_2⟩⟨n_2|H_I(t)|n_1⟩+\text{c.c.}.$

As we work in the interaction picture, we have

$⟨n_1|H_I(s)|n_2⟩⟨n_2|H_I(t)|n_1⟩=\exp (i{\omega }_{12}s)\exp (-i{\omega }_{12}t)|{\lambda }_{n1,n2}{|}^2,$

where we have introduced the transitions frequencies ${\omega }_{12}=\delta \epsilon (n_1/N_1-n_2/N_2)$. Using the statistical properties of the interaction Hamiltonian and performing the integration over $s$ finally yields

${\gamma }_i=\frac{2{\lambda }^2}{N_i}\sum _{n1,n2}\frac{\sin ({\omega }_{12}t)}{{\omega }_{12}}.$

For the next step, we note that the sinc function appearing in the relaxation rates is a representation of the Dirac $\delta$ distribution, i.e.,

$\pi {\delta }_t({\omega }_{12})=\underset{t\to \mathrm{\infty }}{\lim }\frac{\sin ({\omega }_{12}t)}{{\omega }_{12}}.$

Then, at large enough times and small enough couplings , we may approximate the rates by

$\begin{array}{cc}{\gamma }_1& \approx \frac{2\pi {\lambda }^2}{N_1}\sum _{k,l}\delta (E_{n_1}-E_{n_2})=2\pi {\lambda }^2\frac{N_2}{\delta \epsilon }\\ {\gamma }_2& \approx 2\pi {\lambda }^2\frac{N_1}{\delta \epsilon }.\end{array}$

The equation of motion for the excited state probability then reads

${\dot{P}}_e=-({\gamma }_2+{\gamma }_1)P_e+{\gamma }_2.$

Assuming the two-level system is initially in its excited state, the stationary state given by $P_e^s={\gamma }_2/({\gamma }_1+{\gamma }_2)$ is approached in an exponential decay according to

$P_e(t)=P_e^s+\frac{{\gamma }_1}{{\gamma }_1+{\gamma }_2}\exp [-2({\gamma }_1+{\gamma }_2)t].$

This solution is in excellent agreement with numerical simulations of the full Schrödinger equation .

It is important to stress that this efficient description of such a complex system is only possible because our projection operator correctly captures the essential features of the dynamics. If one choses a different projection operator, say, one without correlations between system and environments, according to

$𝒫\rho ={\text{Tr}}_E\{\rho \}\otimes {\rho }_E,$

then the second order TCL master equation no longer correctly describes the dynamics, as higher orders are divergent .