Solution (source code)

= Solution

Write $\rho=\sum_n r_n\sigma_n$ in the <Generalized Bloch representation>. Linearity of the <Lindbladian> gives $\dot r_m=\sum_n\operatorname{Tr}(\sigma_m\mathcal L(\sigma_n))r_n$. For the Hamiltonian term, cyclicity yields
$$
\operatorname{Tr}(\sigma_m[-iH,\sigma_n])=-i\operatorname{Tr}(H\sigma_n\sigma_m-H\sigma_m\sigma_n)=\operatorname{Tr}(iH[\sigma_m,\sigma_n])=L_{mn}.
$$
For one <Lindblad operator> $V_d$, the jump term becomes $\operatorname{Tr}(\sigma_mV_d\sigma_nV_d^\dagger)=\operatorname{Tr}(V_d^\dagger\sigma_mV_d\sigma_n)$, and the two remaining terms give
$$
\operatorname{Tr}(\sigma_mV_d^\dagger V_d\sigma_n)+\operatorname{Tr}(\sigma_m\sigma_nV_d^\dagger V_d)=\operatorname{Tr}(V_d^\dagger V_d\{\sigma_m,\sigma_n\}).
$$
Thus
$$
\boxed{D^{(d)}_{mn}=\operatorname{Tr}(V_d^\dagger\sigma_mV_d\sigma_n)-\frac12\operatorname{Tr}(V_d^\dagger V_d\{\sigma_m,\sigma_n\}),\qquad \dot r=\left(L+\sum_dD^{(d)}\right)r}.
$$
Every coefficient is real because $\mathcal L(\sigma_n)$ is Hermitian. In addition, $L_{mn}=-L_{nm}$, although a dissipative $D^{(d)}$ need not be symmetric.

The Hamiltonian <commutator> has zero <trace>, and each dissipator obeys $\operatorname{Tr}(V\rho V^\dagger)=\operatorname{Tr}(V^\dagger V\rho)$, so $\operatorname{Tr}\mathcal L(\rho)=0$. Consequently $\boxed{\dot r_{N^2}=0,\qquad r_{N^2}=1/\sqrt N}$ for normalized <density operators>. Orthogonality to $I/\sqrt N$ makes every other $\sigma_k$ traceless. Define $s=(r_1,\ldots,r_{N^2-1})^T$. With $M=L+\sum_dD^{(d)}$, <trace> preservation gives the block form
$$
M=\begin{pmatrix}A&b\\0&0\end{pmatrix},\qquad \boxed{\dot s=As+c,\quad A_{mn}=M_{mn},\quad c_m=M_{m,N^2}/\sqrt N}.
$$
Here the indices of $A$ range from $1$ to $N^2-1$. Since $\mathcal L(I)=\sum_d(V_dV_d^\dagger-V_d^\dagger V_d)$, the offset can also be written $c_m=N^{-1}\operatorname{Tr}(\sigma_m\sum_d[V_d,V_d^\dagger])$. This is the <Affine Bloch equation>; its offset is zero for unital dynamics, which preserve the identity.