# Solving a Lindblad differential equation in matrix form

I am trying to solve the following Lindblad equation analytically with Mathematica:

$$\frac{d\rho}{dt}=-i[H,\rho]+\sum\nolimits_{i=0}^1 K_i\rho K_i^+-\frac12\cdot\{K_i^+K_i,\rho\},$$ where $$H=\begin{pmatrix} -1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & -1 \end{pmatrix},$$

$$K_1=\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \end{pmatrix},$$

$$K_2=\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & -1 \end{pmatrix}.$$

My code looks like this:

Clear["Global*"]
H[t_] = {{-1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, -1}};
K1 = {{1/Sqrt[2], 0, 0, 0}, {0, 1/Sqrt[2], 0, 0}, {0, 0, -1/Sqrt[2],
0}, {0, 0, 0, -1/Sqrt[2]}};
K1plus = ConjugateTranspose[K1];
K2 = {{1/Sqrt[2], 0, 0, 0}, {0, -1/Sqrt[2], 0, 0}, {0, 0, 1/Sqrt[2],
0}, {0, 0, 0, -1/Sqrt[2]}};
K2plus = ConjugateTranspose[K2];
eqn = I*D[rho[t], t] == -I*H[t] . rho[t] + I*rho[t] . H[t] +
K1 . rho[t] . K1plus - 1/2*K1plus . K1 . rho[t] -
1/2*rho[t] . K1plus . K1 + K2 . rho[t] . K2plus -
1/2*K2plus . K2 . rho[t] - 1/2*rho[t] . K2plus . K2;

sol = DSolve[eqn, rho, t]


Nevertheless, I don't get the exact solution. Mathematica only prints the input DSolve[...].

Does that mean that it can't find the solution?

DSolve can't automatically determine the dimensionality of your solution. Instead it assumes it to be scalar. To solve it as a matrix, you have to provide it as a matrix $$\rho(t)=\begin{pmatrix}\rho_{1,1}(t) & \rho_{1,2}(t) & \dots\\\rho_{2,1}(t) & \rho_{2,2}(t) & \dots \\\vdots & \vdots & \ddots\end{pmatrix}$$

You can do this easily with Table:

rhoMat = Table[rho[i][j][t], {i, 1, Length@K1}, {j, 1, Length@K1}];


The replace every occurence of rho[t] with rhoMat:

eqn = I*D[rhoMat, t] == -I*H[t] . rhoMat + I*rhoMat . H[t] +
K1 . rhoMat . K1plus - 1/2*K1plus . K1 . rhoMat -
1/2*rhoMat . K1plus . K1 + K2 . rhoMat . K2plus -
1/2*K2plus . K2 . rhoMat - 1/2*rhoMat . K2plus . K2;


Plug this into DSolve and get your solutions:

DSolve[eqn // Flatten, rhoMat // Flatten, t];
First[rhoMat /. DSolve[eqn // Flatten, rhoMat // Flatten, t]]


{{C[1], E^((2 + I) t) C[2], E^((2 + I) t) C[3], E^(2 I t) C[4]}, {E^((-2 + I) t) C[5], C[6], E^(2 I t) C[7], E^((-2 + I) t) C[8]}, {E^((-2 + I) t) C[9], E^(2 I t) C[10], C[11], E^((-2 + I) t) C[12]}, {E^(2 I t) C[13], E^((2 + I) t) C[14], E^((2 + I) t) C[15], C[16]}}

Btw, I think you added an unnecessary I infront of the time derivative. The original Lindbladian does not feature this.

EDIT: When I solve the Lindbladian, I usually do not take DSolve. Instead I reduce the problem to $$\rho'(t)=M\rho(t)$$ Which can easily be done with CoefficientArrays and then solved as a Eigensystem problem. I paste the code for your equation here (ommiting the I whih I think is wrong).

eqn = -I*H[t] . rhoMat + I*rhoMat . H[t] + K1 . rhoMat . K1plus -
1/2*K1plus . K1 . rhoMat - 1/2*rhoMat . K1plus . K1 +
K2 . rhoMat . K2plus - 1/2*K2plus . K2 . rhoMat -
1/2*rhoMat . K2plus . K2;

systemMatrix = Last[CoefficientArrays[Flatten@eqn, Flatten@rhoMat]];

{eigVec, eigSys} = Eigensystem[mat2];
sol = Total[
Table[(Exp[eigVec[[i]]*t]*eigSys[[i]]*C[i]), {i, 1,
Length@eigVec}]];
sol = rhoMat /. Thread[Flatten[rhoMat] -> sol]
`

{{C[16], E^((-1 + 2 I) t) C[4], E^((-1 + 2 I) t) C[3], E^(-2 t) C[12]}, {E^((-1 - 2 I) t) C[8], C[15], E^(-2 t) C[11], E^((-1 - 2 I) t) C[7]}, {E^((-1 - 2 I) t) C[6], E^(-2 t) C[10], C[ 14], E^((-1

• 2 I) t) C[5]}, {E^(-2 t) C[9], E^((-1 + 2 I) t) C[2], E^((-1 + 2 I) t) C[1], C[13]}}

Besides the enumeration of the constants, the solutions are identical. With this approach, one can also easily solve for initial value problems with $$\rho(0)=m_{eig}^TC$$

For which $$C$$ is your sought after vector of constants, $$m_{eig}$$ the eigensystem matrix and $$\rho(0)$$ your intial boundary values.