# The Dissipator of a Two-Level-System as an Explicit Superoperator¶

We derive the analytical dissipation superoperator for decay and dephasing in the two-level system

In :
from qnet.algebra.operator_algebra import *
from qnet.algebra.hilbert_space_algebra import *
from qnet.algebra.state_algebra import *
from qnet.convert.to_sympy_matrix import convert_to_sympy_matrix
import sympy
import numpy as np
from sympy import sqrt, symbols

In :
sympy.init_printing()


We define our symbols as follows:

In :
hs = LocalSpace('TLS', basis=(0,1))

In :
γ1, γ2 = symbols('gamma_1 gamma_2', positive=True)
ρ00, ρ11 = symbols('rho_00, rho_11', positive=True)
ρ01, ρ10 = symbols('rho_01, rho_10')
half = sympy.S(1) / 2

In :
L1 = sqrt(γ1) * LocalSigma(0, 1, hs=hs)
L2 = sqrt(2*γ2) * LocalSigma(1, 1, hs=hs)
L = [L1, L2]

In :
ket0, ket1 = BasisKet(0, hs=hs), BasisKet(1, hs=hs)
bra0, bra1 = Bra(ket0), Bra(ket1)

In :
def ketbra(i, j):
k = {0: ket0, 1: ket1}
return k[i] * Bra(k[j])

In :
basis_rho = [ # concatenation of columns!
ketbra(0, 0),
ketbra(1, 0),
ketbra(0, 1),
ketbra(1, 1)
]
rho = ρ00 * ketbra(0, 0) + ρ01 * ketbra(0, 1) + ρ10 * ketbra(1, 0) + ρ11 * ketbra(1, 1)


The dissipator in Lindblad form is defined as

In :
def D(Ls, rho):
res = 0
for L in Ls:
L_dag_L = L.dag() * L
res += (L * rho * L.dag() - half * (L_dag_L * rho + rho * L_dag_L)).expand()
return res


Applied to a general density matrix this is:

In :
convert_to_sympy_matrix(D(L, rho))

Out:
$$\left[\begin{matrix}\gamma_{1} \rho_{11} & - \frac{\gamma_{1} \rho_{01}}{2} - \gamma_{2} \rho_{01}\\- \frac{\gamma_{1} \rho_{10}}{2} - \gamma_{2} \rho_{10} & - \gamma_{1} \rho_{11}\end{matrix}\right]$$

We now calculate the explicit form of the dissipator as a superoperator

In :
def hs_prod(A, B):
"""Hilbert-Schmidt product"""
return tr((A.dag() * B).expand(), over_space=hs)

In :
def get_DSup():
DSup = sympy.Matrix(np.zeros((4,4)))
for i in range(4):
for j in range(4):
d = hs_prod(basis_rho[i], D(L, basis_rho[j]))
if d != 0:
assert d.term == IdentityOperator
DSup[i, j] = d.coeff
else:
DSup[i, j] = sympy.S(0)
return DSup

In :
DSup = get_DSup()


We now have the desired superoperator

In :
DSup

Out:
$$\left[\begin{matrix}0 & 0 & 0 & \gamma_{1}\\0 & - \frac{\gamma_{1}}{2} - \gamma_{2} & 0 & 0\\0 & 0 & - \frac{\gamma_{1}}{2} - \gamma_{2} & 0\\0 & 0 & 0 & - \gamma_{1}\end{matrix}\right]$$

Let's test this by applying the above matrix to a vectorized form of a general density matrix (concatenation of columns):

In :
rho_vec = sympy.Matrix([
[ρ00,],
[ρ10,],
[ρ01,],
[ρ11,]])

In :
(DSup * rho_vec).expand()

Out:
$$\left[\begin{matrix}\gamma_{1} \rho_{11}\\- \frac{\gamma_{1} \rho_{10}}{2} - \gamma_{2} \rho_{10}\\- \frac{\gamma_{1} \rho_{01}}{2} - \gamma_{2} \rho_{01}\\- \gamma_{1} \rho_{11}\end{matrix}\right]$$

This is equal to the vectorized form of the previous result

In :
convert_to_sympy_matrix(D(L, rho)).transpose().reshape(4,1)

Out:
$$\left[\begin{matrix}\gamma_{1} \rho_{11}\\- \frac{\gamma_{1} \rho_{10}}{2} - \gamma_{2} \rho_{10}\\- \frac{\gamma_{1} \rho_{01}}{2} - \gamma_{2} \rho_{01}\\- \gamma_{1} \rho_{11}\end{matrix}\right]$$

We can also use this to solve the equation of motion:

In :
t = symbols('t', positive=True)

In :
(DSup * t).exp() * rho_vec

Out:
$$\left[\begin{matrix}\rho_{00} + \rho_{11} \left(1 - e^{- \gamma_{1} t}\right)\\\frac{\rho_{10}}{e^{\frac{\gamma_{1} t}{2} + \gamma_{2} t}}\\\frac{\rho_{01}}{e^{\frac{\gamma_{1} t}{2} + \gamma_{2} t}}\\\frac{\rho_{11}}{e^{\gamma_{1} t}}\end{matrix}\right]$$