Skip to content

Spectrum of the Linblad Super Operator

In the study of Markovian open quantum systems, the time-evolution of the system density matrix can, in many interesting scenarios, be given by a master equation in Lindblad form:

\[ \frac{d \rho}{dt} = -i[H,\rho] + \sum_{\mu} \left( L_{\mu} \rho L^\dagger_{\mu} - \frac{1}{2} \{L^\dagger_\mu L_\mu, \rho \} \right), \]

where \(\rho\) is the system density matrix, \(L_\mu\) is a quantum jump operator representing a dissipative process, and \(H\) is the system's Hamiltonian, responsible for generating the coherent part of the dynamics.

The Lindblad spectrum reveals the system's excitations and decay timescales. A key feature is the dissipative gap (the absolute real part of the eigenvalue closest to zero) which determines the slowest relaxation time and signals dissipative phase transitions [1,2]. In this example, we use XDiag to compute the full spectrum of the Linblad superoperator corresponding to a Heisenberg XXZ chain with bulk dephasing.

The system evolves under the Hamiltonian \(\mathcal{H}\) and dephasing jump operators \(L_{n}\):

\[ \mathcal{H} = J \sum_{n} \left( S^+_n S^-_{n+1} + S^+_{n+1} S^-_{n} + \Delta S^z_n S^z_{n+1} \right), \quad L_{n} = \sqrt{\gamma} S^z_n, \]

where \(J\) is the exchange constant and \(\gamma\) is the dissipation strength.

To diagonalize this in XDiag, we first write the Lindblad superoperator by vectorizing the density matrices using the Choi isomorphism:

\[ \rho = \sum_{i,j} \rho_{i,j} |i\rangle \langle j| \rightarrow |\rho\rangle\rangle = \sum_{i,j} \rho_{i,j} |i\rangle \otimes |j\rangle. \]

This introduces an auxiliary spin chain \(\tilde{S}\) acting on the "bra" degrees of freedom. The resulting non-Hermitian operator is:

\[ \mathcal{L} =\sum_{n=0}^{L-1} \left\{ -iJ \left( S_{n}^{+}S_{n+1}^{-} + S_{n+1}^{+}S_{n}^{-} + \Delta S_{n}^{z}S_{n+1}^{z} \right) + iJ \left( \tilde{S}_{n}^{+}\tilde{S}_{n+1}^{-} + \tilde{S}_{n+1}^{+}\tilde{S}_{n}^{-} + \Delta\tilde{S}_{n}^{z}\tilde{S}_{n+1}^{z} \right) + \gamma S_{n}^{z}\tilde{S}_{n}^{z} \right\} -i\frac{\gamma}{4}L. \]

The former operator can be constructed in XDiag as a single OpSum (note that XDiag automatically applies Hermitian conjugation when using certain built-in operators).

Image title

As the system possesses strong U(1) symmetry (\([S^z, \mathcal{H}] = [S^z, L_\mu] = 0\)), we can diagonalize the Liouvillian within sectors of well-defined magnetization and obtain the spectrum and the dissipative gap. The figure above presente these results for \(\Delta=0\): \(a)\) Linblad spectrum for \(S^z=0\), \(\gamma =2.0J\) and \(L=6\), \(b)\) the dissipative gap in the one-body sector as a function of system size for \(\gamma=2.0J\), compared with the analytical result [3] (valid for sizes greater than the critical size); and \(c)\) the dissipative gap in the \(S^z=0\) sector as function of the disspation strength for \(L=6\).

#include <xdiag/all.hpp>
using namespace xdiag;
int main()
try
{
    // system parameters:
    int L = 6; // total number of sites
    double gamma = 2.0; // dissipation strength
    int Nsites_liouville_space = 2 * L; // the effective Hilbert Space is doubled!
    int nup = Nsites_liouville_space / 2; // S^z = 0 sector

    // build OpSum
    arma::cx_mat sx = arma::cx_mat({{0, std::complex<double>(0.5,0.0)},{std::complex<double>(0.5,0.0), 0}});
    arma::cx_mat sy = arma::cx_mat({{0, std::complex<double>(0.0,-0.5)},{std::complex<double>(0.0,0.5), 0}});

    arma::cx_mat sx_sx = arma::kron(sx,sx);
    arma::cx_mat sy_sy = arma::kron(sy,sy);

    //
    Spinhalf block(Nsites_liouville_space, nup);

    OpSum ops;
    for (int i = 0; i < L; ++i) {
        ops += std::complex<double>(0.0,-1.0)* Op("Matrix", {i, (i + 1) % L},sx_sx);
        ops += std::complex<double>(0.0,-1.0)* Op("Matrix", {i, (i + 1) % L},sy_sy);

        ops += std::complex<double>(0.0,1.0)* Op("Matrix", {L+i,L + (i + 1) % L},sx_sx);
        ops += std::complex<double>(0.0,1.0)* Op("Matrix", {L+i,L + (i + 1) % L},sy_sy);
        ops += gamma* Op("SzSz", {i, L + i});
    }

    // convert to matrix to do Full ED:
    arma::cx_vec eigs;

    arma::cx_mat Lmat = matrixC(ops, block); // convert to matrix
    arma::eig_gen(eigs, Lmat);       // perform exact diagonalization

    eigs -= L*gamma/4.0; // correct zero point value

    auto save_fl = FileH5("bulk_dephasing_spectrum.h5", "w!");
    save_fl["eigenvalues"] = eigs;

    return 0;
}
catch (Error e)
{
    error_trace(e);
}

using XDiag
using LinearAlgebra
using Printf
using HDF5

Sx = [0.0 0.5; 0.5 0.0]
Sy = [0.0 -0.5*1im; 0.5*1im 0.0]


function main()
    L = 6 # total number of sites

    gamma = 2.0 # dissipation strength

    Nsites_liouville_space = 2 * L ## the effective Hilbert Space is doubled!

    nup = Nsites_liouville_space ÷ 2

    block = Spinhalf(Nsites_liouville_space, nup)

    #build Hamiltonian
    H = OpSum()
    for i in 1:L
        H += -1.0*1im*Op("Matrix", [i, mod1(i+1, L)],kron(Sx,Sx)) ## normal spin chain
        H += -1.0*1im*Op("Matrix", [i, mod1(i+1, L)],kron(Sy,Sy))

        H += 1.0im*Op("Matrix", [L+i, L+mod1(i+1, L)],kron(Sx,Sx)) ## normal spin chain
        H += 1.0im*Op("Matrix", [L+i, L+mod1(i+1, L)],kron(Sy,Sy))

        H += gamma*Op("SzSz", [i, L+i]) ## recycling term
    end


    Hmat = matrix(H, block)

    evals = eigvals(Hmat)    

    evals .-= gamma*L/4

    filename = @sprintf("./heisenberg_bulkd_dephasing.Nsites.%d.gamma.%.3f.nup_sector.%d.outfile.h5", L,gamma,nup)

    h5open(filename, "w") do file
        write(file, "real_part_eigenvalues", real.(evals))
        write(file, "imag_part_eigenvalues", imag.(evals))
    end
end

main()

references

[1] Fabrizio Minganti, Alberto Biella, Nicola Bartolo, and Cristiano Ciuti, Physical Review A 98, 042118 (2018)

[2] Rafael D. Soares, Matteo Brunelli and Marco Schirò, arXiv:2505.15711 (2025)

[3] Marko Žnidarič, Physical Review E 92, 042143 (2015)