d^{2}\boldsymbol{\rho}_{1}\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho}_{1};\omega)\bar{\bar{G}}^{\dagger}(\boldsymbol{r}^{\prime},\boldsymbol{\rho}_{1};\omega)=$
$\displaystyle\frac{1}{4\bar{\varepsilon}^{2}}\int
d^{2}\boldsymbol{\rho}_{1}\int
d^{2}\boldsymbol{q}\frac{1}{q}\frac{e^{i\boldsymbol{q}\cdot(\boldsymbol{\rho}-\boldsymbol{\rho}_{1})}e^{-qz}}{i\omega\sigma(\omega)q/(2\bar{\varepsilon}\varepsilon_{0}c^{2})+\omega^{2}/c^{2}}\tilde{G}(\boldsymbol{q})$
$\displaystyle\int
d^{2}\boldsymbol{p}\frac{1}{p}\frac{e^{-i\boldsymbol{p}\cdot(\boldsymbol{\rho}^{\prime}-\boldsymbol{\rho}_{1})}e^{-pz^{\prime}}}{-i\omega\sigma^{*}(\omega)p/(2\bar{\varepsilon}\varepsilon_{0}c^{2})+\omega^{2}/c^{2}}\tilde{G}^{\dagger}(\boldsymbol{p})$

$=\frac{1}{4\bar{\varepsilon}^{2}}\int d^{2}\boldsymbol{q}\int
d^{2}\boldsymbol{p}\int
d^{2}\boldsymbol{\rho}_{1}\frac{1}{qp}\frac{e^{i\boldsymbol{q}\cdot\boldsymbol{\rho}+i\left(\boldsymbol{p}-\boldsymbol{q}\right)\cdot\boldsymbol{\rho}_{1}-i\boldsymbol{p}\cdot\boldsymbol{\rho}^{\prime}}e^{-qz}e^{-pz^{\prime}}}{D(q,\omega)D^{*}(p,\omega)}\tilde{G}(\boldsymbol{q})\tilde{G}^{\dagger}(\boldsymbol{p}),$
(433)

where $D(q,\omega)$ is shorthand for the denominator. Using the fact that

$\int
d^{2}\boldsymbol{\rho}_{1}e^{i\left(\boldsymbol{p}-\boldsymbol{q}\right)\cdot\boldsymbol{\rho}_{1}}=(2\pi)^{2}\delta(\boldsymbol{p}-\boldsymbol{q}),$
(434)

we have

$\int
d^{2}\boldsymbol{\rho}_{1}\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho}_{1};\omega)\bar{\bar{G}}^{\dagger}(\boldsymbol{r}^{\prime},\boldsymbol{\rho}_{1};\omega)=\frac{(2\pi)^{2}}{4\bar{\varepsilon}^{2}}\int
d^{2}\boldsymbol{q}\frac{1}{q^{2}}\frac{e^{i\boldsymbol{q}\cdot\left(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime}\right)}e^{-q(z+z^{\prime})}}{\left|i\omega\sigma(\omega)q/(2\bar{\varepsilon}\varepsilon_{0}c^{2})+\omega^{2}/c^{2}\right|^{2}}\tilde{G}(\boldsymbol{q})\tilde{G}^{\dagger}(\boldsymbol{q}).$
(435)

The commutator then reads

$\displaystyle\left[\hat{\boldsymbol{A}}(\boldsymbol{r},t),-\varepsilon_{0}\hat{\boldsymbol{E}}(\boldsymbol{r}^{\prime},t)\right]=$
$\displaystyle\frac{i\mu_{0}\pi\mathcal{C}}{c^{2}\bar{\varepsilon}^{2}}\int
d^{2}\boldsymbol{q}\frac{1}{q^{2}}e^{i\boldsymbol{q}\cdot\left(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime}\right)}e^{-q(z+z^{\prime})}\tilde{G}(\boldsymbol{q})\tilde{G}^{\dagger}(\boldsymbol{q})$
$\displaystyle\int_{0}^{\infty}\frac{d\omega}{2\pi}\frac{\omega\alpha(\omega)}{\left|i\omega\sigma(\omega)q/(2\bar{\varepsilon}\varepsilon_{0}c^{2})+\omega^{2}/c^{2}\right|^{2}}$
$\displaystyle=$
$\displaystyle\frac{i\pi\mathcal{C}}{\varepsilon_{0}\bar{\varepsilon}^{2}}\int
d^{2}\boldsymbol{q}\frac{1}{q^{2}}e^{i\boldsymbol{q}\cdot\left(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime}\right)}e^{-q(z+z^{\prime})}\tilde{G}(\boldsymbol{q})\tilde{G}^{\dagger}(\boldsymbol{q})$
$\displaystyle\int_{0}^{\infty}\frac{d\omega}{2\pi}\frac{\omega\alpha(\omega)}{\left|i\omega\sigma(\omega)q/(2\bar{\varepsilon}\varepsilon_{0})+\omega^{2}\right|^{2}}.$
(436)

Now we can move on to the evaluation of the integral over the frequency

$\mathcal{I}(q)=\int\frac{d\omega}{2\pi}\frac{\omega\alpha(\omega)}{\left|i\omega\sigma(\omega)q/(2\bar{\varepsilon}\varepsilon_{0})+\omega^{2}\right|^{2}},$
(437)

which when considering the complex conductivity
$\sigma(\omega)=\sigma_{1}(\omega)+i\sigma_{2}(\omega)$ as usual becomes

$\mathcal{I}(q)=\int_{0}^{\infty}\frac{d\omega}{2\pi}\frac{\omega\alpha(\omega)}{\left|i\omega\sigma_{1}(\omega)q/(2\bar{\varepsilon}\varepsilon_{0}c^{2})-\omega\sigma_{2}(\omega)q/(2\bar{\varepsilon}\varepsilon_{0}c^{2})+\omega^{2}/c^{2}\right|^{2}}.$
(438)

We can yet again evaluate this integral in the complex plane. In particular we
may look at an auxilliary problem where we consider a function $g(z)$ such
that its contour integral over a large semi-circle in the upper half-plane
vanishes, that is to say

$\oint dzg(z)=\lim_{R\to\infty}\left[\int_{-R}^{R}d\omega
g(\omega)+\int_{0}^{\pi}d\theta iRe^{i\theta}g(Re^{i\theta})\right]=0,$ (439)

which is equivalent to stating that $g(z)$ has no poles in the upper half
plane. We therefore know that the integrals are symmetric and as such, we may
evaluate one to find the other. In particular, we may evaluate pick a $g(z)$
of the following form

$g(z)=\frac{z}{z^{2}+iqz\sigma(z)/(2\bar{\varepsilon}\varepsilon_{0})},$ (440)

noting that if

$\lim_{R\to\infty}g(Re^{i\theta})=\frac{1}{Re^{i\theta}},$ (441)

we have that

$\int_{-\infty}^{\infty}d\omega g(\omega)=-i\pi,$ (442)

which folllows directly from equation (439), and noting that if we assume that
$iqz\sigma(z)/(2\bar{\varepsilon}\varepsilon_{0})$ grows more slowly than
$z^{2}$ (which it does, virtue of the dispersion relation we had found), then
equation (440) obeys the desired limit and thus can be evaluated as in
equation (442). Multiplying and dividing by the conjugate of the denominator
of $g(\omega),$we may write

$\displaystyle\int_{-\infty}^{\infty}d\omega\frac{\omega^{3}-iq\omega^{2}\sigma^{*}(\omega)/(2\bar{\varepsilon}\varepsilon_{0})}{|\omega^{2}+iq\omega\sigma(\omega)/(2\bar{\varepsilon}\varepsilon_{0})|^{2}}$
$\displaystyle=-i\pi$
$\displaystyle\Leftrightarrow\int_{-\infty}^{\infty}d\omega\frac{\omega^{3}-iq\omega^{2}\left[\sigma_{1}(\omega)-i\sigma_{2}(\omega)\right]/(2\bar{\varepsilon}\varepsilon_{0})}{|\omega^{2}+iq\omega\sigma(\omega)/(2\bar{\varepsilon}\varepsilon_{0})|^{2}}$
$\displaystyle=-i\pi,$ (443)

and if we consider the dispersion relation we found for $\sigma_{2}(\omega)$,
written

$\frac{\sigma_{2}(\omega)}{\varepsilon_{0}\omega}=\frac{2\bar{\varepsilon}}{q}$
(444)

we may write
$\sigma_{2}(\omega)=\frac{2\bar{\varepsilon}\varepsilon_{0}\omega}{q}$, and
therefore

$\displaystyle\int_{-\infty}^{\infty}d\omega\frac{\omega^{3}-iq\omega^{2}\left[\sigma_{1}(\omega)-i\sigma_{2}(\omega)\right]/(2\bar{\varepsilon}\varepsilon_{0})}{|\omega^{2}+iq\omega\sigma(\omega)/(2\bar{\varepsilon}\varepsilon_{0})|^{2}}$
$\displaystyle=\frac{-iq}{2\bar{\varepsilon}\varepsilon_{0}}\int_{-\infty}^{\infty}d\omega\frac{\sigma_{1}(\omega)\omega^{2}}{|\omega^{2}+iq\omega\sigma(\omega)/(2\bar{\varepsilon}\varepsilon_{0})|^{2}}$
$\displaystyle=-i\pi.$ (445)

We can now multiply by $-\frac{1}{2\pi i}$ on both sides, to obtain

$\int_{-\infty}^{\infty}d\omega\frac{1}{2\pi}\frac{\sigma_{1}(\omega)\omega^{2}}{|\omega^{2}+iq\omega\sigma(\omega)/(2\bar{\varepsilon}\varepsilon_{0})|^{2}}=\frac{\bar{\varepsilon}\varepsilon_{0}}{q},$

and so, we can set $\alpha(\omega)=\omega\sigma_{1}(\omega)$, and solve the
integral of equation (437) which (given that $\sigma_{1}(\omega)$ is an even
function) will evaluate to

$\mathcal{I}(q)=\frac{\bar{\varepsilon}\varepsilon_{0}}{2q}.$ (446)

Remembering that $\alpha(\omega)$ is chosen such that the commutator becomes a
longitudinal $\delta-$function, this choice of $\alpha(\omega)$ allows us to
write

$\displaystyle\left[\hat{\boldsymbol{A}}(\boldsymbol{r},t),-\varepsilon_{0}\hat{\boldsymbol{E}}(\boldsymbol{r}^{\prime},t)\right]$
$\displaystyle=\frac{i\pi\mathcal{C}}{\varepsilon_{0}\bar{\varepsilon}^{2}}\bar{\varepsilon}\varepsilon_{0}\int
d^{2}\boldsymbol{q}\frac{e^{i\boldsymbol{q}\cdot\left(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime}\right)}}{q^{3}}e^{-q(z+z^{\prime})}\tilde{G}(\boldsymbol{q})\tilde{G}^{\dagger}(\boldsymbol{q})$
$\displaystyle=\frac{i\pi\mathcal{C}}{\bar{\varepsilon}}\int
d^{2}\boldsymbol{q}\frac{e^{i\boldsymbol{q}\cdot\left(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime}\right)}}{q^{3}}e^{-q(z+z^{\prime})}\tilde{G}(\boldsymbol{q})\tilde{G}^{\dagger}(\boldsymbol{q}).$
(447)

Noting that using the general expression for
$\tilde{G}(\boldsymbol{q}),\tilde{G}^{\dagger}(\boldsymbol{q})$, one can
easily verify (in the electrostatic limit) that

$\tilde{G}(\boldsymbol{q})\tilde{G}^{\dagger}(\boldsymbol{q})=2q^{2}\tilde{G}(\boldsymbol{q}),$
(448)

which gives the commutator

$\left[\hat{\boldsymbol{A}}(\boldsymbol{r},t),-\varepsilon_{0}\hat{\boldsymbol{E}}(\boldsymbol{r}^{\prime},t)\right]=\frac{i\pi\mathcal{C}}{\bar{\varepsilon}}\int
d^{2}\boldsymbol{q}\frac{e^{i\boldsymbol{q}\cdot\left(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime}\right)}}{q}e^{-q(z+z^{\prime})}\tilde{G}(\boldsymbol{q}).$
(449)

We can now consider elements of this commutator, in particular, we may look at

$\left[\hat{A}_{i}(\boldsymbol{r},t),-\varepsilon_{0}\hat{E}_{j}(\boldsymbol{r}^{\prime},t)\right]=\frac{i\pi\mathcal{C}}{\bar{\varepsilon}}\int
d^{2}\boldsymbol{q}\frac{q_{i}q_{j}}{q}e^{i\boldsymbol{q}\cdot\left(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime}\right)}e^{-q(z+z^{\prime})},$
(450)

for the longitudinal components of the fields. $i=x,y$. Bringing the
exponentials together, by reintroducing
$\boldsymbol{k}\equiv\boldsymbol{p}_{1}=(q_{x},q_{y},i\kappa_{1})\approx(q_{x},q_{y},iq)$
and $\boldsymbol{r}=(x,y,z)$ and
$\boldsymbol{r}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})$, we can write

$\displaystyle\left[\hat{A}_{i}(\boldsymbol{r},t),-\varepsilon_{0}\hat{E}_{j}(\boldsymbol{r}^{\prime},t)\right]$
$\displaystyle=\frac{i\pi\mathcal{C}}{\bar{\varepsilon}}\int
d^{2}\boldsymbol{q}\frac{q_{i}q_{j}}{q}e^{i\boldsymbol{k}\cdot\left(\boldsymbol{r}-\boldsymbol{r}^{\prime}\right)},$
(451)

and extending the integral to the whole $q-$space as can easily be shown by
integrating in the complex plane

$\displaystyle\int
d^{3}\boldsymbol{k}\frac{k_{i}k_{j}}{q^{2}+k_{z}^{2}}e^{i\boldsymbol{k}\cdot\left(\boldsymbol{r}-\boldsymbol{r}^{\prime}\right)}$
$\displaystyle=\int
d^{2}kk_{i}k_{j}e^{iq(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})}\int
dk_{z}\frac{e^{-k_{z}(z+z^{\prime})}}{(q+ik_{z})(q-ik_{z})}$
$\displaystyle=\int
d^{2}kk_{i}k_{j}e^{iq(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})}2\pi
i\frac{e^{-q(z+z^{\prime})}}{2qi}$ $\displaystyle=\pi\int
d^{2}\boldsymbol{q}\frac{q_{i}q_{j}}{q}e^{iq(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})}.$
(452)

The commutator evaluates to

$\left[\hat{A}_{i}(\boldsymbol{r},t),-\varepsilon_{0}\hat{E}_{j}(\boldsymbol{r}^{\prime},t)\right]=\frac{i\mathcal{C}}{\bar{\varepsilon}}\int
d^{3}\boldsymbol{k}\frac{k_{i}k_{j}}{k^{2}}e^{i\boldsymbol{k}\cdot\left(\boldsymbol{r}-\boldsymbol{r}^{\prime}\right)}.$

This is very similar to the for of the longitudinal delta function
$\delta_{ij}^{\parallel}(\boldsymbol{r}-\boldsymbol{r}^{\prime})$. In fact,
recalling that

$\delta_{ij}^{\parallel}(\boldsymbol{r})\equiv\frac{1}{(2\pi)^{3}}\int
d^{3}\boldsymbol{k}\frac{k_{i}k_{j}}{k^{2}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}},$
(453)

we can recover equation (427) if we let

$\mathcal{C}=\frac{\bar{\varepsilon}\hbar}{(2\pi)^{3}}.$ (454)

Using the form of $\mathcal{C}$ and $\alpha(\omega)$ we discovered, the
commutators for the noise operators can also be found for $i$ and $j=x,y$, we
can write

$\left[\hat{J}_{N,i}(\boldsymbol{\rho};\omega),\hat{J}_{N,j}^{\dagger}(\boldsymbol{\rho}^{\prime};\omega^{\prime})\right]=\frac{\bar{\varepsilon}\hbar}{(2\pi)^{3}}\omega\sigma_{1}(\omega)\delta_{ij}\delta(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})\delta(\omega-\omega^{\prime}),$
(455)

and therefore we can define the operators

$\hat{\text{{f}}}(\boldsymbol{\rho},\omega)=\frac{\hat{\boldsymbol{J}}_{N}(\boldsymbol{\rho},\omega)}{\sqrt{\frac{\bar{\varepsilon}\hbar}{(2\pi)^{3}}\omega\sigma_{1}(\omega)}},$
(456)

which obey bosonic commutation relations

$\left[\hat{\text{f}}_{i}(\boldsymbol{\rho},\omega),\hat{\text{f}}_{i}^{\dagger}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\right]=\delta_{ij}\delta(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})\delta(\omega-\omega^{\prime}).$
(457)

The vector potential can thus be written making use of these operators as

$\displaystyle\hat{\boldsymbol{A}}(\boldsymbol{r},\omega)$
$\displaystyle=\mu_{0}\sqrt{\frac{\bar{\varepsilon}\hbar}{(2\pi)^{3}}\omega\sigma_{1}(\omega)}\int
d^{3}\boldsymbol{r}^{\prime}\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{r}^{\prime},\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\omega)\delta(z^{\prime})$
$\displaystyle\equiv\beta(\omega)\int
d^{2}\boldsymbol{\rho}^{\prime}\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho}^{\prime},z^{\prime}=0;\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\omega).$
(458)

The electric field operator can then be written at the cost of the vector
potential operator in the Schrodinger picture, by Fourier transforming at
$t=0$, giving the quantized field

$\hat{\boldsymbol{E}}(\boldsymbol{r})=i\int_{0}^{\infty}d\omega\omega\hat{\boldsymbol{A}}(\boldsymbol{r},\omega)-i\int_{0}^{\infty}d\omega\omega\hat{\boldsymbol{A}}^{\dagger}(\boldsymbol{r},\omega).$
(459)

This result concludes the quantization of the SPP field. In the next few
sections we will aim to use these results, especially the bosonic creation and
annihilation operators of equation (456), to derive the dynamics of two level-
systems coupled to graphene plasmons.

##  Appendix H Details of the calculation of the dynamics of one qubit coupled
to the plasmonic bath

###  H.1 Obtaining the dynamics via the Schrodinger equation

In the previous appendices we went over all the building blocks to set up a
mathematical description of our system, in particular the plasmonic bath. We
now move on to a quantitative description of the interaction between a two-
level system and the plasmonic bath. The Hamiltonian of a two-level system,
field and atom-field interaction is given, in the Schrodinger picture by

$\displaystyle\hat{H}=$ $\displaystyle\hbar\omega_{0}\sigma_{+}\sigma_{-}+\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\hbar\omega\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho},\omega)$
$\displaystyle-\sigma_{+}\boldsymbol{d}_{eg}\cdot\int_{0}^{\infty}d\omega\hat{\boldsymbol{E}}(\boldsymbol{r}_{A},\omega)-\sigma_{-}\boldsymbol{d}_{eg}\cdot\int_{0}^{\infty}d\omega\hat{\boldsymbol{E}}(\boldsymbol{r}_{A},\omega).$
(460)

In the rotating wave approximation, we keep only the resonant terms making up
the electric field. The Hamiltonian will therefore be

$\displaystyle\hat{H}=$ $\displaystyle\hbar\omega_{0}\sigma_{+}\sigma_{-}+\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\hbar\omega\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho},\omega)$
$\displaystyle-\sigma_{+}\boldsymbol{d}_{eg}\cdot
i\int_{0}^{\infty}d\omega\omega\hat{\boldsymbol{A}}(\boldsymbol{r}_{0},\omega)+i\int_{0}^{\infty}d\omega\omega\hat{\boldsymbol{A}}^{\dagger}(\boldsymbol{r}_{0},\omega)\cdot\boldsymbol{d}_{eg}\sigma_{-},$
(461)

or at the cost of only the creation and destruction operators

$\displaystyle\hat{H}=$ $\displaystyle\hbar\omega_{0}\sigma_{+}\sigma_{-}+\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\hbar\omega\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho},\omega)$
$\displaystyle-\sigma_{+}\boldsymbol{d}_{eg}\cdot
i\int_{0}^{\infty}d\omega\omega\beta(\omega)\int
d^{2}\boldsymbol{\rho}\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho}^{\prime},z^{\prime}=0;\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\omega)$
$\displaystyle+i\int_{0}^{\infty}d\omega\omega\beta^{*}(\omega)\int
d^{2}\boldsymbol{\rho}^{\prime}\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho}^{\prime},z^{\prime}=0;\omega)\cdot\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\omega)\cdot\boldsymbol{d}_{eg}\sigma_{-}.$
(462)

,We therefore have a zeroth order Hamiltonian of the form

$\hat{H}_{0}=\hbar\omega_{0}\sigma_{+}\sigma_{-}+\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\hbar\omega\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho},\omega),$
(463)

and an interaction term which can be written as

$\displaystyle\hat{V}=$ $\displaystyle-\sigma_{+}\boldsymbol{d}_{eg}\cdot
i\int_{0}^{\infty}d\omega\omega\beta(\omega)\int
d^{2}\boldsymbol{\rho}^{\prime}\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho}^{\prime},z^{\prime}=0;\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\omega)$
$\displaystyle+i\int_{0}^{\infty}d\omega\omega\beta(\omega)\int
d^{2}\boldsymbol{\rho}^{\prime}\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho}^{\prime},z^{\prime}=0;\omega)\cdot\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\omega)\cdot\boldsymbol{d}_{eg}\sigma_{-}.$
(464)

We can now aim to solve the Schrodinger equation

$i\hbar\frac{\partial}{\partial
t}\left|\psi(t)\right\rangle=H\left|\psi(t)\right\rangle.$ (465)

To solve the Schrodinger equation we make the following ansatz

$\left|\psi(t)\right\rangle=C_{e}(t)e^{-i\omega_{0}t}\left|0\right\rangle\left|e\right\rangle+\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega e^{-i\omega
t}\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t)\left|1(\boldsymbol{\rho},\omega,t)\right\rangle\left|g\right\rangle.$
(466)

Substituting onto the wave equation, on the left we have

$\displaystyle i\hbar\frac{\partial}{\partial
t}\left|\psi(t)\right\rangle=i\hbar\dot{C}_{e}(t)e^{-i\omega_{0}t}\left|0\right\rangle\left|e\right\rangle+\hbar\omega_{0}C_{e}(t)e^{-i\omega_{0}t}\left|0\right\rangle\left|e\right\rangle$
$\displaystyle+\int d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\hbar\omega
e^{-i\omega
t}\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t)\left|1(\boldsymbol{\rho},\omega)\right\rangle\left|g\right\rangle$
$\displaystyle+\int d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega i\hbar
e^{-i\omega
t}\dot{\boldsymbol{C}}_{g}(\boldsymbol{\rho},\omega,t)\left|1(\boldsymbol{\rho},\omega)\right\rangle\left|g\right\rangle.$
(467)

We now note that we have written $\left|\psi(t)\right\rangle$ as a function of
tensor product states made up of
$\left|0\right\rangle,\left|1(\boldsymbol{\rho},\omega)\right\rangle$ which
describe the plasmonic vacuum and the presence of a single plasmon in the
graphene and the states of the atom $\left|e\right\rangle$ and
$\left|g\right\rangle$, which are respectively the excited and ground state
for the two-level atom. This way, the operators which make up the Hamiltonian
act in the following ways

$\begin{cases}\sigma_{+}\left|i\right\rangle\left|e\right\rangle=0,\\\
\sigma_{+}\left|i\right\rangle\left|g\right\rangle=\left|i\right\rangle\left|e\right\rangle,\\\
\sigma_{-}\left|i\right\rangle\left|e\right\rangle=\left|i\right\rangle\left|g\right\rangle,\\\
\sigma_{-}\left|i\right\rangle\left|g\right\rangle=0,\end{cases}$ (468)
$\begin{cases}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\left|1(\boldsymbol{\rho},\omega)\right\rangle\left|j\right\rangle=\left|2(\boldsymbol{\rho},\omega)\right\rangle\left|j\right\rangle,\\\
\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\left|0\right\rangle\left|j\right\rangle=\left|1(\boldsymbol{\rho},\omega)\right\rangle\left|j\right\rangle,\\\
\hat{\text{{f}}}(\boldsymbol{\rho},\omega)\left|1(\boldsymbol{\rho},\omega)\right\rangle\left|j\right\rangle=\left|0\right\rangle\left|j\right\rangle,\\\
\hat{\text{{f}}}(\boldsymbol{\rho},\omega)\left|0\right\rangle\left|j\right\rangle=0,\end{cases},$
(469)

where $i=0,1$ and $j=g,e$. Note that since the operators are matched such that
$\hat{\text{{f}}}(\boldsymbol{\rho},\omega)$ is always paired with
$\sigma_{+}$ and vice-versa, the state
$\left|\boldsymbol{2}(\boldsymbol{\rho}_{1},\boldsymbol{\rho}_{2},\omega_{1},\omega_{2})\right\rangle\left|j\right\rangle$
is never actually produced. Acting with the zeroth order Hamiltonian on the
state therefore gives

$\hat{H}_{0}\left|\psi(t)\right\rangle=\hbar\omega_{0}C_{e}(t)e^{-i\omega_{0}t}+\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\hbar\omega\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t)e^{-i\omega
t}\left|1(\boldsymbol{\rho}^{\prime},\nu,t)\right\rangle\left|g\right\rangle,$
(470)

where we have also used that
$\hat{\text{{f}}}(\boldsymbol{\rho},\omega)\left|1(\boldsymbol{\rho}^{\prime},\nu,t)\right\rangle=\left|0\right\rangle\delta(\omega-\nu)\delta(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})$.
Using the same properties the interaction Hamiltonian gives

$\displaystyle\hat{V}\left|\psi(t)\right\rangle=$ $\displaystyle-i\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\omega\beta(\omega)\boldsymbol{d}_{eg}\cdot\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot
e^{-i\omega
t}\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t)\left|0\right\rangle\left|e\right\rangle$
$\displaystyle+i\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\omega\beta(\omega)\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot\boldsymbol{d}_{eg}C_{e}(t)e^{-i\omega_{0}t}\left|1(\boldsymbol{\rho},\omega,t)\right\rangle\left|g\right\rangle.$

Acting on the left with the bra $\left\langle e\right|\left\langle 0\right|$
on the left and hand side of the Schrodinger equation, we have

$\displaystyle i\hbar\dot{C}_{e}(t)e^{-i\omega_{0}t}=$ $\displaystyle-i\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\omega\beta(\omega)\boldsymbol{d}_{eg}\cdot\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)e^{-i\omega
t}\cdot\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t),$

or

$\dot{C}_{e}(t)=-\frac{1}{\hbar}\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\omega\beta(\omega)\boldsymbol{d}_{eg}\cdot\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)e^{-i(\omega-\omega_{0})t}\cdot\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t),$
(471)

while on the other hand, acting on the left with $\left\langle
g\right|\left\langle 1(\boldsymbol{\rho},\omega)\right|$ we have

$\displaystyle\hbar\omega e^{-i\omega
t}\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t)+i\hbar e^{-i\omega
t}\dot{\boldsymbol{C}}_{g}(\boldsymbol{\rho},\omega,t)=$
$\displaystyle\hbar\omega\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t)e^{-i\omega
t}$
$\displaystyle+i\omega\beta(\omega)\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot\boldsymbol{d}_{0}C_{e}(t)e^{-i\omega_{0}t}.$
(472)

$\Leftrightarrow\dot{\boldsymbol{C}}_{g}(\boldsymbol{\rho},\omega,t)=\frac{1}{\hbar}\omega\beta(\omega)\boldsymbol{d}_{eg}\cdot\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot
e^{-i(\omega_{0}-\omega)t}\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t).$
(473)

Integrating the previous equations with respect to time, choosing $C_{e}(0)=1$
and $\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,0)=\boldsymbol{0}$ we have

$C_{e}(t)=1-\frac{1}{\hbar}\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\omega\beta(\omega)\boldsymbol{d}_{eg}\cdot\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot\int_{0}^{t}dt_{1}e^{-i(\omega-\omega_{0})t_{1}}\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t_{1}),$
(474)

$\boldsymbol{C}_{g}(\boldsymbol{\rho},\omega,t)=\frac{1}{\hbar}\omega\beta(\omega)\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot\boldsymbol{d}_{eg}\int_{0}^{t}dt_{1}e^{-i(\omega_{0}-\omega)t_{1}}C_{e}(t_{1}).$
(475)

We can in particular, substitute equation (475) into equation (471), to obtain
the differential equation

$\displaystyle\dot{C}_{e}(t)=$ $\displaystyle-\frac{1}{\hbar}\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\omega\beta(\omega)\boldsymbol{d}_{eg}\cdot\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)e^{-i(\omega-\omega_{0})t}$
$\displaystyle\cdot\frac{1}{\hbar}\omega\beta(\omega)\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot\boldsymbol{d}_{eg}\int_{0}^{t}dt_{1}e^{-i(\omega_{0}-\omega)t_{1}}C_{e}(t_{1}).$
(476)

We can simplify the previous expression as

$\displaystyle\dot{C}_{e}(t)=$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{\infty}d\omega\omega^{2}\beta^{2}(\omega)\boldsymbol{d}_{eg}\cdot\int
d^{2}\boldsymbol{\rho}\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)$
$\displaystyle\cdot\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot\boldsymbol{d}_{eg}e^{-i(\omega-\omega_{0})t}\int_{0}^{t}dt_{1}e^{-i(\omega_{0}-\omega)t_{1}}C_{e}(t_{1}),$
(477)

and we can further introduce the kernel (see appendix F.11)
$\bar{\mathcal{K}}(\boldsymbol{r},\omega)$ defined by

$\displaystyle\bar{\mathcal{K}}(\boldsymbol{r},\omega)$
$\displaystyle=\varepsilon_{0}\mu_{0}\omega^{2}\frac{\sigma_{1}(\omega)}{\varepsilon_{0}\omega}\int
d^{2}\boldsymbol{\rho}\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)$
$\displaystyle=\frac{\omega^{2}}{c^{2}}\frac{\sigma_{1}(\omega)}{\varepsilon_{0}\omega}\int
d^{2}\boldsymbol{\rho}\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)$
$\displaystyle=\text{Im}[\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{r};\omega)],$
(478)

We shall discuss the evaluation of this integral in due time, but first we
move on to write

$\dot{C}_{e}(t)=-\frac{\bar{\varepsilon}\mu_{0}}{(2\pi)^{3}\hbar}\int_{0}^{\infty}d\omega\omega^{2}\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}\int_{0}^{t}dt_{1}e^{-i(\omega_{0}-\omega)(t_{1}-t)}C_{e}(t_{1}).$
(479)

We can now make the first Markov approximation to write $C_{e}(t_{1})\approx
C_{e}(t)$. We then have

$\displaystyle\dot{C}_{e}(t)$
$\displaystyle=-\frac{\bar{\varepsilon}\mu_{0}}{(2\pi)^{3}\hbar}\int_{0}^{\infty}d\omega\omega^{2}\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}C_{e}(t)\int_{0}^{t}dt_{1}e^{-i(\omega_{0}-\omega)(t_{1}-t)}$
$\displaystyle=-\frac{\bar{\varepsilon}\mu_{0}}{(2\pi)^{3}\hbar}\int_{0}^{\infty}d\omega\omega^{2}\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}C_{e}(t).$
(480)

By the Sokhotski–Plemelj identity the integral is

$\mathcal{I}_{SP}(\omega)=\pi\delta(\omega-\omega_{0})+i\mathcal{P}\frac{1}{\omega-\omega_{0}}$
(481)

Using this result, we have

$\dot{C_{e}}(t)=iC_{e}(t)\left[-\Delta+i\Gamma\right],$ (482)

where we have introduced the Lamb Shift

$\Delta=\frac{\bar{\varepsilon}\mu_{0}}{(2\pi)^{3}\hbar}\mathcal{P}\int
d\omega\omega^{2}\frac{\boldsymbol{d}_{eg}\cdot\text{Im}[\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{r};\omega)]\cdot\boldsymbol{d}_{eg}}{\omega-\omega_{0}},$
(483)

and the decay rate

$\displaystyle\Gamma$
$\displaystyle=-\pi\frac{\bar{\varepsilon}\mu_{0}}{(2\pi)^{3}\hbar}\int
d\omega\omega^{2}\boldsymbol{d}_{eg}\cdot\text{Im}[\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{r};\omega)]\cdot\boldsymbol{d}_{eg}\delta(\omega-\omega_{0})$
$\displaystyle=-\pi\frac{\bar{\varepsilon}}{(2\pi)^{3}\hbar\varepsilon_{0}c^{2}}\omega_{0}^{2}\boldsymbol{d}_{eg}\cdot\text{Im}[\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{r};\omega_{0})]\cdot\boldsymbol{d}_{eg}.$
(484)

We will now compute this rate using the dispersion relation

####  H.1.1 Computation of the decay rate of a single qubit

We now aim to calculate the decay rate for our two-level system close to the
Graphene sheet. We therefore will make use of the Green’s function we found
previously in equation (35), to calculate

$\text{Im}[\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{r};\omega_{0})]=c^{2}\frac{\pi^{2}}{\bar{\varepsilon}}e^{-2q_{\text{spp}}z}\frac{q_{\text{spp}}^{2}}{A(\omega_{0})}\text{Im}[\bar{M}(q_{\text{ssp}},0,0)],$
(485)

where $A$ is in general, a function of $\omega_{0}$. The imaginary part of the
evaluated matrix gives

$\begin{cases}\text{Im}\left[M_{11}\right]=-J_{0}\left(q_{\text{spp}}R\right)+\frac{J_{1}\left(q_{\text{spp}}R\right)}{q_{\text{spp}}R},\\\
\text{Im}\left[M_{22}\right]=-\frac{J_{1}\left(q_{\text{spp}}R\right)}{q_{\text{spp}}R},\\\
\text{Im}\left[M_{33}\right]=-J_{0}\left(q_{\text{spp}}R\right),\end{cases}$
(486)

and

$\begin{cases}\text{Im}\left[M_{12}\right]=\text{Im}\left[M_{21}\right]=0,\\\
\text{Im}\left[M_{13}\right]=\text{Im}\left[M_{31}\right]=J_{1}\left(q_{\text{spp}}R\right),\\\
\text{Im}\left[M_{23}\right]=\text{Im}\left[M_{32}\right]=0.\end{cases}.$
(487)

Therefore, we have, if we let $a\equiv q_{\text{spp}}R$ for the sake of
simplifying the notation, we find

$\Gamma=-\pi\frac{\bar{\varepsilon}}{(2\pi)^{3}\hbar\varepsilon_{0}c^{2}}\omega_{0}^{2}\boldsymbol{d}_{eg}\cdot\text{Im}[\bar{G}_{1}(\boldsymbol{r},\boldsymbol{r};\omega_{0})]\cdot\boldsymbol{d}_{eg}$

$\displaystyle\Gamma$
$\displaystyle=-\frac{\pi\bar{\varepsilon}\mu_{0}}{(2\pi)^{3}\hbar}\omega_{0}^{2}\left(-c^{2}\frac{\pi^{2}}{\bar{\varepsilon}}e^{-2q_{\text{spp}}z}\frac{q_{\text{spp}}^{2}}{A(\omega)}\right)$
$\displaystyle\times\boldsymbol{d}_{eg}\cdot\left[\lim_{a\to
0}\begin{pmatrix}J_{0}\left(a\right)-\frac{J_{1}\left(a\right)}{a}&0&-J_{1}\left(a\right)\\\
0&\frac{J_{1}\left(a\right)}{a}&0\\\
-J_{1}\left(a\right)&0&J_{0}\left(a\right)\end{pmatrix}\right]\cdot\boldsymbol{d}_{eg}.$
(488)

Introducing the notation

$\displaystyle\Gamma$
$\displaystyle=c^{2}\frac{\mu_{0}}{8\hbar}\omega_{0}^{2}e^{-2q_{\text{spp}}z}\frac{q_{\text{spp}}^{2}}{A(\omega_{0})}\boldsymbol{d}_{eg}\cdot\begin{pmatrix}1/2&0&0\\\
0&1/2&0\\\ 0&0&1\end{pmatrix}\cdot\boldsymbol{d}_{eg}$
$\displaystyle=\frac{1}{16\hbar\varepsilon_{0}}\frac{\omega_{0}^{2}}{A(\omega_{0})}e^{-2q_{\text{spp}}z}q_{\text{spp}}^{2}\boldsymbol{d}_{eg}\cdot\begin{pmatrix}1&0&0\\\
0&1&0\\\ 0&0&2\end{pmatrix}\cdot\boldsymbol{d}_{eg}$
$\displaystyle=\frac{1}{16\hbar\varepsilon_{0}}q_{\text{spp}}^{3}e^{-2q_{\text{spp}}z}(d_{eg}^{\parallel
2}+2d_{eg}^{\perp 2}),$ (489)

and writing

$q_{\text{spp}}^{3}=\frac{(2\pi)^{3}}{\lambda_{\text{spp}}^{3}},$ (490)

we have a result compatible with the literature

$\Gamma=\frac{1}{16\hbar\varepsilon_{0}}\frac{(2\pi)^{3}}{\lambda_{\text{spp}}^{3}}e^{-2q_{\text{spp}}z}(d_{eg}^{\parallel
2}+2d_{eg}^{\perp 2})$ (491)

or by writing $\Gamma=\gamma/2$, we find the rate

$\gamma=\frac{1}{8\hbar\varepsilon_{0}}\frac{(2\pi)^{3}}{\lambda_{\text{spp}}^{3}}e^{-2q_{\text{spp}}z}(d_{eg}^{\parallel
2}+2d_{eg}^{\perp 2}).$ (492)

This is the final result for the rate $\gamma$ that characterizes the
spontaneous decay of a two-level system coupled to an SPP field.

###  H.2 Deriving the Lindblad equation and computing the dynamics of a single
qubit

We now aim to present a microscopic derivation of the Lindblad equation for
our system. This will allow us to extend these calculations to the case where
both qubits are excited, but we first present a derivation for a single qubit
coupled to the SPP field. As with the case of the thermal EM field, we start
with noting that for interaction Hamiltonians of the form

$\hat{V}=\sum_{i}\beta_{i}\otimes\Gamma_{i},$ (493)

it is possible to write the Redfield equation in terms of the reservoir
correlation functions as

$\displaystyle\frac{d\hat{\rho}_{S}}{dt}=$
$\displaystyle-\frac{1}{\hbar^{2}}\sum_{ij}\int_{t_{0}}^{t}d\tau\left(\beta_{i}(t)\beta_{j}(\tau)\hat{\rho}_{S}(t)-\beta_{j}(\tau)\hat{\rho}_{S}(t)\beta_{i}(t)\right)\left\langle\Gamma_{i}(t)\Gamma_{j}(\tau)\right\rangle_{R}$
$\displaystyle+\left(\hat{\rho}_{S}(t)\beta_{j}(\tau)\beta_{i}(t)-\beta_{i}(t)\hat{\rho}_{S}(t)\beta_{j}(\tau)\right)\left\langle\Gamma_{j}(\tau)\Gamma_{i}(t)\right\rangle_{R}.$
(494)

These reservoir correlation functions are defined as

$\begin{cases}\left\langle\Gamma_{i}(t)\Gamma_{j}(\tau)\right\rangle\equiv\text{tr}_{R}\left\\{\hat{\rho}_{R}\Gamma_{i}(t)\Gamma_{j}(\tau)\right\\},\\\
\left\langle\Gamma_{j}(t)\Gamma_{i}(\tau)\right\rangle\equiv\text{tr}_{R}\left\\{\hat{\rho}_{R}\Gamma_{j}(t)\Gamma_{i}(\tau)\right\\}.\end{cases}$
(495)

The Hamiltonian for a single qubit are given by equation (603) if the indices
$\alpha$ are all set to $1$. Since the commutation relations for SPP creation
and destruction operators follow directly from equation (469), we can find the
operators in the interaction picture via

$\displaystyle\hat{\text{{f}}}_{I}^{\dagger}(\boldsymbol{\rho},\omega)$
$\displaystyle=e^{i\hat{H}_{0}t/\hbar}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)e^{-i\hat{H}_{0}t/\hbar}.$

We make use of the commutator property $Ae^{B}=e^{c}e^{B}A$, where
$\left[A,B\right]=cA$. In this case, we have

$e^{i\hat{H}_{0}t/\hbar}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)e^{-i\hat{H}_{0}t/\hbar}=e^{c}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega),$
(496)

where

$\displaystyle\left[\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega),-i\hat{H}_{0}t/\hbar\right]$
$\displaystyle=\left[\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega),-it/\hbar\int
d^{2}\boldsymbol{\rho}^{\prime}\int_{0}^{\infty}d\nu\hbar\nu\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\nu)\right]$
$\displaystyle=-it\int
d^{2}\boldsymbol{\rho}^{\prime}\int_{0}^{\infty}d\nu\nu\left[\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega),\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\nu)\right].$
(497)

The commutator can be expanded as

$\displaystyle\left[\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega),\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\nu)\right]$
$\displaystyle=\left[\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega),\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\right]\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\nu)+\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\left[\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega),\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\nu)\right]$
$\displaystyle=-\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\delta(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})\delta(\omega-\nu),$
(498)

from which we have

$\displaystyle\left[\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega),-i\hat{H}_{0}t/\hbar\right]$
$\displaystyle=it\int
d^{2}\boldsymbol{\rho}^{\prime}\int_{0}^{\infty}d\nu\nu\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\delta(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})\delta(\omega-\nu)$
$\displaystyle=i\omega t\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega).$
(499)

Thus, we can use the privously mentioned identity, with $c=i\omega t$. The
result is the expected

$\hat{\text{{f}}}_{I}^{\dagger}(\boldsymbol{\rho},\omega)=\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)e^{i\omega
t}.$ (500)

A similar calculation yields

$\left[\hat{\text{{f}}}(\boldsymbol{\rho},\omega),-i\hat{H}_{0}t/\hbar\right]=\int
d^{2}\boldsymbol{\rho}^{\prime}\int_{0}^{\infty}d\nu\nu\left[\hat{\text{{f}}}(\boldsymbol{\rho},\omega),\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\nu)\right],$
(501)

where the commutator yields

$\displaystyle\left[\hat{\text{{f}}}(\boldsymbol{\rho},\omega),\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\nu)\right]$
$\displaystyle=\left[\hat{\text{{f}}}(\boldsymbol{\rho},\omega),\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\right]\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\nu)+\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\nu)\left[\hat{\text{{f}}}(\boldsymbol{\rho},\omega),\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\nu)\right]$
$\displaystyle=\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\nu)\delta(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})\delta(\omega-\nu),$
(502)

from where follows

$\left[\hat{\text{{f}}}(\boldsymbol{\rho},\omega),-i\hat{H}_{0}t/\hbar\right]=-i\omega
t\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega),$ (503)

which gives the representation

$\hat{\text{{f}}}_{I}(\boldsymbol{\rho},\omega)=\hat{\text{{f}}}(\boldsymbol{\rho},\omega)e^{-i\omega
t}.$ (504)

As we have seen before, we also have for the qubit raising and lowering
operators

$\begin{cases}\sigma_{+}^{I}=\sigma_{+}e^{i\omega_{0}t},\\\
\sigma_{-}^{I}=\sigma_{-}e^{-i\omega_{0}t},\end{cases}$ (505)

We can now identify the Hamiltonian of SPP field coupled to the qubit as one
of the form (493), by making the definitions

$\begin{cases}\hbar\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\equiv\omega\beta(\omega)\boldsymbol{d}\cdot\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega),\\\
\hbar\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega)\equiv\omega\beta^{*}(\omega)\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot\boldsymbol{d}^{\dagger},\end{cases}$

and then

$\begin{cases}\beta_{1}=\sigma_{-},\\\ \beta_{2}=\sigma_{+},\\\
\Gamma_{1}=\Gamma^{\dagger}=i\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\cdot\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega),\\\
\Gamma_{2}=\Gamma=-i\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho},\omega),\end{cases}$
(506)

with which the interaction Hamiltonian takes a form identical to that we had
previously obtained when treating the thermal E-M field

$\hat{V}=\hbar\left(\sigma_{-}\Gamma^{\dagger}+\sigma_{+}\Gamma\right).$

The procedure to derive the Lindblad equation is therefore very similar to the
case of the thermal E-M field. We can write all the operators in the
interaction picture

$\begin{cases}\tilde{\beta}_{1}(t)=\sigma_{-}e^{-i\omega_{0}t},\\\
\tilde{\beta}_{2}(t)=\sigma_{+}e^{i\omega_{0}t},\\\
\tilde{\Gamma}_{1}(t)=i\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\cdot\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega)e^{i\omega
t},\\\ \tilde{\Gamma}_{2}(t)=-i\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho},\omega)e^{-i\omega
t}.\end{cases}$ (507)

The summation present in the Redfield equation is over $i=1,2$ and $j=1,2$ and
therefore, we can evaluate all the reservoir correlation functions in the
interaction picture

$\displaystyle\left\langle\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle$
$\displaystyle=\text{tr}_{R}\left\\{\hat{\rho}_{R}\tilde{\Gamma}(t)\tilde{\Gamma}(\tau)\right\\}$
$\displaystyle=\text{tr}_{R}\left\\{-\hat{\rho}_{R}\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\cdot\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega)e^{i\omega
t}\right.$ $\displaystyle\left.\int_{0}^{\infty}d\omega^{\prime}\int
d^{2}\boldsymbol{\rho}^{\prime}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\cdot\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho}^{\prime};\omega^{\prime})e^{i\omega^{\prime}\tau}\right\\}$
$\displaystyle=0,$ (508)

in the rotating wave approximation. Similarly

$\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{2}(\tau)\right\rangle=0.$
(509)

The story becomes much more complicated, however, for the remaining reservoir
correlation functions
$\left\langle\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{2}(\tau)\right\rangle$ and
$\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle$. We
will therefore consider first the case of zero temperature, which is much
simpler and will allow us to recover the results obtained in the Schrodinger
picture. Then, we will move on to a description at finite temperature.

###  H.3 Evaluating the correlation functions at temperature $T=0$: Density
matrix $\rho_{R}=|0\rangle\langle 0|$

In this appendix we provide the aforementioned zero temperature description of
the The reason for choosing the temperature $T=0$ is that the density matrix
of the reservoir will be given by

$\rho_{R}=|\psi\rangle\langle\psi|=|0\rangle\langle 0|,$ (510)

since it has not been perturbed by the decay of the qubit. We can also work in
the subspace of the reservoir with basis
$\left\\{\left|0\right\rangle,\left|\mathds{1}(\boldsymbol{\rho},\omega)\right\rangle\right\\}$,
since at most there will be a single SPP present. In this subspace, the trace
will in principle be calculated over two terms, however we see that one of
these terms immediately vanishes due to the form of the density matrix. This
results in

$\displaystyle\left\langle\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{2}(\tau)\right\rangle=$
$\displaystyle\text{tr}_{R}\left\\{\hat{\rho}_{R}\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{2}(\tau)\right\\}$
$\displaystyle=\braket{0}{0}\left\langle 0\right|\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\cdot\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega)e^{i\omega
t}$ $\displaystyle\int_{0}^{\infty}d\omega^{\prime}\int
d^{2}\boldsymbol{\rho}^{\prime}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho}^{\prime};\omega^{\prime})\cdot\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\omega^{\prime})e^{i\omega^{\prime}\tau}\left|0\right\rangle$
$\displaystyle=0.$

The other correlation function is given by

$\displaystyle\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle$
$\displaystyle=\text{tr}_{R}\left\\{\hat{\rho}_{R}\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\\}$
$\displaystyle=\braket{0}{0}\left\langle 0\right|\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho},\omega)e^{i\omega
t}\int_{0}^{\infty}d\omega^{\prime}\int
d^{2}\boldsymbol{\rho}^{\prime}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\cdot\mathcal{G}^{+}(\boldsymbol{r}^{\prime},\boldsymbol{\rho}^{\prime};\omega^{\prime})e^{i\omega^{\prime}\tau}\left|0\right\rangle$
$\displaystyle=\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\cdot\left\langle\mathds{1}(\boldsymbol{\rho},\omega)\right|e^{-i\omega
t}\int_{0}^{\infty}d\omega^{\prime}\int
d^{2}\boldsymbol{\rho}^{\prime}\left|\mathds{1}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\right\rangle\cdot\mathcal{G}^{+}(\boldsymbol{r}^{\prime},\boldsymbol{\rho}^{\prime};\omega^{\prime})e^{i\omega^{\prime}\tau}$
$\displaystyle=\int_{0}^{\infty}d\omega^{\prime}\int
d^{2}\boldsymbol{\rho}^{\prime}\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\cdot\mathcal{G}^{+}(\boldsymbol{r}^{\prime},\boldsymbol{\rho}^{\prime};\omega^{\prime})\delta(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})\delta(\omega-\omega^{\prime})e^{i\omega^{\prime}\tau}e^{-i\omega
t}$ $\displaystyle=\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\cdot\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega)e^{i\omega(\tau-t)}.$
(511)

Similarly

$\left\langle\tilde{\Gamma}_{2}(\tau)\tilde{\Gamma}_{1}(t)\right\rangle=\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle^{*}.$

Substituting into the Redfield equation (539), all terms but those
proportional to this correlation function vanish, and therefore we are left
with

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{t_{0}}^{t}d\tau\left(\tilde{\beta}_{2}(t)\tilde{\beta}_{1}(\tau)\hat{\rho}_{S}^{I}(t)-\tilde{\beta}_{1}(\tau)\hat{\rho}_{S}^{I}(t)\tilde{\beta}_{2}(t)\right)\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle_{R}$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\tilde{\beta}_{2}(\tau)\tilde{\beta}_{1}(t)-\tilde{\beta}_{1}(t)\hat{\rho}_{S}^{I}(t)\tilde{\beta}_{2}(\tau)\right)\left\langle\tilde{\Gamma}_{2}(\tau)\tilde{\Gamma}_{1}(t)\right\rangle_{R},$
(512)

and substituting the operators, we obtain

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{t}d\tau\left(\sigma_{+}\sigma_{-}e^{i\omega_{0}(t-\tau)}\hat{\rho}_{S}^{I}(t)-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}e^{i\omega_{0}(t-\tau)}\right)\left\langle\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{2}(\tau)\right\rangle$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}e^{i\omega_{0}(\tau-t)}-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}e^{i\omega_{0}(\tau-t)}\right)\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle^{*}.$
(513)

If we introduce the notation

$\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}\equiv\frac{1}{\hbar^{2}}\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\cdot\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega).$
(514)

We can perform the Markov approximation, extending the integral in time to
infinity and writing $\hat{\rho}_{S}^{I}(t_{1})\approx\hat{\rho}_{S}^{I}(t)$.
This yields

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\left(\sigma_{+}\sigma_{-}\hat{\rho}_{S}^{I}(t)-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)\int_{0}^{\infty}d\omega\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}\int_{0}^{\infty}dt_{1}e^{i(\omega-\omega_{0})(t_{1}-t)}$
$\displaystyle-\left(\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)\int_{0}^{\infty}d\omega\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}\int_{0}^{\infty}dt_{1}e^{-i(\omega-\omega_{0})(t_{1}-t)},$
(515)

and we obtain

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\left(\sigma_{+}\sigma_{-}\hat{\rho}_{S}^{I}(t)-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)\int_{0}^{\infty}d\omega\bar{\mathcal{K}}(\boldsymbol{r},\omega)\mathcal{I}(\omega)$
$\displaystyle-\left(\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)\int_{0}^{\infty}d\omega\bar{\mathcal{K}}(\boldsymbol{r},\omega)\mathcal{I}^{*}(\omega),$
(516)

where

$\mathcal{I}(\omega)=\pi\delta(\omega-\omega_{\alpha^{\prime}})-i\mathcal{P}\frac{1}{\omega-\omega_{\alpha^{\prime}}}.$
(517)

We can therefore define

$\int_{0}^{\infty}d\omega\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}\mathcal{I}(\omega)\equiv\Gamma+i\Delta,$
(518)

which lets us write

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\left(\sigma_{+}\sigma_{-}\hat{\rho}_{S}^{I}(t)-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)(\Gamma+i\Delta)$
$\displaystyle-\left(\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)(\Gamma-i\Delta).$
(519)

We can write this more compactly as

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\left(\sigma_{+}\sigma_{-}\hat{\rho}_{S}^{I}(t)-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}+\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)\Gamma$
$\displaystyle-\left(\sigma_{+}\sigma_{-}\hat{\rho}_{S}^{I}(t)-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}+\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)i\Delta$
$\displaystyle=$
$\displaystyle-\Gamma\left(\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right\\}-2\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)-i\Delta\left(\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right]\right).$

We now want to return to the Shcrodinger picture, and therefore use the
definition of the operator in the interaction picture

$\hat{\rho}_{S}^{I}(t)=e^{iH_{0}t/\hbar}\hat{\rho}_{S}e^{-iH_{0}t/\hbar},$
(520)

with which we can calculate the derivative

$\frac{d\hat{\rho}_{S}^{I}}{dt}=e^{iH_{0}t/\hbar}\frac{d\hat{\rho}_{S}}{dt}e^{-iH_{0}t/\hbar}+\frac{i}{\hbar}e^{iH_{0}t/\hbar}\left[H_{0},\hat{\rho}_{S}(t)\right]e^{-iH_{0}t/\hbar},$
(521)

and can further rewrite

$\frac{d\hat{\rho}_{S}}{dt}=e^{-iH_{0}t/\hbar}\frac{d\hat{\rho}_{S}^{I}}{dt}e^{iH_{0}t/\hbar}-\frac{i}{\hbar}\left[H_{0},\hat{\rho}_{S}(t)\right].$
(522)

We can now convert from the interaction picture to the Schrodinger picture, by
writing

$\displaystyle\frac{d\hat{\rho}_{S}}{dt}=$
$\displaystyle-\Gamma\left(e^{-iH_{0}t/\hbar}\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right\\}e^{iH_{0}t/\hbar}-e^{-iH_{0}t/\hbar}2\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}e^{iH_{0}t/\hbar}\right)$
$\displaystyle-i\Delta\left(e^{-iH_{0}t/\hbar}\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right]e^{iH_{0}t/\hbar}\right)$
$\displaystyle-\frac{i}{\hbar}\left[H_{0},\hat{\rho}_{S}(t)\right]$
$\displaystyle=\Gamma\left(2\sigma_{-}\hat{\rho}_{S}(t)\sigma_{+}-\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right\\}\right)-i\Delta\left(\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right]\right)-\frac{i}{\hbar}\left[H_{0},\hat{\rho}_{S}(t)\right].$
(523)

We can now define

$\begin{cases}\hat{L}=\sqrt{\Gamma}\sigma_{-},\\\
\hat{L}^{\dagger}=\sqrt{\Gamma}\sigma_{+},\end{cases}$ (524)

which yields

$\displaystyle\frac{d\hat{\rho}_{S}}{dt}$
$\displaystyle=2\hat{L}\hat{\rho}_{S}(t)\hat{L}^{\dagger}-\left\\{\hat{L}^{\dagger}\hat{L},\hat{\rho}_{S}(t)\right\\}-\frac{i}{\hbar}\left[(\omega+\Delta)\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right].$
(525)

This is an equation in Lindblad form, and we now check that the dynamics
reproduces those found via Schrodinger’s equation. We write the previous
equation in the convenient form

$\frac{d\hat{\rho}_{S}}{dt}=\Gamma\left(2\sigma_{-}\hat{\rho}_{S}(t)\sigma_{+}-\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right\\}\right)-\frac{i}{\hbar}\left[\hbar(\omega_{0}+\Delta)\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right].$

We now expand the density matrix as

$\displaystyle\hat{\rho}_{S}$
$\displaystyle=\sum_{ij}\rho_{ij}(t)\left|i\right\rangle\left\langle
j\right|,$ (526)

but look only at the populations $P_{i}=\rho_{ii}$. To do so, we substitute
this onto the differential equation above and act on the left and right with
the bra and ket corresponding to the state $\left|e\right\rangle$. This yields

$\displaystyle\dot{P}_{e}(t)=$ $\displaystyle\Gamma\left(\left\langle
e\right|2\sigma_{-}\hat{\rho}_{S}(t)\sigma_{+}\left|e\right\rangle-\left\langle
e\right|\sigma_{+}\sigma_{-}\hat{\rho}_{S}(t)\left|e\right\rangle-\left\langle
e\right|\hat{\rho}_{S}(t)\sigma_{+}\sigma_{-}\left|e\right\rangle\right)$
$\displaystyle-\frac{i}{\hbar}(\omega_{0}+\Delta)\left(\left\langle
e\right|\sigma_{+}\sigma_{-}\hat{\rho}_{S}(t)\left|e\right\rangle-\left\langle
e\right|\hat{\rho}_{S}(t)\sigma_{+}\sigma_{-}\left|e\right\rangle\right)$
$\displaystyle=$
$\displaystyle\Gamma\left(-P_{e}-P_{e}\right)-\frac{i}{\hbar}(\omega_{0}+\Delta)(P_{e}-P_{e})$
$\displaystyle=$ $\displaystyle-2\Gamma P_{e}(t).$ (527)

While for the ground state

$\displaystyle\dot{P}_{g}(t)=$ $\displaystyle\Gamma\left(\left\langle
g\right|2\sigma_{-}\hat{\rho}_{S}(t)\sigma_{+}\left|g\right\rangle-\left\langle
g\right|\sigma_{+}\sigma_{-}\hat{\rho}_{S}(t)\left|g\right\rangle-\left\langle
g\right|\hat{\rho}_{S}(t)\sigma_{+}\sigma_{-}\left|g\right\rangle\right)$
$\displaystyle-\frac{i}{\hbar}(\omega_{0}+\Delta)\left(\left\langle
g\right|\sigma_{+}\sigma_{-}\hat{\rho}_{S}(t)\left|g\right\rangle-\left\langle
g\right|\hat{\rho}_{S}(t)\sigma_{+}\sigma_{-}\left|g\right\rangle\right)$
$\displaystyle=$ $\displaystyle 2\Gamma P_{e}(t).$ (528)

This yields the same dynamics as we had obtained via Schrodinger’s equation,
namely

$P_{e}(t)=P_{e}(0)e^{-2\Gamma t}.$ (529)

and for the ground state

$\displaystyle P_{g}(t)$ $\displaystyle=2\Gamma\int_{0}^{t}P_{e}(t^{\prime})$
$\displaystyle=2\Gamma P_{e}(0)\frac{e^{-2\Gamma t}-1}{-2\Gamma}$
$\displaystyle=P_{e}(0)\left(1-e^{-2\Gamma t}\right)$ (530)

Setting $\gamma=2\Gamma$ we obtain

$P_{g}(t)=P_{e}(0)\left(1-e^{-\gamma t}\right).$ (531)

Note that $\gamma$ has exactly of the same form as derived in appendix H.1.1.

###  H.4 Evaluating the correlation functions at a finite temperature $T$

The Lindblad approach has the advantage of being able to derive these results
for a finite temperature aswell, where the SPP field will be described by a
bosonic distribution

$\bar{n}_{\text{spp}}(\omega)=\frac{1}{e^{\hbar\omega/k_{B}T}-1}.$ (532)

We will use the following results found in Scheel and Buhmann [9].

$\begin{cases}\left\langle\hat{\text{{f}}}(\boldsymbol{\rho},\omega)\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\right\rangle_{R}=0,\\\
\left\langle\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\right\rangle_{R}=0,\\\
\left\langle\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\right\rangle_{R}=\bar{n}_{\text{spp}}(\omega)\delta(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})\delta(\omega-\omega^{\prime}),\\\
\left\langle\hat{\text{{f}}}(\boldsymbol{\rho},\omega)\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\right\rangle_{R}=\left[\bar{n}_{\text{spp}}(\omega)+1\right]\delta(\boldsymbol{\rho}-\boldsymbol{\rho}^{\prime})\delta(\omega-\omega^{\prime}).\end{cases}.$
(533)

The derivation of the Lindblad equation proceeds exactly as described in
section B.2, and as before we have only to derive the reservoir correlation
functions

$\displaystyle\left\langle\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle$
$\displaystyle=\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho}^{\prime};\omega^{\prime})\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega)e^{-i\omega^{\prime}\tau}e^{i\omega
t}\left\langle\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\right\rangle_{R}$
$\displaystyle=0.$ (534)

Similarly, using
$\left\langle\hat{\text{{f}}}(\boldsymbol{\rho},\omega)\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\right\rangle=0$,
we find

$\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{2}(\tau)\right\rangle=0.$
(535)

Therefore, we have to calculate only

$\displaystyle\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle$
$\displaystyle=\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega^{\prime}\int
d^{2}\boldsymbol{\rho}^{\prime}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho}^{\prime};\omega^{\prime})e^{-i\omega
t}e^{i\omega^{\prime}\tau}\left\langle\hat{\text{{f}}}(\boldsymbol{\rho},\omega)\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\right\rangle_{R}$
$\displaystyle=\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega)\left[\bar{n}_{\text{spp}}(\omega)+1\right]e^{i\omega(\tau-t)}.$
(536)

In the same manner, we can find

$\displaystyle\left\langle\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{2}(\tau)\right\rangle$
$\displaystyle=\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega^{\prime}\int
d^{2}\boldsymbol{\rho}^{\prime}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho}^{\prime};\omega^{\prime})e^{i\omega
t}e^{-i\omega^{\prime}\tau}\left\langle\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\hat{\text{{f}}}(\boldsymbol{\rho}^{\prime},\omega^{\prime})\right\rangle_{R}$
$\displaystyle=\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega)\bar{n}_{\text{spp}}(\omega)e^{-i\omega(\tau-t)}.$
(537)

We also have

$\begin{cases}\left\langle\tilde{\Gamma}_{2}(\tau)\tilde{\Gamma}_{1}(t)\right\rangle_{R}=\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle_{R}^{*},\\\
\left\langle\tilde{\Gamma}_{1}(\tau)\tilde{\Gamma}_{2}(t)\right\rangle_{R}=\left\langle\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{2}(\tau)\right\rangle_{R}^{*}.\end{cases}$
(538)

$\displaystyle\frac{d\hat{\rho}_{S}}{dt}=$
$\displaystyle-\frac{1}{\hbar^{2}}\sum_{ij}\int_{t_{0}}^{t}d\tau\left(\beta_{i}(t)\beta_{j}(\tau)\hat{\rho}_{S}(t)-\beta_{j}(\tau)\hat{\rho}_{S}(t)\beta_{i}(t)\right)\left\langle\Gamma_{i}(t)\Gamma_{j}(\tau)\right\rangle_{R}$
$\displaystyle+\left(\hat{\rho}_{S}(t)\beta_{j}(\tau)\beta_{i}(t)-\beta_{i}(t)\hat{\rho}_{S}(t)\beta_{j}(\tau)\right)\left\langle\Gamma_{j}(\tau)\Gamma_{i}(t)\right\rangle_{R}.$
(539)

Therefore, the Redfield equation reads,

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{t}d\tau\left(\tilde{\beta}_{2}(t)\tilde{\beta}_{1}(\tau)\hat{\rho}_{S}^{I}(t)-\tilde{\beta}_{1}(\tau)\hat{\rho}_{S}^{I}(t)\tilde{\beta}_{2}(t)\right)\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle_{R}$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\tilde{\beta}_{2}(\tau)\tilde{\beta}_{1}(t)-\tilde{\beta}_{1}(t)\hat{\rho}_{S}^{I}(t)\tilde{\beta}_{2}(\tau)\right)\left\langle\tilde{\Gamma}_{2}(\tau)\tilde{\Gamma}_{1}(t)\right\rangle_{R}$
$\displaystyle+\left(\tilde{\beta}_{1}(t)\tilde{\beta}_{2}(\tau)\hat{\rho}_{S}^{I}(t)-\tilde{\beta}_{2}(\tau)\hat{\rho}_{S}^{I}(t)\tilde{\beta}_{1}(t)\right)\left\langle\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{2}(\tau)\right\rangle_{R}$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\tilde{\beta}_{1}(\tau)\tilde{\beta}_{2}(t)-\tilde{\beta}_{2}(t)\hat{\rho}_{S}^{I}(t)\tilde{\beta}_{1}(\tau)\right)\left\langle\tilde{\Gamma}_{1}(\tau)\tilde{\Gamma}_{2}(t)\right\rangle_{R},$
(540)

and using the previous definitions for the operators $\beta_{i}$ and
$\Gamma_{i}$ we write it as

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{t}d\tau\left(\sigma_{+}\sigma_{-}\hat{\rho}_{S}^{I}(t)-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)e^{-i\omega_{0}(\tau-t)}\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle_{R}$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)e^{i\omega_{0}(\tau-t)}\left\langle\tilde{\Gamma}_{2}(t)\tilde{\Gamma}_{1}(\tau)\right\rangle_{R}^{*}$
$\displaystyle+\left(\sigma_{-}\sigma_{+}\hat{\rho}_{S}^{I}(t)-\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}\right)e^{i\omega_{0}(\tau-t)}\left\langle\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{2}(\tau)\right\rangle_{R}$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{-}\sigma_{+}-\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}\right)e^{-i\omega_{0}(\tau-t)}\left\langle\tilde{\Gamma}_{1}(t)\tilde{\Gamma}_{2}(\tau)\right\rangle_{R}^{*}.$
(541)

Substituting the previously evaluated reservoir correlation functions, and
using the notation

$\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}\equiv\frac{1}{\hbar^{2}}\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\mathcal{G}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega),$
(542)

we write

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\int_{0}^{\infty}d\omega\left(\sigma_{+}\sigma_{-}\hat{\rho}_{S}^{I}(t)-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)\bar{\mathcal{K}}(\boldsymbol{r},\omega)\left[\bar{n}_{\text{spp}}(\omega)+1\right]\int_{0}^{t}d\tau
e^{i(\omega-\omega_{0})(\tau-t)}$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}\right)\bar{\mathcal{K}}(\boldsymbol{r},\omega)\left[\bar{n}_{\text{spp}}(\omega)+1\right]\int_{0}^{t}d\tau
e^{-i(\omega-\omega_{0})(\tau-t)}$
$\displaystyle+\left(\sigma_{-}\sigma_{+}\hat{\rho}_{S}^{I}(t)-\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}\right)\bar{\mathcal{K}}(\boldsymbol{r},\omega)\bar{n}_{\text{spp}}(\omega)\int_{0}^{t}d\tau
e^{-i(\omega-\omega_{0})(\tau-t)}$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{-}\sigma_{+}-\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}\right)\bar{\mathcal{K}}(\boldsymbol{r},\omega)\bar{n}_{\text{spp}}(\omega)\int_{0}^{t}d\tau
e^{i(\omega-\omega_{0})(\tau-t)}.$ (543)

Letting the integral go to infinity and using the Sokhotski-Plemelj, we can
define

$\begin{cases}\Gamma\equiv\pi\int_{0}^{\infty}d\omega\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}\delta(\omega-\omega_{0}),\\\
\Delta\equiv\mathcal{P}\int_{0}^{\infty}d\omega\frac{\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}}{\omega_{0}-\omega},\\\
\Delta^{\prime}\equiv\mathcal{P}\int_{0}^{\infty}d\omega\frac{\boldsymbol{d}_{eg}\cdot\bar{\mathcal{K}}(\boldsymbol{r},\omega)\cdot\boldsymbol{d}_{eg}\bar{n}_{\text{spp}}(\omega)}{\omega_{0}-\omega}.\end{cases}$
(544)

to obtain

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle\left(\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\sigma_{+}\sigma_{-}\hat{\rho}_{S}^{I}(t)\right)\left[\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]+i(\Delta+\Delta^{\prime})\right]$
$\displaystyle+\left(\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}\right)\left[\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]-i(\Delta+\Delta^{\prime})\right]$
$\displaystyle\text{+}\left(\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}-\sigma_{-}\sigma_{+}\hat{\rho}_{S}^{I}(t)\right)\left[\Gamma\bar{n}_{\text{spp}}(\omega_{0})-i\Delta^{\prime}\right]$
$\displaystyle+\left(\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}-\hat{\rho}_{S}^{I}(t)\sigma_{-}\sigma_{+}\right)\left[\Gamma\bar{n}_{\text{spp}}(\omega_{0})+i\Delta^{\prime}\right].$
(545)

Grouping together all the terms in , we obtain

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\sigma_{+}\sigma_{-}\hat{\rho}_{S}^{I}(t)+\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}\right)$
$\displaystyle+i(\Delta+\Delta^{\prime})\left(\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\sigma_{+}\sigma_{-}\hat{\rho}_{S}^{I}(t)-\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}+\hat{\rho}_{S}^{I}(t)\sigma_{+}\sigma_{-}\right)$
$\displaystyle+\Gamma\bar{n}_{\text{spp}}(\omega_{0})\left(\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}-\sigma_{-}\sigma_{+}\hat{\rho}_{S}^{I}(t)+\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}-\hat{\rho}_{S}^{I}(t)\sigma_{-}\sigma_{+}\right)$
$\displaystyle-i\Delta^{\prime}\left(\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}-\sigma_{-}\sigma_{+}\hat{\rho}_{S}^{I}(t)-\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}+\hat{\rho}_{S}^{I}(t)\sigma_{-}\sigma_{+}\right),$
(546)

which can be compacted as

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(2\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle-i(\Delta+\Delta^{\prime})\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right]$
$\displaystyle+\Gamma\bar{n}_{\text{spp}}(\omega_{0})\left(2\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}-\left\\{\sigma_{-}\sigma_{+},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle+i\Delta^{\prime}\left[\sigma_{-}\sigma_{+},\hat{\rho}_{S}^{I}(t)\right].$
(547)

We now note that the $\sigma_{-}\sigma_{+}$ can be written as

$\sigma_{-}\sigma_{+}=\sigma_{+}\sigma_{-}+\left[\sigma_{-},\sigma_{+}\right],$
(548)

where this commutator can be calculated employing

$\begin{cases}\sigma_{-}=\frac{1}{2}(\sigma_{x}-i\sigma_{y}),\\\
\sigma_{+}=\frac{1}{2}(\sigma_{x}+i\sigma_{y}),\end{cases}$ (549)

and the commutator relations of the Pauli matrices

$\left[\sigma_{i},\sigma_{j}\right]=2i\varepsilon_{ijk}\sigma_{k}.$ (550)

This yields

$\displaystyle\left[\sigma_{-},\sigma_{+}\right]$
$\displaystyle=\frac{1}{4}\left[\sigma_{x}-i\sigma_{y},\sigma_{x}+i\sigma_{y}\right]$
$\displaystyle=\frac{1}{4}\left(\left[\sigma_{x},\sigma_{x}\right]+i\left[\sigma_{x},\sigma_{y}\right]-i\left[\sigma_{y},\sigma_{x}\right]+\left[\sigma_{y},\sigma_{y}\right]\right)$
$\displaystyle=\frac{i}{2}\left[\sigma_{x},\sigma_{y}\right]$
$\displaystyle=-\sigma_{z}.$ (551)

Substituting this into our dynamical equation, we obtain

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(2\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle-\Gamma\bar{n}_{\text{spp}}(\omega_{0})\left(2\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}-\left\\{\sigma_{-}\sigma_{+},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle-i(\Delta+\Delta^{\prime})\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right]$
$\displaystyle+i\Delta^{\prime}\left[\sigma_{+}\sigma_{-}-\sigma_{z},\hat{\rho}_{S}^{I}(t)\right].$
(552)

We can also use the fact that

$\sigma_{+}\sigma_{-}=\frac{1}{2}\left(\mathds{1}+\sigma_{z}\right),$ (553)

or rather, that

$\sigma_{z}=2\sigma_{+}\sigma_{-}-\mathds{1}.$ (554)

This yields

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(2\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle-\Gamma\bar{n}_{\text{spp}}(\omega_{0})\left(2\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}-\left\\{\sigma_{-}\sigma_{+},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle-i(\Delta+\Delta^{\prime})\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right]$
$\displaystyle+i\Delta^{\prime}\left[\sigma_{+}\sigma_{-}-2\sigma_{+}\sigma_{-}+\mathds{1},\hat{\rho}_{S}^{I}(t)\right],$
(555)

and since $\left[\mathds{1},\hat{\rho}_{S}^{I}(t)\right]=0$, we obtain

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(2\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle+\Gamma\bar{n}_{\text{spp}}(\omega_{0})\left(2\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}-\left\\{\sigma_{-}\sigma_{+},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle-i(\Delta+\Delta^{\prime})\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right]$
$\displaystyle-i\Delta^{\prime}\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right],$
(556)

which finally yields

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(2\sigma_{-}\hat{\rho}_{S}^{I}(t)\sigma_{+}-\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle+\Gamma\bar{n}_{\text{spp}}(\omega_{0})\left(2\sigma_{+}\hat{\rho}_{S}^{I}(t)\sigma_{-}-\left\\{\sigma_{-}\sigma_{+},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle-i(\Delta+2\Delta^{\prime})\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}^{I}(t)\right].$
(557)

We can convert this into the Schrodinger picture, and the procedure is similar
to that carried out in detail for a single qubit. In the end, we simply get
rid of the $I$ index and add a term corresponding to the coherent evolution of
the system. This procedure results in

$\displaystyle\frac{d\hat{\rho}_{S}}{dt}=$
$\displaystyle\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(2\sigma_{-}\hat{\rho}_{S}(t)\sigma_{+}-\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right\\}\right)$
$\displaystyle+\Gamma\bar{n}_{\text{spp}}(\omega_{0})\left(2\sigma_{+}\hat{\rho}_{S}(t)\sigma_{-}-\left\\{\sigma_{-}\sigma_{+},\hat{\rho}_{S}(t)\right\\}\right)$
$\displaystyle-i(\Delta+2\Delta^{\prime})\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right]-\frac{i}{\hbar}\left[H_{0},\hat{\rho}_{S}(t)\right].$
(558)

Since $H_{0}=\hbar\omega_{0}\sigma_{+}\sigma_{-}$, we can write compactly

$\displaystyle\frac{d\hat{\rho}_{S}}{dt}=$
$\displaystyle\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(2\sigma_{-}\hat{\rho}_{S}(t)\sigma_{+}-\left\\{\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right\\}\right)$
$\displaystyle+\Gamma\bar{n}_{\text{spp}}(\omega_{0})\left(2\sigma_{+}\hat{\rho}_{S}(t)\sigma_{-}-\left\\{\sigma_{-}\sigma_{+},\hat{\rho}_{S}(t)\right\\}\right)$
$\displaystyle-i(\omega_{0}+\Delta+2\Delta^{\prime})\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right].$
(559)

Defining the operators

$\begin{cases}\hat{L}_{1}=\sqrt{\Gamma\bar{n}_{\text{spp}}(\omega_{0})}\sigma_{\text{+}},\\\
\hat{L}_{2}=\sqrt{\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]}\sigma_{-},\end{cases}$
(560)

we can write the equation in Lindblad form

$\displaystyle\frac{d\hat{\rho}_{S}}{dt}=$
$\displaystyle\left(2\hat{L}_{2}\hat{\rho}_{S}(t)\hat{L}_{2}^{\dagger}-\left\\{\hat{L}_{2}^{\dagger}\hat{L}_{2},\hat{\rho}_{S}(t)\right\\}\right)$
$\displaystyle+\left(2\hat{L}_{1}\hat{\rho}_{S}(t)\hat{L}_{1}^{\dagger}-\left\\{\hat{L}_{1}^{\dagger}\hat{L}_{1},\hat{\rho}_{S}(t)\right\\}\right)$
$\displaystyle-i(\omega_{0}+\Delta+2\Delta^{\prime})\left[\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right],$
(561)

or introducing a sum over $k=1,2$, the result

$\frac{d\hat{\rho}_{S}}{dt}=-\frac{i}{\hbar}\left[\hbar\omega_{0}^{\prime}\sigma_{+}\sigma_{-},\hat{\rho}_{S}(t)\right]+\sum_{k=1,2}\left(2\hat{L}_{k}\hat{\rho}_{S}(t)\hat{L}_{k}^{\dagger}-\left\\{\hat{L}_{k}^{\dagger}\hat{L}_{k},\hat{\rho}_{S}(t)\right\\}\right),$
(562)

with $\omega_{0}^{\prime}=\omega_{0}+\Delta+\Delta^{\prime}$. We can convert
this into a dynamical equation for the probability of finding each state by
proceeding analogously to the case of a single qubit. We write

$\displaystyle\hat{\rho}_{S}$
$\displaystyle=\rho_{ee}\left|e\right\rangle\left\langle
e\right|+\rho_{eg}\left|e\right\rangle\left\langle
g\right|+\rho_{ge}\left|g\right\rangle\left\langle
e\right|+\rho_{gg}\left|g\right\rangle\left\langle g\right|$
$\displaystyle=\begin{bmatrix}\rho_{ee}&\rho_{eg}\\\
\rho_{ge}&\rho_{gg}\end{bmatrix}.$ (563)

We can write the Lindblad jump operators in matrix form

$\begin{cases}\hat{L}_{1}&=\sqrt{\Gamma\bar{n}_{\text{spp}}(\omega_{0})}\begin{bmatrix}0&1\\\
0&0\end{bmatrix},\\\
\hat{L}_{2}&=\sqrt{\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]}\begin{bmatrix}0&0\\\
1&0\end{bmatrix}.\end{cases}$ (564)

We can now substitute these objects into the Lindblad equation and setting
$\gamma=2\Gamma$, we find

$\begin{bmatrix}\dot{\rho}_{ee}&\dot{\rho}_{eg}\\\
\dot{\rho}_{ge}&\dot{\rho}_{gg}\end{bmatrix}=-i4\omega_{0}\begin{bmatrix}0&\rho_{eg}\\\
\rho_{ge}&0\end{bmatrix}+\begin{bmatrix}-\gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{ee}+\gamma\bar{n}_{\text{spp}}(\omega_{0})\rho_{gg}&-\frac{\gamma}{2}\left(\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]+\bar{n}_{\text{spp}}(\omega_{0})\right)\rho_{eg}\\\
-\frac{\gamma}{2}\left(\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]+\bar{n}_{\text{spp}}(\omega_{0})\right)\rho_{ge}&\gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{ee}-\gamma\bar{n}_{\text{spp}}(\omega_{0})\rho_{gg}\end{bmatrix}.$
(565)

Substituting into equation (559) and looking at the diagonal elements which
correspond to the populations, we obtain a coupled set of differential
equations for the probabilities

$\begin{cases}\dot{P}_{e}=-\gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]P_{e}+\gamma\bar{n}_{\text{spp}}(\omega_{0})P_{g},\\\
\dot{P}_{g}=\gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]P_{e}-\gamma\bar{n}_{\text{spp}}(\omega_{0})P_{g}.\end{cases}$
(566)

This is slightly more complicated then the case at zero temperature since the
qubits may also be excited by thermal fluctuations of the SPP field. These
differential equations have a general solution given by

$\displaystyle P_{e}(t)=$ $\displaystyle
C_{1}\frac{\gamma\left(\bar{n}_{\text{spp}}(\omega_{0})+\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]e^{-\gamma(2\bar{n}_{\text{spp}}(\omega_{0})+1)t}\right)}{\gamma\left[2\bar{n}_{\text{spp}}(\omega_{0})+1\right]}$
$\displaystyle+C_{2}\frac{\gamma\bar{n}_{\text{spp}}(\omega_{0})\left(-1+e^{-\gamma(2\bar{n}_{\text{spp}}(\omega_{0})+1)t}\right)}{\gamma\left[2\bar{n}_{\text{spp}}(\omega_{0})+1\right]},$
(567)

$\displaystyle P_{g}(t)=$ $\displaystyle-
C_{1}\frac{\gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(-1+e^{-\gamma(2\bar{n}_{\text{spp}}(\omega_{0})+1)t}\right)}{2\Gamma\left[2\bar{n}_{\text{spp}}(\omega_{0})+1\right]}$
$\displaystyle+C_{2}\frac{\gamma\left(\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]+\bar{n}_{\text{spp}}(\omega_{0})e^{-\gamma(2\bar{n}_{\text{spp}}(\omega_{0})+1)t}\right)}{\gamma\left[2\bar{n}_{\text{spp}}(\omega_{0})+1\right]}.$
(568)

If we choose initial conditions such that $P_{g}(0)=0$ and $P_{e}(0)=1$, we
can simplifiy the previous equations since we find

$C_{2}\frac{\left(\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]+\bar{n}_{\text{spp}}(\omega_{0})\right)}{\left[2\bar{n}_{\text{spp}}(\omega_{0})+1\right]}=0,$
(569)

which gives $C_{2}=0$. Therefore, we can write

$P_{e}(0)=C_{1}\frac{2\bar{n}_{\text{spp}}(\omega_{0})+1}{2\bar{n}_{\text{spp}}(\omega_{0})+1}=1,$
(570)

which gives $C_{1}=1$. We therefore have

$P_{e}(t)=\frac{\left(\bar{n}_{\text{spp}}(\omega_{0})+\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]e^{-\gamma(2\bar{n}_{\text{spp}}(\omega_{0})+1)t}\right)}{2\bar{n}_{\text{spp}}(\omega_{0})+1}.$
(571)

We can even attempt to further simplify this, using

$\bar{n}_{\text{spp}}(\omega)=\frac{1}{e^{\hbar\omega/k_{B}T}-1},$ (572)

which yields

$\displaystyle\bar{n}_{\text{spp}}(\omega)+1$
$\displaystyle=\frac{1}{e^{\hbar\omega/k_{B}T}-1}+1$
$\displaystyle=\frac{e^{\hbar\omega/k_{B}T}}{e^{\hbar\omega/k_{B}T}-1},$

and

$\displaystyle 2\bar{n}_{\text{spp}}(\omega)+1$
$\displaystyle=\frac{2}{e^{\hbar\omega/k_{B}T}-1}+1$
$\displaystyle=\frac{e^{\hbar\omega/k_{B}T}+1}{e^{\hbar\omega/k_{B}T}-1}.$
(573)

We can see that for $\hbar\omega_{0}\gg k_{B}T$, both terms approach $1$,
while $\bar{n}_{\text{spp}}(\omega)\to 0$. In this limit we recover

$P_{e}(t)\approx e^{-\gamma t},$ (574)

which is the result for zero temperature. The results thus appear to be
consistent. In the other limit for large temperature, we find
$\bar{n}_{\text{spp}}(\omega)+1\approx\bar{n}_{\text{spp}}(\omega)$. And thus,
we find

$\displaystyle P_{e}(t)$
$\displaystyle=\frac{\bar{n}_{\text{spp}}(\omega_{0})\left(1+e^{-2\gamma\bar{n}_{\text{spp}}(\omega_{0})t}\right)}{2\bar{n}_{\text{spp}}(\omega_{0})}$
$\displaystyle=\frac{1}{2}\left(1+e^{-2\gamma\bar{n}_{\text{spp}}(\omega_{0})t}\right),$
(575)

while

$\displaystyle P_{g}(t)$
$\displaystyle=-\frac{\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(-1+e^{-\gamma(2\bar{n}_{\text{spp}}(\omega_{0})+1)t}\right)}{\left[2\bar{n}_{\text{spp}}(\omega_{0})+1\right]}$
$\displaystyle\approx-\frac{1}{2}\left(-1+e^{-2\gamma\bar{n}_{\text{spp}}(\omega_{0})t}\right)$
$\displaystyle\approx\frac{1}{2}\left(1-e^{-2\gamma\bar{n}_{\text{spp}}(\omega_{0})t}\right),$
(576)

and both states will tend to equilibrium at equal probability.

##  Appendix I Thermal averages of $\hat{\boldsymbol{\text{f}}}$-operators

In this appendix we present a small derivation of the thermal averages of the
surface-plasmon-polariton creation and anihilation operators and their
product, based on discretizing the integrals that they are made up of. Let us
consider the Hamiltonian of the free electromagnetic field described by the
operators
$\hat{\mathbf{f}}(\boldsymbol{r},\omega)=\left\\{f_{x}(\mathbf{r},\omega),f_{y}(\mathbf{r},\omega),f_{z}(\mathbf{r},\omega)\right\\}$

$\displaystyle H_{F}$ $\displaystyle=\int d^{3}\boldsymbol{r}\int
d\omega\hbar\omega\hat{\mathbf{f}}^{\dagger}(\boldsymbol{r},\omega)\cdot\hat{\mathbf{f}}(\boldsymbol{r},\omega)$
$\displaystyle=\sum_{i}\int d^{3}\boldsymbol{r}\int d\omega\hbar\omega
f_{i}^{\dagger}(\boldsymbol{r},\omega)f_{i}(\boldsymbol{r},\omega).$ (577)

Let us discretize the integrals over $\mathbf{r}$ and $\omega$ to write the
Hamiltonian as

$H_{F}=\sum_{i}\sum_{n}\Delta_{\boldsymbol{r}}\sum_{\ell}\Delta_{\omega}\hbar\omega_{\ell}f_{i}^{\dagger}(\boldsymbol{r}_{n},\omega_{\ell})f_{i}(\boldsymbol{r}_{n},\omega_{\ell}).$
(578)

Now let us define the operators

$\displaystyle\tilde{f}_{i,n,\ell}$
$\displaystyle=\sqrt{\Delta_{\boldsymbol{r}}\Delta_{\omega}}f_{i}(\boldsymbol{r}_{n},\omega_{\ell}),$
(579) $\displaystyle\tilde{f}_{i,n,\ell}^{\dagger}$
$\displaystyle=\sqrt{\Delta_{\boldsymbol{r}}\Delta_{\omega}}f_{i}^{\dagger}(\boldsymbol{r}_{n},\omega_{\ell}).$
(580)

Let us assume that

$\left[\tilde{f}_{i,n,\ell},\tilde{f}_{i^{\prime},n^{\prime},\ell^{\prime}}^{\dagger}\right]=\delta_{i,i^{\prime}}\delta_{n,n^{\prime}}\delta_{\ell,\ell^{\prime}}$
(581)

This leads us to

$\displaystyle\left[f_{i}(\boldsymbol{r}_{n},\omega_{\ell}),f_{i^{\prime}}^{\dagger}(\boldsymbol{r}_{n^{\prime}},\omega_{\ell^{\prime}})\right]$
$\displaystyle=\frac{1}{\Delta_{\boldsymbol{r}}\Delta_{\omega}}\left[\tilde{f}_{i,n,\ell},\tilde{f}_{i^{\prime},n^{\prime},\ell^{\prime}}^{\dagger}\right]$
$\displaystyle=\frac{1}{\Delta_{\boldsymbol{r}}\Delta_{\omega}}\delta_{i,i^{\prime}}\delta_{n,n^{\prime}}\delta_{\ell,\ell^{\prime}}$
$\displaystyle=\delta_{i,i^{\prime}}\frac{1}{\Delta_{\boldsymbol{r}}}\delta_{n,n^{\prime}}\frac{1}{\Delta_{\omega}}\delta_{\ell,\ell^{\prime}}$
(582)

In the limit of $\Delta_{\mathbf{r}}\rightarrow 0$ and
$\Delta_{\omega}\rightarrow 0$, we obtain

$\displaystyle\lim_{\Delta_{\mathbf{r}}\rightarrow
0}\frac{1}{\Delta_{\mathbf{r}}}\delta_{n,n^{\prime}}$
$\displaystyle=\delta(\boldsymbol{r}_{n}-\boldsymbol{r}_{n^{\prime}})$ (583)
$\displaystyle\lim_{\Delta_{\omega}\rightarrow
0}\frac{1}{\Delta_{\omega}}\delta_{\ell,\ell^{\prime}}$
$\displaystyle=\delta(\omega_{\ell}-\omega_{\ell^{\prime}})$ (584)

such that

$\left[f_{i}(\boldsymbol{r}_{n},\omega_{\ell}),f_{i^{\prime}}^{\dagger}(\boldsymbol{r}_{n^{\prime}},\omega_{\ell^{\prime}})\right]=\delta_{i,i^{\prime}}\delta(\boldsymbol{r}_{n}-\boldsymbol{r}_{n^{\prime}})\delta(\omega_{\ell}-\omega_{\ell^{\prime}})$
(585)

as intended. Therefore, we can write

$H_{F}=\sum_{i}\sum_{n}\sum_{\ell}\hbar\omega_{\ell}\tilde{f}_{i,n,\ell}^{\dagger}\tilde{f}_{i,n,\ell}.$
(586)

An eigenstate of this Hamiltonian is characterized by the occupation number of
each state $\left(i,n,\ell\right)$, $n_{\lambda,i,n,\ell}$. We represent such
state as

$\left|\left\\{n_{i_{j},n_{j},\ell_{j}}\right\\}\right\rangle=\left|n_{i_{1},n_{1},\ell_{1}},...,n_{i_{N},n_{N},\ell_{N}}\right\rangle,$
(587)

which has $n_{i_{j},n_{j},\ell_{j}}$ particles in state
$\left(i_{j},n_{j},\ell_{j}\right)$. Now let us compute

$\left\langle
f_{i^{\prime}}^{\dagger}(\boldsymbol{r}_{n^{\prime}},\omega_{\ell^{\prime}})f_{i}(\boldsymbol{r}_{n},\omega_{\ell})\right\rangle=\frac{1}{\Delta_{\boldsymbol{r}}\Delta_{\omega}}\left\langle\tilde{f}_{\lambda^{\prime}i^{\prime},n^{\prime},\ell^{\prime}}^{\dagger}\tilde{f}_{\lambda^{\prime}i^{\prime},n^{\prime},\ell^{\prime}}\right\rangle.$
(588)

Now we evaluate
$\left\langle\tilde{f}_{i,n,\ell}^{\dagger}\tilde{f}_{i,n,\ell}\right\rangle$.
For simplicity we group all the state labels into a superindex
$I=\left(i,n,\ell\right)$. Therefore, we want to evaluate

$\displaystyle\left\langle\tilde{f}_{I}^{\dagger}\tilde{f}_{I}\right\rangle$
$\displaystyle=\frac{\sum_{\left\\{n_{I^{\prime}}\right\\}}\left\langle\left\\{n_{I^{\prime}}\right\\}\right|e^{-\beta\sum_{J^{\prime}}\hbar\omega_{J^{\prime}}\tilde{f}_{J^{\prime}}^{\dagger}\tilde{f}_{J^{\prime}}}\tilde{f}_{I}^{\dagger}\tilde{f}_{I}\left|\left\\{n_{I^{\prime}}\right\\}\right\rangle}{\sum_{\left\\{n_{I^{\prime}}\right\\}}\left\langle\left\\{n_{I^{\prime}}\right\\}\right|e^{-\beta\sum_{J^{\prime}}\hbar\omega_{J^{\prime}}\tilde{f}_{J^{\prime}}^{\dagger}\tilde{f}_{J^{\prime}}}\left|\left\\{n_{I^{\prime}}\right\\}\right\rangle},$
(589)

where $\beta=1/k_{B}T$. Now we write, assuming that the system is populated by
$N$ bosons, that

$\displaystyle\sum_{\left\\{n_{I^{\prime}}\right\\}}\left\langle\left\\{n_{I^{\prime}}\right\\}\right|e^{-\beta\sum_{J^{\prime}}\hbar\omega_{J^{\prime}}\tilde{f}_{J^{\prime}}^{\dagger}\tilde{f}_{J^{\prime}}}\tilde{f}_{I}^{\dagger}\tilde{f}_{I}\left|\left\\{n_{I^{\prime}}\right\\}\right\rangle$
$\displaystyle=$
$\displaystyle\sum_{n_{I_{1}}=0}^{+\infty}...\sum_{n_{I}=0}^{+\infty}...\sum_{n_{I_{N}}=0}^{+\infty}\left\langle
n_{I_{1}},...,n_{I},...,n_{I_{N}}\right|e^{-\beta\hbar\omega_{I_{1}}\tilde{f}_{I_{1}}^{\dagger}\tilde{f}_{I_{1}}}...e^{-\beta\hbar\omega_{I}\tilde{f}_{I}^{\dagger}\tilde{f}_{I}}\tilde{f}_{I}^{\dagger}\tilde{f}_{I}...e^{-\beta\hbar\omega_{I_{1}}\tilde{f}_{I_{N}}^{\dagger}\tilde{f}_{I_{N}}}\left|n_{I_{1}},...,n_{I},...,n_{I_{N}}\right\rangle$
$\displaystyle=$
$\displaystyle\sum_{n_{I_{1}}=0}^{+\infty}...\sum_{n_{I}=0}^{+\infty}...\sum_{n_{I_{N}}=0}^{+\infty}e^{-\beta\hbar\omega_{I_{1}}n_{I_{1}}}...e^{-\beta\hbar\omega_{I}n_{I}}n_{I}...e^{-\beta\hbar\omega_{I_{1}}n_{I_{N}}}$
$\displaystyle=$
$\displaystyle\sum_{n_{I}=0}^{+\infty}e^{-\beta\hbar\omega_{I}n_{I}}n_{I}\prod_{n_{I^{\prime}}\neq
n_{I}}\sum_{n_{I^{\prime}}=0}^{+\infty}e^{-\beta\hbar\omega_{I^{\prime}}n_{I^{\prime}}}.$
(590)

Likewise

$\displaystyle\sum_{\left\\{n_{I^{\prime}}\right\\}}\left\langle\left\\{n_{I^{\prime}}\right\\}\right|e^{-\beta\sum_{J^{\prime}}\hbar\omega_{J^{\prime}}\tilde{f}_{J^{\prime}}^{\dagger}\tilde{f}_{J^{\prime}}}\left|\left\\{n_{I^{\prime}}\right\\}\right\rangle$
$\displaystyle=$
$\displaystyle\sum_{n_{I_{1}}=0}^{+\infty}...\sum_{n_{I}=0}^{+\infty}...\sum_{n_{I_{N}}=0}^{+\infty}\left\langle
n_{I_{1}},...,n_{I},...,n_{I_{N}}\right|e^{-\beta\hbar\omega_{I_{1}}\tilde{f}_{I_{1}}^{\dagger}\tilde{f}_{I_{1}}}...e^{-\beta\hbar\omega_{I}\tilde{f}_{I}^{\dagger}\tilde{f}_{I}}...e^{-\beta\hbar\omega_{I_{1}}\tilde{f}_{I_{N}}^{\dagger}\tilde{f}_{I_{N}}}\left|n_{I_{1}},...,n_{I},...,n_{I_{N}}\right\rangle$
$\displaystyle=$
$\displaystyle\sum_{n_{I_{1}}=0}^{+\infty}...\sum_{n_{I}=0}^{+\infty}...\sum_{n_{I_{N}}=0}^{+\infty}e^{-\beta\hbar\omega_{I_{1}}n_{I_{1}}}...e^{-\beta\hbar\omega_{I}n_{I}}...e^{-\beta\hbar\omega_{I_{1}}n_{I_{N}}}$
$\displaystyle=$
$\displaystyle\sum_{n_{I}=0}^{+\infty}e^{-\beta\hbar\omega_{I}n_{I}}\prod_{n_{I^{\prime}}\neq
n_{I}}\sum_{n_{I^{\prime}}=0}^{+\infty}e^{-\beta\hbar\omega_{I^{\prime}}n_{I^{\prime}}}.$
(591)

Therefore

$\displaystyle\left\langle\tilde{f}_{I}^{\dagger}\tilde{f}_{I}\right\rangle$
$\displaystyle=\frac{\sum_{n_{I}=0}^{+\infty}e^{-\beta\hbar\omega_{I}n_{I}}n_{I}\prod_{n_{I^{\prime}}\neq
n_{I}}\sum_{n_{I^{\prime}}=0}^{+\infty}e^{-\beta\hbar\omega_{I^{\prime}}n_{I^{\prime}}}}{\sum_{n_{I}=0}^{+\infty}e^{-\beta\hbar\omega_{I}n_{I}}\prod_{n_{I^{\prime}}\neq
n_{I}}\sum_{n_{I^{\prime}}=0}^{+\infty}e^{-\beta\hbar\omega_{I^{\prime}}n_{I^{\prime}}}}$
$\displaystyle=\frac{\sum_{n_{I}=0}^{+\infty}e^{-\beta\hbar\omega_{I}n_{I}}n_{I}}{\sum_{n_{I}=0}^{+\infty}e^{-\beta\hbar\omega_{I}n_{I}}}$
$\displaystyle=\frac{1}{\sum_{n_{I}=0}^{+\infty}e^{-\beta\hbar\omega_{I}n_{I}}}\left(-\frac{\partial}{\partial\beta\hbar\omega_{I}}\sum_{n_{I}=0}^{+\infty}e^{-\beta\hbar\omega_{I}n_{I}}\right).$
(592)

Since

$\displaystyle\sum_{n_{I}=0}^{+\infty}e^{-\beta\hbar\omega_{I}n_{I}}$
$\displaystyle=1+\left(e^{-\beta\hbar\omega_{I}}\right)+\left(e^{-\beta\hbar\omega_{I}}\right)^{2}+...$
$\displaystyle=\frac{1}{1-e^{-\beta\hbar\omega_{I}}}.$ (593)

Now

$\frac{\partial}{\partial
x}\left(\frac{1}{1-e^{-x}}\right)=-\frac{e^{-x}}{\left(1-e^{-x}\right)^{2}}.$
(594)

Therefore

$\displaystyle\left\langle\tilde{f}_{I}^{\dagger}\tilde{f}_{I}\right\rangle$
$\displaystyle=\left(1-e^{-\beta\hbar\omega_{I}}\right)\left(\frac{e^{-\beta\hbar\omega_{I}}}{\left(1-e^{-\beta\hbar\omega_{I}}\right)^{2}}\right)$
$\displaystyle=\frac{e^{-\beta\hbar\omega_{I}}}{1-e^{-\beta\hbar\omega_{I}}}$
$\displaystyle=\frac{1}{e^{\beta\hbar\omega_{I}}-1}$
$\displaystyle=\bar{n}_{\text{spp}}(\omega_{I}).$ (595)

If $I\neq J$, then we have

$\left\langle\tilde{f}_{I}^{\dagger}\tilde{f}_{J}\right\rangle=0,I\neq J.$
(596)

Therefore, we can write

$\left\langle\tilde{f}_{i^{\prime},n^{\prime},\ell^{\prime}}^{\dagger}\tilde{f}_{i,n,\ell}\right\rangle=\delta_{i,i^{\prime}}\delta_{n,n^{\prime}}\delta_{\ell,\ell^{\prime}}\bar{n}_{\text{spp}}(\omega_{\ell}).$
(597)

Therefore

$\left\langle
f_{i^{\prime}}^{\dagger}(\boldsymbol{r}_{n^{\prime}},\omega_{\ell^{\prime}})f_{i}(\mathbf{r}_{n},\omega_{\ell})\right\rangle=\frac{1}{\Delta_{\mathbf{r}}\Delta_{\omega}}\delta_{i,i^{\prime}}\delta_{n,n^{\prime}}\delta_{\ell,\ell^{\prime}}\bar{n}_{\text{spp}}(\omega_{\ell})$
(598)

As before, in the limit of $\Delta_{\mathbf{r}},\Delta_{\omega}\rightarrow 0$,
we obtain

$\left\langle
f_{i^{\prime}}^{\dagger}(\mathbf{r}_{n^{\prime}},\omega_{\ell^{\prime}})f_{i}(\mathbf{r}_{n},\omega_{\ell})\right\rangle=\delta_{i,i^{\prime}}\delta(\boldsymbol{r}_{n}-\boldsymbol{r}_{n^{\prime}})\delta(\omega_{\ell}-\omega_{\ell^{\prime}})\bar{n}_{\text{spp}}(\omega_{\ell}),$
(599)

or

$\left\langle
f_{i^{\prime}}^{\dagger}(\boldsymbol{r}^{\prime},\omega^{\prime})f_{i}(\boldsymbol{r},\omega)\right\rangle=\delta_{i,i^{\prime}}\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})\delta(\omega-\omega^{\prime})\bar{n}_{\text{spp}}(\omega).$
(600)

This is the result presented in the main text without further justification,
but this brief derivation shows that it is indeed valid.

##  Appendix J Details of the calculation of the dynamics of two qubits
coupled to the plasmonic bath

###  J.1 Obtaining the dynamics via the Schrodinger equation

The Hamiltonian for this system can be constructed analogously to the previous
results. We consider the Hamiltonian corresponding to each atom $H_{0,1}$ and
$H_{0,2}$ as well as the SPP field $H_{\text{spp}}.$ In addition, we must
consider the interaction terms between the plasmon field and both qubits in
the dipole approximation $H_{d,1}$ and $H_{d,2}$. The Hamiltonian is written
as

$H=\sum_{\alpha=1,2}H_{0,\alpha}+H_{\text{spp}}+\sum_{\alpha=1,2}H_{d,\alpha}.$
(601)

To simplify the following equations, we introduce the notation

$\begin{cases}\hbar\mathcal{G}_{\alpha}^{-}(\boldsymbol{r},\boldsymbol{\rho};\omega)\equiv\omega\beta(\omega)\boldsymbol{d}_{\alpha}\cdot\bar{\bar{G}}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega),\\\
\hbar\mathcal{G}_{\alpha}^{+}(\boldsymbol{r},\boldsymbol{\rho};\omega)\equiv\omega\beta^{*}(\omega)\bar{\bar{G}}^{\dagger}(\boldsymbol{r},\boldsymbol{\rho},z=0;\omega)\cdot\boldsymbol{d}_{\alpha}^{\dagger}.\end{cases}$
(602)

The components of the Hamiltonian are therefore, in the rotating wave
approximation, given by

$\begin{cases}H_{0,\alpha}=\hbar\omega_{0}\mathds{1_{\text{spp}}}\sigma_{\alpha}^{+}\sigma_{\alpha}^{-},\\\
H_{\text{spp}}=\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\hbar\omega\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\hat{\text{{f}}}(\boldsymbol{\rho},\omega),\\\
H_{d,\alpha}=i\hbar\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\left[\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\cdot\mathcal{G}_{\alpha}^{+}(\boldsymbol{r}_{\alpha},\boldsymbol{\rho};\omega)\sigma_{\alpha}^{-}-\mathcal{G}_{\alpha}^{-}(\boldsymbol{r}_{\alpha},\boldsymbol{\rho};\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho},\omega)\sigma_{\alpha}^{+}\right].\end{cases}.$
(603)

As in the previous notes, we want to solve the Schrodinger equation

$i\hbar\frac{\partial}{\partial
t}\left|\psi(t)\right\rangle=H\left|\psi(t)\right\rangle.$ (604)

We considet that the system is coupled to a single excitation of the field.
This means that initially one of the atoms (label it as atom $\alpha$) will be
in the excited state while the other will be in the ground state and there
will be no SPPs in the graphene. We write an ansatz for the wave-function as

$\left|\psi(t)\right\rangle=\sum_{\alpha=1,2}C_{\alpha}(t)e^{-i\omega_{\alpha}t}\left|\alpha\right\rangle\left|0\right\rangle+\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\boldsymbol{C}_{g}(\boldsymbol{\rho},t;\omega)e^{-i\omega
t}\cdot\left|g\right\rangle\left|\mathds{1}(\boldsymbol{\rho},\omega)\right\rangle.$
(605)

We now substitute this into the Schrodinger equation on both sides. The left-
hand side, involving the time derivative, yields

$\displaystyle i\hbar\frac{\partial}{\partial t}\left|\psi(t)\right\rangle=$
$\displaystyle\sum_{\alpha^{\prime}=1,2}\left[\hbar\omega_{\alpha^{\prime}}C_{\alpha^{\prime}}(t)e^{-i\omega_{\alpha^{\prime}}t}+i\hbar\dot{C}_{\alpha^{\prime}}(t)e^{-i\omega_{\alpha^{\prime}}t}\right]\left|\alpha^{\prime}\right\rangle\left|0\right\rangle$
$\displaystyle+\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\left[\hbar\omega\boldsymbol{C}_{g}(\boldsymbol{\rho},t;\omega)e^{-i\omega
t}+i\hbar\dot{\boldsymbol{C}}_{g}(\boldsymbol{\rho},t;\omega)e^{-i\omega
t}\right]\cdot\left|g\right\rangle\left|\mathds{1}(\boldsymbol{\rho};\omega)\right\rangle.$
(606)

On the other hand, for the right hand side, we have

$\displaystyle\hat{H}\left|\psi(t)\right\rangle=$
$\displaystyle\sum_{\alpha^{\prime}=1,2}\hbar\omega_{\alpha^{\prime}}C_{\alpha^{\prime}}(t)e^{-i\omega_{\alpha^{\prime}}t}$
$\displaystyle+\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\hbar\omega\boldsymbol{C}_{g}(t)e^{-i\omega
t}\left|g\right\rangle\left|\mathds{1}(\boldsymbol{\rho};\omega)\right\rangle$
$\displaystyle+i\hbar\sum_{\alpha^{\prime}=1,2}\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}_{\alpha^{\prime}}(t)e^{-i\omega_{\alpha^{\prime}}t}\left|g\right\rangle\left|\mathds{1}(\boldsymbol{\rho};\omega)\right\rangle\cdot\mathcal{G}_{\alpha^{\prime}}^{+}(\boldsymbol{r}_{\alpha^{\prime}},\boldsymbol{\rho};\omega)$
$\displaystyle-i\hbar\sum_{\alpha^{\prime}=1,2}\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}_{\alpha^{\prime}}^{-}(\boldsymbol{r}_{\alpha^{\prime}},\boldsymbol{\rho};\omega)\cdot\boldsymbol{C}_{g}(\boldsymbol{\rho},t;\omega)e^{-i\omega
t}\left|\alpha^{\prime}\right\rangle\left|0\right\rangle.$

Thus, we see that, equating both sides, many terms cancel, leaving us with

$\displaystyle
i\hbar\sum_{\alpha^{\prime}=1,2}\dot{C}_{\alpha^{\prime}}(t)e^{-i\omega_{\alpha^{\prime}}t}\left|\alpha^{\prime}\right\rangle\left|0\right\rangle+i\hbar\int
d^{2}\boldsymbol{\rho}\int_{0}^{\infty}d\omega\dot{\boldsymbol{C}}_{g}(\boldsymbol{\rho},t;\omega)e^{-i\omega
t}\cdot\left|g\right\rangle\left|\mathds{1}(\boldsymbol{\rho};\omega)\right\rangle$
$\displaystyle=i\hbar\sum_{\alpha^{\prime}=1,2}\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}_{\alpha^{\prime}}(t)e^{-i\omega_{\alpha}t}\left|g\right\rangle\left|\mathds{1}(\boldsymbol{\rho};\omega)\right\rangle\cdot\mathcal{G}_{\alpha^{\prime}}^{+}(\boldsymbol{r}_{\alpha^{\prime}},\boldsymbol{\rho};\omega)$
$\displaystyle-i\hbar\sum_{\alpha^{\prime}=1,2}\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}_{\alpha^{\prime}}^{-}(\boldsymbol{r}_{\alpha^{\prime}},\boldsymbol{\rho};\omega)\cdot\boldsymbol{C}_{g}(\boldsymbol{\rho},t;\omega)\left|\alpha^{\prime}\right\rangle\left|0\right\rangle.$
(607)

Acting on the left on both sides of the equation with the state
$\left\langle\alpha\right|\left\langle 0\right|$, we get

$\dot{C}_{\alpha}(t)e^{-i\omega_{\alpha}t}=-\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}_{\alpha}^{-}(\boldsymbol{r}_{\alpha},\boldsymbol{\rho};\omega)\cdot\boldsymbol{C}_{g}(\boldsymbol{\rho},t;\omega)e^{-i\omega
t},$ (608)

while if we act only on the state with $\left\langle
g\right|\left\langle\mathds{1}(\boldsymbol{\rho};\omega)\right|$, we obtain

$\dot{\boldsymbol{C}}_{g}(\boldsymbol{\rho},t;\omega)e^{-i\omega
t}=\sum_{\alpha^{\prime}=1,2}C_{\alpha^{\prime}}(t)e^{-i\omega_{\alpha^{\prime}}t}\mathcal{G}_{\alpha^{\prime}}^{+}(\boldsymbol{r}_{\alpha^{\prime}},\boldsymbol{\rho};\omega).$
(609)

We can now integrate this final expression in time, and via performing the
first Markov approximation, isolate the integral over the exponentials by
writing $C_{\alpha^{\prime}}(t_{1})\approx C_{\alpha}(t)$.

$\boldsymbol{C}_{g}(\boldsymbol{\rho},t;\omega)=\sum_{\alpha^{\prime}=1,2}C_{\alpha^{\prime}}(t)\mathcal{G}_{\alpha^{\prime}}^{+}(\boldsymbol{r}_{\alpha^{\prime}},\boldsymbol{\rho};\omega)\int_{0}^{t}dt_{1}e^{i\left(\omega-\omega_{\alpha^{\prime}}\right)t_{1}}.$
(610)

We now take this expression and substitute it back into the equation for
$C_{\alpha}(t)$ to obtain

$\dot{C}_{\alpha}(t)=-\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}_{\alpha}^{-}(\boldsymbol{r}_{\alpha},\boldsymbol{\rho};\omega)e^{-i(\omega-\omega_{\alpha})t}\cdot\left[\sum_{\alpha^{\prime}=1,2}C_{\alpha^{\prime}}(t)\mathcal{G}_{\alpha^{\prime}}^{+}(\boldsymbol{r}_{\alpha^{\prime}},\boldsymbol{\rho};\omega)\int_{0}^{t}dt_{1}e^{-i(\omega_{\alpha^{\prime}}-\omega)t_{1}}\right].$
(611)

Grouping the terms together everything together and introducing the notation

$\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\alpha^{\prime}}(\omega)\cdot\boldsymbol{d}_{\alpha^{\prime}}=\int
d^{2}\boldsymbol{\rho}\mathcal{G}_{\alpha}^{-}(\boldsymbol{r}_{\alpha},\boldsymbol{\rho};\omega)\mathcal{G}_{\alpha^{\prime}}^{+}(\boldsymbol{r}_{\alpha^{\prime}},\boldsymbol{\rho};\omega),$

we obtain

$\dot{C}_{\alpha}(t)=-\int_{0}^{\infty}d\omega\sum_{\alpha^{\prime}=1,2}\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\alpha^{\prime}}(\omega)\cdot\boldsymbol{d}_{\alpha^{\prime}}\int_{0}^{t}dt_{1}e^{-i(\omega_{\alpha^{\prime}}-\omega)t_{1}-i(\omega-\omega_{\alpha})t}.$
(612)

For simplicity we admit that the two-level atoms are identical, i.e.,
$\omega_{1}=\omega_{2}$. This allows us to write

$\dot{C}_{\alpha}(t)=-\int_{0}^{\infty}d\omega\sum_{\alpha^{\prime}=1,2}C_{\alpha^{\prime}}(t)\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\alpha^{\prime}}(\omega)\cdot\boldsymbol{d}_{\alpha^{\prime}}\int_{0}^{t}dt_{1}e^{-i(\omega_{\alpha^{\prime}}-\omega)\left(t_{1}-t\right)},$
(613)

The integral over time by making the change of variables $\tau=t-t_{1}$ as
well as using the Sokhotski-Plemelj identity, where in the third line we
extended the region of integration via the second Markov approximation, and
then applied the Sokhtski-Plemelj theorem. We are left with the differential
equation

$\displaystyle\dot{C}_{\alpha}(t)$
$\displaystyle=-\sum_{\alpha^{\prime}=1,2}C_{\alpha^{\prime}}(t)\left[-i\mathcal{P}\int_{0}^{\infty}d\omega\frac{\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\alpha^{\prime}}(\omega)\cdot\boldsymbol{d}_{\alpha^{\prime}}}{\omega-\omega_{\alpha}}+\pi\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\alpha^{\prime}}(\omega_{\alpha})\cdot\boldsymbol{d}_{\alpha^{\prime}}\right]$
$\displaystyle\equiv-\sum_{\alpha^{\prime}=1,2}C_{\alpha^{\prime}}(t)\left[-i\Delta_{\alpha\alpha^{\prime}}+\Gamma_{\alpha\alpha^{\prime}}\right].$
(614)

Explicitly, for our two-level system, we have

$\begin{cases}\dot{C}_{1}(t)=-C_{1}(t)\left[\Gamma_{11}-i\Delta_{11}\right]-C_{2}(t)\left[\Gamma_{12}-i\Delta_{12}\right],\\\
\dot{C}_{2}(t)=-C_{1}(t)\left[\Gamma_{21}-i\Delta_{21}\right]-C_{2}(t)\left[\Gamma_{22}-i\Delta_{22}\right].\end{cases}.$
(615)

We can construct the states

$\begin{cases}\left|+\right\rangle=\frac{1}{\sqrt{2}}\left(\left|1\right\rangle+\left|2\right\rangle\right),\\\
\left|-\right\rangle=\frac{1}{\sqrt{2}}\left(\left|1\right\rangle-\left|2\right\rangle\right).\end{cases},$

with amplitudes

$\begin{cases}C_{+}(t)=\frac{1}{\sqrt{2}}\left[C_{1}(t)+C_{2}(t)\right],\\\
C_{-}(t)=\frac{1}{\sqrt{2}}\left[C_{1}(t)-C_{2}(t)\right].\end{cases}$ (616)

If we now admit that
$\Delta_{11}=\Delta_{22}\equiv\Delta,$$\Gamma_{11}=\Gamma_{22}\equiv\Gamma$
and $\Delta_{21}=\Delta_{12}\equiv
g_{12},\Gamma_{12}=\Gamma_{21}\equiv\Gamma_{12}$, we have

$\begin{cases}\dot{C}_{1}(t)=-C_{1}(t)\left[\Gamma-i\Delta\right]-C_{2}(t)\left[\Gamma_{12}-ig_{12}\right],\\\
\dot{C}_{2}(t)=-C_{1}(t)\left[\Gamma_{12}-ig_{12}\right]-C_{2}(t)\left[\Gamma-i\Delta\right].\end{cases}.$
(617)

$\begin{cases}\dot{C}_{+}(t)=C_{+}(t)\left[i\left(\Delta+g_{12}\right)-\left(\Gamma+\Gamma_{12}\right)\right],\\\
\dot{C}_{-}(t)=C_{-}(t)\left[i\left(\Delta-
g_{12}\right)-\left(\Gamma-\Gamma_{12}\right)\right].\end{cases}$ (618)

We note also that the decay rates are given by the dyadic Green’s function as

$\displaystyle\Gamma_{\alpha\alpha^{\prime}}$
$\displaystyle=\pi\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\alpha^{\prime}}(\omega)\cdot\boldsymbol{d}_{\alpha^{\prime}}$
$\displaystyle=\frac{\pi\bar{\varepsilon}}{(2\pi)^{3}\hbar\varepsilon_{0}c^{2}}\omega_{\alpha}^{2}\boldsymbol{d}_{\alpha}\cdot\text{Im}[\bar{\bar{G}}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\alpha^{\prime}};\omega_{\alpha})]\cdot\boldsymbol{d}_{\alpha\text{'}}.$
(619)

This yields the solutions

$\begin{cases}C_{+}(t)=C_{+}^{0}e^{i\left(\Delta+g_{12}\right)t-\left(\Gamma+\Gamma_{12}\right)t}\\\
C_{-}(t)=C_{-}^{0}e^{i\left(\Delta-
g_{12}\right)t-\left(\Gamma-\Gamma_{12}\right)t}\end{cases}$ (620)

and the probability is given by

$\begin{cases}P_{+}(t)=P_{+}^{0}e^{-\left(\gamma+\gamma_{12}\right)t},\\\
P_{-}(t)=P_{-}^{0}e^{-\left(\gamma-\gamma_{12}\right)t}.\end{cases}$ (621)

We can now imagine that initially $C_{1}(t=0)=1$ and $C_{2}(t=0)=0$ in order
to look at the dynamics of the system. In particular, we can write

$\begin{cases}C_{1}(t)=\frac{1}{\sqrt{2}}\left[C_{+}^{0}e^{-i\left(\Delta+g_{12}\right)t-\left(\Gamma+\Gamma_{12}\right)t}+C_{-}^{0}e^{-i\left(\Delta-
g_{12}\right)t-\left(\Gamma-\Gamma_{12}\right)t}\right],\\\
C_{2}(t)=\frac{1}{\sqrt{2}}\left[C_{+}^{0}e^{-i\left(\Delta+g_{12}\right)t-\left(\Gamma+\Gamma_{12}\right)t}-C_{-}^{0}e^{-i\left(\Delta-
g_{12}\right)t-\left(\Gamma-\Gamma_{12}\right)t}\right].\end{cases}$

or using the initial conditions

$C_{2}(t=0)=0\Rightarrow C_{+}^{0}=C_{-}^{0}\equiv C_{1}(t=0).$ (622)

Therefore, grouping the necessary terms together, yields

$\begin{cases}C_{1}(t)=\frac{1}{\sqrt{2}}C_{1}(0)e^{-i\Delta-\Gamma}\left(e^{-i\left(g_{12}-i\Gamma_{12}\right)t}+e^{i\left(g_{12}-i\Gamma_{12}\right)t}\right),\\\
C_{2}(t)=\frac{1}{\sqrt{2}}C_{2}(0)e^{-i\Delta-\Gamma}\left(e^{-i\left(g_{12}-i\Gamma_{12}\right)t}-e^{i\left(g_{12}-i\Gamma_{12}\right)t}\right).\end{cases}$
(623)

and the corresponding probabilities are

$\begin{cases}P_{1}(t)&=P_{1}^{0}e^{-\gamma}\left(\cos\left(2g_{12}t\right)+\cosh\left(\gamma_{12}t\right)\right),\\\
P_{2}(t)&=P_{2}^{0}e^{-\gamma}\left(\cos\left(2g_{12}t\right)-\cosh\left(\gamma_{12}t\right)\right).\end{cases}$
(624)

Plots of the corresponding probabilities can be found in the main text in
figures 5a and 5b in the main text.

###  J.2 Deriving the Lindblad equation for two qubits

To obtain the dynamics of two coupled qubits the procedure is yet again very
similar to the previous calculations. We start with the Hamiltonian, as usual,
which can be written as

$\hat{H}=\sum_{\alpha=1,2}\hat{H}_{\alpha}+\hat{H}_{\text{spp}}+\hat{V},$
(625)

with each component given in equation (603). We then have the Redfield
equation and with

$\begin{cases}\beta_{\alpha,1}=\sigma_{\alpha}^{-},\\\
\beta_{\alpha,2}=\sigma_{\alpha}^{-},\\\
\Gamma_{\alpha,1}=\Gamma_{\alpha}^{\dagger}=i\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\hat{\text{{f}}}^{\dagger}(\boldsymbol{\rho},\omega)\cdot\mathcal{G}^{+}(\boldsymbol{r}_{\alpha},\boldsymbol{\rho};\omega),\\\
\Gamma_{\alpha,2}=\Gamma_{\alpha}=-i\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r}_{\alpha},\boldsymbol{\rho};\omega)\cdot\hat{\text{{f}}}(\boldsymbol{\rho},\omega).\end{cases},$
(626)

which we can write as

$\displaystyle\frac{d\hat{\rho}_{S}}{dt}=$
$\displaystyle-\frac{1}{\hbar^{2}}\sum_{\alpha\beta}\sum_{ij}\int_{t_{0}}^{t}d\tau\left(\beta_{\alpha,i}(t)\beta_{\beta,j}(\tau)\hat{\rho}_{S}(t)-\beta_{\beta,j}(\tau)\hat{\rho}_{S}(t)\beta_{\alpha,i}(t)\right)\left\langle\Gamma_{\alpha,i}(t)\Gamma_{\beta,j}(\tau)\right\rangle_{R}$
$\displaystyle+\left(\hat{\rho}_{S}(t)\beta_{\beta,j}(\tau)\beta_{\alpha,i}(t)-\beta_{\alpha,i}(t)\hat{\rho}_{S}(t)\beta_{\beta,j}(\tau)\right)\left\langle\Gamma_{\beta,j}(\tau)\Gamma_{\alpha,i}(t)\right\rangle_{R}.$
(627)

By analogy with the previous results, the reservoir correlation functions are

$\begin{cases}\left\langle\Gamma_{\alpha,1}(t)\Gamma_{\beta,1}(\tau)\right\rangle_{R}=\left\langle\Gamma_{\alpha,2}(t)\Gamma_{\beta,2}(\tau)\right\rangle_{R}=0,\\\
\left\langle\Gamma_{\alpha,1}(t)\Gamma_{\beta,2}(\tau)\right\rangle_{R}=\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r}_{\beta},\boldsymbol{\rho};\omega)\mathcal{G}^{+}(\boldsymbol{r}_{\alpha},\boldsymbol{\rho};\omega)\bar{n}_{\text{spp}}(\omega)e^{-i\omega(\tau-t)},\\\
\left\langle\Gamma_{\alpha,2}(t)\Gamma_{\beta,1}(\tau)\right\rangle_{R}=\int_{0}^{\infty}d\omega\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r}_{\alpha},\boldsymbol{\rho};\omega)\mathcal{G}^{+}(\boldsymbol{r}_{\beta},\boldsymbol{\rho};\omega)\left[\bar{n}_{\text{spp}}(\omega)+1\right]e^{i\omega(\tau-t)},\end{cases}$
(628)

and if we define, for ease of notation the kernel

$\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\beta}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\equiv\frac{1}{\hbar^{2}}\int
d^{2}\boldsymbol{\rho}\mathcal{G}^{-}(\boldsymbol{r}_{\alpha},\boldsymbol{\rho};\omega)\cdot\mathcal{G}^{+}(\boldsymbol{r}_{\beta},\boldsymbol{\rho};\omega),$
(629)

and the integral

$\mathcal{I}_{SP}(\omega)=\int_{0}^{\infty}d\tau
e^{i(\omega-\omega_{0})(\tau-t)}$ (630)

we obtain, evaluating the sums of equation (627) and in the Markov
approximation, the result

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\int_{0}^{\infty}d\omega\sum_{\alpha\beta}\left(\sigma_{\alpha}^{+}\sigma_{\beta}^{-}\hat{\rho}_{S}^{I}(t)-\sigma_{\beta}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{+}\right)\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\beta}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\left[\bar{n}_{\text{spp}}(\omega)+1\right]\mathcal{I}(\omega)$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{\text{+}}\sigma_{\alpha}^{-}-\sigma_{\alpha}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{+}\right)\boldsymbol{d}_{\beta}\cdot\bar{\mathcal{K}}_{\beta\alpha}(\boldsymbol{r}_{\beta},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\left[\bar{n}_{\text{spp}}(\omega)+1\right]\mathcal{I}^{*}(\omega)$
$\displaystyle+\left(\sigma_{\alpha}^{-}\sigma_{\beta}^{+}\hat{\rho}_{S}^{I}(t)-\sigma_{\beta}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{-}\right)\boldsymbol{d}_{\alpha}\bar{\mathcal{K}}_{\alpha\beta}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\bar{n}_{\text{spp}}(\omega)\mathcal{I}^{*}(\omega)$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{-}\sigma_{\alpha}^{+}-\sigma_{\alpha}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{-}\right)\boldsymbol{d}_{\beta}\cdot\bar{\mathcal{K}}_{\beta\alpha}(\boldsymbol{r}_{\beta},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\bar{n}_{\text{spp}}(\omega)\mathcal{I}(\omega).$
(631)

Due to the summation, this equation will contain terms in $\alpha=\beta=1,2$
coming from the evolution of each qubit as well as cross terms
$\alpha=1,\beta=2$ or $\alpha=2,\beta=1$ coming from the interaction between
the qubits. The equation will be written as

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\int_{0}^{\infty}d\omega\sum_{\alpha\beta}\left(\sigma_{\alpha}^{+}\sigma_{\beta}^{-}\hat{\rho}_{S}^{I}(t)-\sigma_{\beta}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{+}\right)\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\beta}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\left[\bar{n}_{\text{spp}}(\omega)+1\right]\mathcal{I}(\omega)$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{\text{+}}\sigma_{\alpha}^{-}-\sigma_{\alpha}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{+}\right)\boldsymbol{d}_{\beta}\cdot\bar{\mathcal{K}}_{\beta\alpha}(\boldsymbol{r}_{\beta},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\left[\bar{n}_{\text{spp}}(\omega)+1\right]\mathcal{I}^{*}(\omega)$
$\displaystyle+\left(\sigma_{\alpha}^{-}\sigma_{\beta}^{+}\hat{\rho}_{S}^{I}(t)-\sigma_{\beta}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{-}\right)\boldsymbol{d}_{\alpha}\bar{\mathcal{K}}_{\alpha\beta}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\bar{n}_{\text{spp}}(\omega)\mathcal{I}^{*}(\omega)$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{-}\sigma_{\alpha}^{+}-\sigma_{\alpha}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{-}\right)\boldsymbol{d}_{\beta}\cdot\bar{\mathcal{K}}_{\beta\alpha}(\boldsymbol{r}_{\beta},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\bar{n}_{\text{spp}}(\omega)\mathcal{I}(\omega).$
(632)

Notice that we can split the sums into terms where $\alpha=\beta$ and where
$\alpha\neq\beta$. We wrtite the following

$\displaystyle\frac{d\hat{\rho}_{S}^{I}}{dt}=$
$\displaystyle-\int_{0}^{\infty}d\omega\sum_{\alpha}\left(\sigma_{\alpha}^{+}\sigma_{\alpha}^{-}\hat{\rho}_{S}^{I}(t)-\sigma_{\alpha}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{+}\right)\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\alpha}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\left[\bar{n}_{\text{spp}}(\omega)+1\right]\mathcal{I}(\omega)$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{\text{+}}\sigma_{\alpha}^{-}-\sigma_{\alpha}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{+}\right)\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\alpha}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\left[\bar{n}_{\text{spp}}(\omega)+1\right]\mathcal{I}^{*}(\omega)$
$\displaystyle+\left(\sigma_{\alpha}^{-}\sigma_{\alpha}^{+}\hat{\rho}_{S}^{I}(t)-\sigma_{\alpha}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{-}\right)\boldsymbol{d}_{\alpha}\bar{\mathcal{K}}_{\alpha\alpha}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\bar{n}_{\text{spp}}(\omega)\mathcal{I}^{*}(\omega)$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{-}\sigma_{\alpha}^{+}-\sigma_{\alpha}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{-}\right)\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\alpha}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\bar{n}_{\text{spp}}(\omega)\mathcal{I}(\omega)$
$\displaystyle-\int_{0}^{\infty}d\omega\sum_{\alpha\neq\beta}\left(\sigma_{\alpha}^{+}\sigma_{\beta}^{-}\hat{\rho}_{S}^{I}(t)-\sigma_{\beta}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{+}\right)\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\beta}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\left[\bar{n}_{\text{spp}}(\omega)+1\right]\mathcal{I}(\omega)$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{\text{+}}\sigma_{\alpha}^{-}-\sigma_{\alpha}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{+}\right)\boldsymbol{d}_{\beta}\cdot\bar{\mathcal{K}}_{\beta\alpha}(\boldsymbol{r}_{\beta},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\left[\bar{n}_{\text{spp}}(\omega)+1\right]\mathcal{I}^{*}(\omega)$
$\displaystyle+\left(\sigma_{\alpha}^{-}\sigma_{\beta}^{+}\hat{\rho}_{S}^{I}(t)-\sigma_{\beta}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{-}\right)\boldsymbol{d}_{\alpha}\bar{\mathcal{K}}_{\alpha\beta}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\bar{n}_{\text{spp}}(\omega)\mathcal{I}^{*}(\omega)$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{-}\sigma_{\alpha}^{+}-\sigma_{\alpha}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{-}\right)\boldsymbol{d}_{\beta}\cdot\bar{\mathcal{K}}_{\beta\alpha}(\boldsymbol{r}_{\beta},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\bar{n}_{\text{spp}}(\omega)\mathcal{I}(\omega).$
(633)

Performing the same analysis we had done for the case of the single qubit, we
find that the first summation becomes, with the jump operators

$\begin{cases}\hat{L}_{\alpha,1}=\sqrt{\Gamma\bar{n}_{\text{spp}}(\omega_{0})}\sigma_{\alpha}^{+},\\\
\hat{L}_{\alpha,2}=\sqrt{\Gamma\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]}\sigma_{\alpha}^{-}.\end{cases}$
(634)

$\displaystyle\sum_{\alpha}\left(2L_{\alpha,2}\hat{\rho}_{S}(t)L_{\alpha,2}^{\dagger}-\left\\{L_{\alpha,2}^{\dagger}L_{\alpha,2},\hat{\rho}_{S}(t)\right\\}\right)$
$\displaystyle+\left(2L_{\alpha,1}\hat{\rho}_{S}(t)L_{\alpha,1}^{\dagger}-\left\\{L_{\alpha,1}^{\dagger}L_{\alpha,1},\hat{\rho}_{S}(t)\right\\}\right)-i(\omega_{0}+\Delta_{\alpha}+2\Delta^{\prime}_{\alpha})\left[\sigma_{\alpha}^{+}\sigma_{\alpha}^{-},\hat{\rho}_{S}(t)\right],$
(635)

whereas the for the second expression we can split it apart and write

$\displaystyle\sum_{\alpha\neq\beta}\int_{0}^{\infty}d\omega\left(\sigma_{\alpha}^{+}\sigma_{\beta}^{-}\hat{\rho}_{S}^{I}(t)-\sigma_{\beta}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{+}\right)\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\beta}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\left[\bar{n}_{\text{spp}}(\omega)+1\right]\mathcal{I}(\omega)$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{\text{+}}\sigma_{\alpha}^{-}-\sigma_{\alpha}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{+}\right)\boldsymbol{d}_{\beta}\cdot\bar{\mathcal{K}}_{\beta\alpha}(\boldsymbol{r}_{\beta},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\left[\bar{n}_{\text{spp}}(\omega)+1\right]\mathcal{I}^{*}(\omega)$
$\displaystyle+\left(\sigma_{\alpha}^{-}\sigma_{\beta}^{+}\hat{\rho}_{S}^{I}(t)-\sigma_{\beta}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{-}\right)\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}_{\alpha\beta}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\bar{n}_{\text{spp}}(\omega)\mathcal{I}^{*}(\omega)$
(636)
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{-}\sigma_{\alpha}^{+}-\sigma_{\alpha}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{-}\right)\boldsymbol{d}_{\beta}\cdot\bar{\mathcal{K}}_{\beta\alpha}(\boldsymbol{r}_{\beta},\boldsymbol{r}_{\alpha},\omega)\cdot\boldsymbol{d}_{\alpha}\bar{n}_{\text{spp}}(\omega)\mathcal{I}(\omega).$

We can make our usual definitions

$\begin{cases}\Gamma_{\alpha\beta}\equiv\pi\int_{0}^{\infty}d\omega\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\delta(\omega-\omega_{0}),\\\
\Delta_{\alpha\beta}\equiv\mathcal{P}\int_{0}^{\infty}d\omega\frac{\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}}{\omega_{0}-\omega},\\\
\Delta^{\prime}_{\alpha\beta}\equiv\mathcal{P}\int_{0}^{\infty}d\omega\frac{\boldsymbol{d}_{\alpha}\cdot\bar{\mathcal{K}}(\boldsymbol{r}_{\alpha},\boldsymbol{r}_{\beta},\omega)\cdot\boldsymbol{d}_{\beta}\bar{n}_{\text{spp}}(\omega)}{\omega_{0}-\omega}.\end{cases}$
(637)

and write

$\displaystyle-\sum_{\alpha\neq\beta}\left(\sigma_{\alpha}^{+}\sigma_{\beta}^{-}\hat{\rho}_{S}^{I}(t)-\sigma_{\beta}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{+}\right)\left[\Gamma_{\alpha\beta}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]+i\left(\Delta_{\alpha\beta}+\Delta^{\prime}_{\alpha\beta}\right)\right]$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{\text{+}}\sigma_{\alpha}^{-}-\sigma_{\alpha}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{+}\right)\left[\Gamma_{\beta\alpha}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]-i\left(\Delta_{\beta\alpha}+\Delta^{\prime}_{\beta\alpha}\right)\right]$
$\displaystyle+\left(\sigma_{\alpha}^{-}\sigma_{\beta}^{+}\hat{\rho}_{S}^{I}(t)-\sigma_{\beta}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\alpha}^{-}\right)\left[\Gamma_{\alpha\beta}\left(\bar{n}_{\text{spp}}(\omega_{0})+1\right)-i\Delta^{\prime}_{\alpha\beta}\right]$
$\displaystyle+\left(\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{-}\sigma_{\alpha}^{+}-\sigma_{\alpha}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{\beta}^{-}\right)\left[\Gamma_{\beta\alpha}\left(\bar{n}_{\text{spp}}(\omega_{0})+1\right)+i\Delta^{\prime}_{\beta\alpha}\right].$
(638)

We can group together common terms noting that we can particularize for our
two qubit system by defining
$\Gamma_{\alpha\beta}=\Gamma_{\beta\alpha}\equiv\Gamma_{12}$ as well as
$\Delta_{\alpha\beta}=\Delta_{\beta\alpha}\equiv g_{12}$ and
$\Delta_{\alpha\beta}^{\prime}=\Delta_{\beta\alpha}^{\prime}\equiv\Delta_{12}^{\prime}$.
The summation yields

$\displaystyle-\left(\sigma_{1}^{+}\sigma_{2}^{-}\hat{\rho}_{S}^{I}(t)-\sigma_{2}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{1}^{+}+\sigma_{2}^{+}\sigma_{1}^{-}\hat{\rho}_{S}^{I}(t)-\sigma_{1}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{2}^{+}\right)\left[\Gamma_{12}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]+i\left(\Delta_{12}+\Delta^{\prime}_{12}\right)\right]$
$\displaystyle-\left(\hat{\rho}_{S}^{I}(t)\sigma_{2}^{\text{+}}\sigma_{1}^{-}-\sigma_{1}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{2}^{+}+\hat{\rho}_{S}^{I}(t)\sigma_{1}^{\text{+}}\sigma_{2}^{-}-\sigma_{2}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{1}^{+}\right)\left[\Gamma_{12}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]-i\left(\Delta_{12}+\Delta^{\prime}_{12}\right)\right]$
$\displaystyle-\left(\sigma_{1}^{-}\sigma_{2}^{+}\hat{\rho}_{S}^{I}(t)-\sigma_{2}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{1}^{-}+\sigma_{2}^{-}\sigma_{1}^{+}\hat{\rho}_{S}^{I}(t)-\sigma_{1}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{2}^{-}\right)\left[\Gamma_{12}\bar{n}_{\text{spp}}(\omega_{0})-i\Delta^{\prime}_{12}\right]$
$\displaystyle-\left(\hat{\rho}_{S}^{I}(t)\sigma_{2}^{-}\sigma_{1}^{+}-\sigma_{1}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{2}^{-}+\hat{\rho}_{S}^{I}(t)\sigma_{1}^{-}\sigma_{2}^{+}-\sigma_{2}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{1}^{-}\right)\left[\Gamma_{12}\bar{n}_{\text{spp}}(\omega_{0})+i\Delta^{\prime}_{12}\right],$
(639)

and what follows is

$\displaystyle\Gamma_{12}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(2\sigma_{2}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{1}^{+}+2\sigma_{1}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{2}^{+}-\left\\{\sigma_{1}^{+}\sigma_{2}^{-},\hat{\rho}_{S}^{I}(t)\right\\}-\left\\{\sigma_{2}^{+}\sigma_{1}^{-},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle+\Gamma_{12}\bar{n}_{\text{spp}}(\omega_{0})\left(2\sigma_{2}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{1}^{-}+2\sigma_{1}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{2}^{-}-\left\\{\sigma_{1}^{-}\sigma_{2}^{+},\hat{\rho}_{S}^{I}(t)\right\\}-\left\\{\sigma_{2}^{-}\sigma_{1}^{+},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle+i\left(g_{12}+\Delta^{\prime}_{12}\right)\left(\left[\sigma_{1}^{+}\sigma_{2}^{-}+\sigma_{2}^{+}\sigma_{1}^{-},\hat{\rho}_{S}^{I}(t)\right]\right)$
$\displaystyle-i\Delta^{\prime}_{12}\left(\left[\sigma_{1}^{-}\sigma_{2}^{+}+\sigma_{2}^{-}\sigma_{1}^{+},\hat{\rho}_{S}^{I}(t)\right]\right),$
(640)

since the rasing and lowering operators act in different qubits we can write
the commutator

$\left[\sigma_{i}^{+},\sigma_{j}^{-}\right]=0,$ (641)

for $i\neq j$, using which we obtain

$\displaystyle\Gamma_{12}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(2\sigma_{2}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{1}^{+}+2\sigma_{1}^{-}\hat{\rho}_{S}^{I}(t)\sigma_{2}^{+}-\left\\{\sigma_{1}^{+}\sigma_{2}^{-},\hat{\rho}_{S}^{I}(t)\right\\}-\left\\{\sigma_{2}^{+}\sigma_{1}^{-},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle+\Gamma_{12}\bar{n}_{\text{spp}}(\omega_{0})\left(2\sigma_{2}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{1}^{-}+2\sigma_{1}^{+}\hat{\rho}_{S}^{I}(t)\sigma_{2}^{-}-\left\\{\sigma_{1}^{-}\sigma_{2}^{+},\hat{\rho}_{S}^{I}(t)\right\\}-\left\\{\sigma_{2}^{-}\sigma_{1}^{+},\hat{\rho}_{S}^{I}(t)\right\\}\right)$
$\displaystyle+ig_{12}\left(\left[\sigma_{1}^{+}\sigma_{2}^{-}+\sigma_{2}^{+}\sigma_{1}^{-},\hat{\rho}_{S}^{I}(t)\right]\right).$
(642)

We can now define additional jump operators

$\begin{cases}\hat{K}_{\alpha,1}=\sqrt{\Gamma_{12}\bar{n}_{\text{spp}}(\omega_{0})}\sigma_{\alpha}^{+},\\\
\hat{K}_{\alpha,2}=\sqrt{\Gamma_{12}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]}\sigma_{\alpha}^{-}.\end{cases}$
(643)

Converting back to the Schrodinger picture, we have, with the previously
defined jump operators

$\displaystyle\sum_{\alpha\neq\beta}\left(2\hat{K}_{\alpha,2}\hat{\rho}_{S}(t)\hat{K}_{\beta,2}^{\dagger}-\left\\{\hat{K}_{\alpha,2}^{\dagger}\hat{K}_{\beta,2},\hat{\rho}_{S}(t)\right\\}\right)$
$\displaystyle+\left(2\hat{K}_{\alpha,1}\hat{\rho}_{S}(t)\hat{K}_{\beta,1}^{\dagger}-\left\\{\hat{K}_{\alpha,1}^{\dagger}\hat{K}_{\beta,1},\hat{\rho}_{S}(t)\right\\}\right)-i(\omega_{0}+g_{12})\left[\sigma_{\alpha}^{+}\sigma_{\beta}^{-},\hat{\rho}_{S}(t)\right],$
(644)

To all intents and purposes we can now write all terms including the coherent
shift $g_{12}$ and the operators $\hat{L}_{\alpha,k}$ and $\hat{K}_{\alpha,k}$
in the master equation, and therefore we have as a final result

$\displaystyle\frac{d\hat{\rho}_{S}}{dt}=$
$\displaystyle-\frac{i}{\hbar}\left[\hbar\omega_{12}\sum_{\alpha\neq\beta}\sigma_{\alpha}^{+}\sigma_{\beta}^{-},\hat{\rho}_{S}(t)\right]-\frac{i}{\hbar}\left[\hbar\omega_{0}^{\prime}\sum_{\alpha}\sigma_{\alpha}^{+}\sigma_{\alpha}^{-},\hat{\rho}_{S}(t)\right]$
$\displaystyle+\sum_{\alpha}\left(\sum_{k}2\hat{L}_{\alpha,k}\hat{\rho}_{S}(t)\hat{L}_{\alpha,k}^{\dagger}-\left\\{\hat{L}_{\alpha,k}^{\dagger}\hat{L}_{\alpha,k},\hat{\rho}_{S}(t)\right\\}\right)$
$\displaystyle+\sum_{\alpha\neq\beta}\left(\sum_{k}2\hat{K}_{\alpha,k}\hat{\rho}_{S}(t)\hat{K}_{\beta,k}^{\dagger}-\left\\{\hat{K}_{\alpha,k}^{\dagger}\hat{K}_{\beta,k},\hat{\rho}_{S}(t)\right\\}\right).$
(645)

with
$\omega_{0}^{\prime}=\omega_{0}+\Delta_{\alpha}+2\Delta_{\alpha}^{\prime}$ and
$\omega_{12}=\omega_{0}+g_{12}$. This brings our equation to Lindblad form. To
bring this into an equation for the probability we can proceed as before,
expanding the density matrix in the respective terms

$\displaystyle\hat{\rho}_{S}$
$\displaystyle=\sum_{ij}\sum_{kl}\rho_{ij,kl}\left|ij\right\rangle\left\langle
kl\right|,$
$\displaystyle=\begin{bmatrix}\rho_{ee,ee}&\rho_{ee,eg}&\rho_{ee,ge}&\rho_{ee,gg}\\\
\rho_{eg,ee}&\rho_{eg,eg}&\rho_{eg,ge}&\rho_{eg,gg}\\\
\rho_{ge,ee}&\rho_{ge,eg}&\rho_{ge,ge}&\rho_{ge,gg}\\\
\rho_{gg,ee}&\rho_{gg,eg}&\rho_{gg,ge}&\rho_{gggg}\end{bmatrix}.$ (646)

We define for simplicity the following super-operators

$\mathbb{I}_{\alpha\beta}\left[\hat{\rho}_{S}(t)\right]=-\frac{i}{\hbar}\left[\hbar\omega_{12}\sigma_{\alpha}^{+}\sigma_{\beta}^{-},\hat{\rho}_{S}(t)\right],$
(647)

corresponding to the term in the dynamical equation that gives the interaction
between the two-level systems

$\mathbb{C}_{\alpha}\left[\hat{\rho}_{S}(t)\right]=-\frac{i}{\hbar}\left[\hbar\omega_{0}^{\prime}\sigma_{\alpha}^{+}\sigma_{\alpha}^{-},\hat{\rho}_{S}(t)\right],$
(648)

corresponding to the term that gives the coherent evolution of each qubit,

$\mathbb{L}_{\alpha}\left[\hat{\rho}_{S}(t)\right]=\sum_{k}2\hat{L}_{\alpha,k}\hat{\rho}_{S}(t)\hat{L}_{\alpha,k}^{\dagger}-\left\\{\hat{L}_{\alpha,k}^{\dagger}\hat{L}_{\alpha,k},\hat{\rho}_{S}(t)\right\\},$
(649)

which gives the incoherent evolution of both qubits due to the exchange of
plasmons with themselves and finally

$\mathbb{K}_{\alpha\beta}\left[\hat{\rho}_{S}(t)\right]=\sum_{k}2\hat{K}_{\alpha,k}\hat{\rho}_{S}(t)\hat{K}_{\beta,k}^{\dagger}-\left\\{\hat{K}_{\alpha,k}^{\dagger}\hat{K}_{\beta,k},\hat{\rho}_{S}(t)\right\\},$
(650)

which gives the incoherent evolution of both qubits due to the exchange of
plasmons between eachother. The information regarding this exchange is
contained in the jump operators $K_{\alpha,k}$. We have therefore

$\frac{d\hat{\rho}_{S}}{dt}=\left[\sum_{\alpha}\mathbb{C}_{\alpha}+\sum_{\alpha\neq\beta}\mathbb{I}_{\alpha\beta}+\sum_{\alpha}\mathbb{L}_{\alpha}+\sum_{\alpha\neq\beta}\mathbb{K}_{\alpha\beta}\right]\left[\hat{\rho}_{S}(t)\right].$
(651)

We make these definitions in order to simplify the following calculation. Note
that we can write the jump operators in matrix form since we have

$\sigma_{1}^{+}=\sigma^{+}\otimes\mathds{1}=\begin{bmatrix}0&0&1&0\\\
0&0&0&1\\\ 0&0&0&0\\\ 0&0&0&0\end{bmatrix}.$ (652)

$\sigma_{2}^{+}=\mathds{1}\otimes\sigma^{+}=\begin{bmatrix}0&1&0&0\\\
0&0&0&0\\\ 0&0&0&1\\\ 0&0&0&0\end{bmatrix}.$ (653)

We write each term of the dynamical equation in matrix form using these
operators. We are left with

$\sum_{\alpha}\mathbb{C}_{\alpha}=i\omega_{0}^{\prime}\begin{bmatrix}0&-\rho_{eg,ee}&-\rho_{ge,ee}&-2\rho_{gg,ee}\\\
\rho_{eg,ee}&0&0&-\rho_{gg,eg}\\\ \rho_{ge,ee}&0&0&-\rho_{gg,ge}\\\
2\rho_{gg,ee}&\rho_{gg,eg}&\rho_{gg,ge}&0\end{bmatrix},$ (654)

$\sum_{\alpha\neq\beta}\mathbb{I}_{\alpha\beta}=i\omega_{12}\begin{bmatrix}0&\rho_{eg,ee}&\rho_{eg,ee}&0\\\
-\rho_{ge,ee}&0&\rho_{eg,eg}-\rho_{ge,ge}&-\rho_{gg,ge}\\\
-\rho_{eg,ee}&\rho_{ge,ge}-\rho_{eg,eg}&0&-\rho_{gg,eg}\\\
0&\rho_{gg,ge}&\rho_{gg,eg}&0\end{bmatrix}.$ (655)

The whole matrix can easily be calculated for the remaining terms as well
$\sum_{\alpha}\mathbb{L}_{\alpha}$ and
$\sum_{\alpha\neq\beta}\mathbb{K}_{\alpha\beta}$, however, since we are
intrested only in the populations, that is to say, the diagonal elements of
the density matrix, we find the neither $\mathbb{C}_{\alpha}$ nor
$\mathbb{I}_{\alpha\beta}$ contribute, since they both have zeros along the
diagonals. The diagonal contributions of $\mathbb{L_{\alpha}}$ are, if
organized in a vector

$\text{diag}\left[\sum_{\alpha}\mathbb{L}_{\alpha}\right]=\gamma\begin{bmatrix}-2\rho_{ee,ee}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]+\bar{n}_{\text{spp}}(\omega_{0})(\rho_{eg,eg}+\rho_{ge,ge})\\\
\rho_{ee,ee}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]-\rho_{eg,eg}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]-\bar{n}_{\text{spp}}(\omega_{0})\rho_{eg,eg}+\bar{n}_{\text{spp}}(\omega_{0})\rho_{gg,gg}\\\
\rho_{ee,ee}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]-\rho_{ge,ge}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]-\bar{n}_{\text{spp}}(\omega_{0})\rho_{ge,ge}+\bar{n}_{\text{spp}}(\omega_{0})\rho_{gg,gg}\\\
-2\rho_{gg,gg}\bar{n}_{\text{spp}}(\omega_{0})+\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right](\rho_{eg,eg}+\rho_{ge,ge})\end{bmatrix},$
(656)

and the contributions from $\mathbb{K}_{\alpha\beta}$ read

$\text{diag}\left[\sum_{\alpha\neq\beta}\mathbb{K}_{\alpha\beta}\right]=\gamma_{12}\begin{bmatrix}2\bar{n}_{\text{spp}}(\omega_{0})\rho_{ge,eg}\\\
-\bar{n}_{\text{spp}}(\omega_{0})\rho_{ge,eg}-\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{ge,eg}\\\
-\bar{n}_{\text{spp}}(\omega_{0})\rho_{ge,eg}-\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{ge,eg}\\\
2\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{ge,eg}\end{bmatrix}.$
(657)

As expected, the coherent evolutions do not contribute to the populations
since they are responsible at most for a phase shift, however in evaluating
the evolution of the diagonal terms of the density matrix we also find a
dependance on off-diagonal terms. This suggests that a change of basis is in
order, and in fact, we can look at the basis we already used for the
Schrodinger equation approach. In this case, we have the states
$\left|3\right\rangle=\left|ee\right\rangle,\left|+\right\rangle=\left(\left|eg\right\rangle+\left|ge\right\rangle\right)/\sqrt{2},\left|-\right\rangle=\left(\left|eg\right\rangle-\left|ge\right\rangle\right)/\sqrt{2}$
and $\left|0\right\rangle=\left|gg\right\rangle$. We can construct a change of
basis matrix by writing

$U=\begin{bmatrix}1&0&0&0\\\ 0&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\\
0&\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}&0\\\ 0&0&0&1\end{bmatrix},$ (658)

and construct the raising and lowering operators in this basis by writing

$\sigma^{\prime}{}_{\alpha}^{\pm}=U^{\dagger}\sigma_{\alpha}^{\pm}U.$ (659)

We find, for instance

$\sigma^{\prime}{}_{1}^{+}=\begin{bmatrix}0&\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}&0\\\
0&0&0&\frac{1}{\sqrt{2}}\\\ 0&0&0&\frac{1}{\sqrt{2}}\\\ 0&0&0&0\end{bmatrix},$
(660)

$\sigma^{\prime}{}_{2}^{+}=\begin{bmatrix}0&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\\
0&0&0&\frac{1}{\sqrt{2}}\\\ 0&0&0&-\frac{1}{\sqrt{2}}\\\
0&0&0&0\end{bmatrix}.$ (661)

In this basis the density matrix also takes the shape

$\hat{\rho}_{S}=\begin{bmatrix}\rho_{3,3}&\rho_{+,+}&\rho_{+,-}&\rho_{3,0}\\\
\rho_{+,3}&\rho_{+,+}&\rho_{+,-}&\rho_{+,0}\\\
\rho_{-,3}&\rho_{-,+}&\rho_{-,-}&\rho_{-,0}\\\
\rho_{0,3}&\rho_{0,+}&\rho_{0,-}&\rho_{0,0}\end{bmatrix},$ (662)

which means that we can calculate the several components of the dynamical
equation, speciffically, the populations by taking the diagonal elements. As
before, there is no contribution from the coherent evolution given by
$\mathbb{C}_{\alpha}$ or $\mathbb{I}_{\alpha\beta}$, however we can look at
the diagonal elements of the remaining parts

$\text{diag}\left[\sum_{\alpha}\mathbb{L}_{\alpha}\right]=\gamma\begin{bmatrix}-2\rho_{3,3}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]-\bar{n}_{\text{spp}}(\omega_{0})\rho_{3,3}+\bar{n}_{\text{spp}}(\omega_{0})\left(\rho_{-,-}+\rho_{+,+}\right)\\\
\rho_{3,3}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]+\bar{n}_{\text{spp}}(\omega_{0})\rho_{0,0}-\bar{n}_{\text{spp}}(\omega_{0})\rho_{+.+}-\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{+,+}\\\
\rho_{3,3}\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]+\bar{n}_{\text{spp}}(\omega_{0})\rho_{0,0}-\bar{n}_{\text{spp}}(\omega_{0})\rho_{-,-}-\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{-,-}\\\
-2\rho_{0,0}\bar{n}_{\text{spp}}(\omega_{0})+\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(\rho_{-,-}+\rho_{+,+}\right)\end{bmatrix},$
(663)

$\text{diag}\left[\sum_{\alpha}\mathbb{K}_{\alpha}\right]=\gamma_{12}\begin{bmatrix}\bar{n}_{\text{spp}}(\omega_{0})\left(\rho_{+,+}-\rho_{-,-}\right)\\\
\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{3,3}+\bar{n}_{\text{spp}}(\omega_{0})\rho_{0,0}-\bar{n}_{\text{spp}}(\omega_{0})\rho_{+,+}-\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{+,+}\\\
-\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{3,3}+\bar{n}_{\text{spp}}(\omega_{0})\rho_{0,0}+\bar{n}_{\text{spp}}(\omega_{0})\rho_{-,-}+\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\rho_{-,-}\\\
-\left[\bar{n}_{\text{spp}}(\omega_{0})+1\right]\left(\rho_{-,-}-\rho_{+,+}\right)\end{bmatrix}.$
(664)

We see that the diagonal elements of the density matrix in this basis are
self-contained, that is to say, they depend only on each other, as opposed to
the natural basis we analyzed previously. This allows us to evaluate the
dynamics of the density matrix looking only at the diagonal elements. Solving
the dynamical equation can be done numerically, and we obtain the figures
presented in the main text. Note that if we let the temperature go to 0, which
means $\bar{n}_{\text{spp}}(\omega_{0})\to 0$, we find

$\text{diag}\left[\sum_{\alpha}\mathbb{L}_{\alpha}\right]=\gamma\begin{bmatrix}-2\rho_{3,3}\\\
\rho_{3,3}-\rho_{+,+}\\\ \rho_{3,3}-\rho_{-,-}\\\
\rho_{-,-}+\rho_{+,+}\end{bmatrix},$ (665)

$\text{diag}\left[\sum_{\alpha}\mathbb{K}_{\alpha}\right]=\gamma_{12}\begin{bmatrix}0\\\
\rho_{3,3}-\rho_{+,+}\\\ -\rho_{3,3}+\rho_{-,-}\\\
\rho_{+,+}-\rho_{-,-}\end{bmatrix},$ (666)

and therefore, writing the populations as $\rho_{3,3}\equiv P_{3}$,
$\rho_{\pm,\pm}\equiv P_{\pm}$ and $\rho_{0,0}\equiv P_{0}$, we find the
dinamical equation

$\begin{bmatrix}\dot{P}_{3}\\\ \dot{P}_{\text{+}}\\\ \dot{P}_{-}\\\
\dot{P}_{0}\end{bmatrix}=\gamma\begin{bmatrix}-2P_{3}\\\ P_{3}-P_{+}\\\
P_{3}-P_{-}\\\ P_{-}+P_{+}\end{bmatrix}+\gamma_{12}\begin{bmatrix}0\\\
P_{3}-P_{+}\\\ -P_{3}+P_{-}\\\ P_{+}-P_{-}\end{bmatrix},$ (667)

which paired with the initial conditions $P_{3}(0)=0$, $P_{+}(0)=1/2$,
$P_{-}(0)=1/2$, $P_{0}(0)=0$ gives the dynamics

$\begin{cases}P_{3}(t)=0,\\\
P_{+}(t)=\frac{1}{2}e^{-(\gamma+\gamma_{12})t},\\\
P_{-}(t)=\frac{1}{2}e^{-(\gamma-\gamma_{12})t},\\\
P_{0}(t)=1-\frac{e^{-(\gamma-\gamma_{12})t}+e^{-(\gamma+\gamma_{12})t}}{2}.\end{cases}$
(668)

which matches exactly those obtained via the Schrodinger equation.

## References

  * Gonçalves and Peres [2016] P. A. D. Gonçalves and N. M. R. Peres,  _An Introduction to Graphene Plasmonics_ (World Scientific, 2016).
  * Gonzalez-Tudela  _et al._ [2011] A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno, L. Martin-Moreno, C. Tejedor, and F. J. Garcia-Vidal, Phys. Rev. Lett. 106, 020501 (2011).
  * Nechepurenko  _et al._ [2018] I. A. Nechepurenko, E. S. Andrianov, A. A. Zyablovsky, A. V. Dorofeenko, A. A. Pukhov, and Y. E. Lozovik, Phys. Rev. B 98, 075411 (2018).
  * Manzano [2020] D. Manzano, AIP Advances 10, 025106 (2020).
  * Hohenester [2020] U. Hohenester,  _Nano and Quantum Optics: An Introduction to Basic Principles and Theory_ (Spinger, 2020).
  * Carmichael and de Bruxelles [1993] H. Carmichael and U. de Bruxelles,  _An Open Systems Approach to Quantum Optics: Lectures Presented at the Université Libre de Bruxelles, October 28 to November 4, 1991_ , vol. 18 (Springer, 1993).
  * Carmichael [2013] H. Carmichael,  _Statistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations_ , Theoretical and Mathematical Physics (Springer Berlin Heidelberg, 2013).
  * Hecht [2000] K. Hecht, Natural line width: Wigner-weisskopf treatment (2000) pp. 645–654.
  * Scheel and Buhmann [2009] S. Scheel and S. Y. Buhmann, Macroscopic qed - concepts and applications (2009), arXiv:0902.3586 [quant-ph] .
  * Buhmann [2013] S. Buhmann,  _Dispersion Forces I: Macroscopic Quantum Electrodynamics and Ground-State Casimir, Casimir-Polder and van der Waals Forces_ , Springer Tracts in Modern Physics (Springer Berlin Heidelberg, 2013).