Superoperators and Vectorization

Vectorization turns a linear equation for a matrix into an ordinary linear system on a longer vector. Its benefit is access to sparse eigensolvers and exponential actions; its danger is that a transpose error can preserve dimensions while changing the physics. Here columns are stacked, and all tensor factors follow that choice.

Helpful background. Density matrices in code supplies state diagnostics; tensor products explains Kronecker ordering.

Column stacking and multiplication

For a $d\times d$ matrix, define $\operatorname{vec}(X)_{i+dj}=X_{ij}$ with zero-based indices. Then

$$\operatorname{vec}(AXB)=(B^{\mathsf T}\otimes A)\operatorname{vec}(X).$$

To derive the identity, the coefficient of $X_{kl}$ in $(AXB)_{ij}$ is $A_{ik}B_{lj}$. The Kronecker entry at row $i+dj$, column $k+dl$ is precisely $B^{\mathsf T}_{jl}A_{ik}$. Right multiplication therefore needs a transpose, not an adjoint.

For row stacking the factors reverse. Both conventions are valid; mixing them is not. In NumPy, column stacking requires reshape/ravel with order=“F”, independently of how the original array was allocated.

This stacking order is different from the default row-major tensor reshape used to label physical subsystems. Keep the two maps explicit: a vectorized matrix index $i+dj$ labels the pair $(j,i)$ in a Kronecker product with the first factor varying slowly. It does not introduce an extra physical subsystem. After applying a superoperator, reconstruct the matrix using the same column-stacking rule. A round-trip test alone is insufficient: two consistently wrong reshapes can cancel, so also compare the map’s action with direct matrix multiplication.

Hamiltonian and dissipative generators

For Hermitian $H$ and collapse operators $L_\mu$ containing their square-root rates, the column-stacked Lindblad generator is

$$\mathsf L={}-\frac{i}{\hbar}(I\otimes H-H^{\mathsf T}\otimes I) +\sum_\mu\left[ L_\mu^*\otimes L_\mu -\frac12 I\otimes L_\mu^\dagger L_\mu -\frac12(L_\mu^\dagger L_\mu)^{\mathsf T}\otimes I \right].$$

This follows by applying the multiplication identity term by term. The conceptual assumptions behind this generator belong to the Lindblad–GKSL equation. Constructing a matrix does not validate those assumptions.

The trace functional is the row $\operatorname{vec}(I)^\dagger$. Trace preservation means

$$\operatorname{vec}(I)^\dagger\mathsf L=0.$$

It does not generally imply $\mathsf L\operatorname{vec}(I)=0$; that stronger property is unitality. The distinction catches a frequent erroneous “steady state” test.

A two-level convention test

Take $H=\hbar\omega|1\rangle\langle1|$ and $L=\sqrt{\gamma}|0\rangle\langle1|$. In the vector order $(\rho_{00},\rho_{10},\rho_{01},\rho_{11})$,

$$\mathsf L= \begin{pmatrix} 0&0&0&\gamma\\ 0&-\gamma/2-i\omega&0&0\\ 0&0&-\gamma/2+i\omega&0\\ 0&0&0&-\gamma \end{pmatrix}.$$

The population lost from $\rho_{11}$ reappears in $\rho_{00}$; conjugate coherences rotate in opposite directions. Testing only diagonal states misses the frequency-sign error.

For a stronger general test, apply both the direct matrix formula and $\mathsf L$ to every matrix unit $E_{ij}$ in a small dimension. These units need not be physical states: they form a basis for the entire linear map. Also test a complex Hermitian matrix to verify Hermiticity preservation.

Storage and solver choices

A density vector has $d^2$ entries; a dense superoperator has $d^4$. For local sparse $H$ and $L_\mu$, assemble sparse Kronecker products or implement the action on reshaped matrices without storing $\mathsf L$. Matrix-free action saves storage but still requires a $d^2$ state and careful work-buffer management.

A Liouvillian is generally non-Hermitian and can be nonnormal. Do not use a Hermitian eigensolver merely because its input density matrices are Hermitian. Near-zero eigenvalues require residual checks, and decay rates inferred from eigenvalues can fail to describe short transient amplification of perturbations.

Specify the norm in this last diagnostic. Vector Euclidean norm corresponds to the Hilbert–Schmidt norm of the reshaped operator; its induced propagator norm need not decrease for a nonunital channel. By contrast, a completely positive trace-preserving propagator contracts the trace norm of the Hermitian difference of two physical states. Apparent growth in a different norm is not by itself a failure of the channel, while growing trace distance between two states propagated by the same fixed channel signals a problem.

Exercise

Why is the transpose map positive but not completely positive, and what is its column-stacked matrix?

Solution

Its matrix is the permutation $K$ satisfying $K\operatorname{vec}(X)=\operatorname{vec}(X^{\mathsf T})$, so $K_{i+dj,j+di}=1$. Transposing a density matrix preserves its nonnegative spectrum. Acting on only one half of a Bell density matrix gives an eigenvalue $-1/2$, however. This is a useful test that a superoperator need not represent a physical channel merely because it maps single-system test states to positive matrices.

References

• QuTiP developers, “Superoperators and Vectorized Operators”, QuTiP documentation — column stacking and superoperator representations.
• Watrous, J., The Theory of Quantum Information, Cambridge University Press (2018), author’s text — linear maps, vectorization, and complete positivity.