Variational MPO as ground state

[4psireal2]

I'm doing open dynamics simulations with the density matrix represented as vectorized MPO as described in this publication. The main idea is that the density matrix is an mps with 4 legs (2 physical legs pointing down and 2 horizontal virtual legs) and the Hamiltonian (or Linbladian) is an mpo with 6 legs (2 incoming legs, 2 outgoing legs, and 2 virtual legs). I have implemented this approach with TensorKit for my system together with DMRG1 and DMRG2. It works fine with DMRG1 but the simulation doesn't converge for DMRG2. I'm hoping to test this out with other DMRG2 implementation.

I'd greatly appreciate if you could please outline how this setup might be achieved with MPSKit.

[maartenvd]

There is no "real" support for linbladians, in the sense that all implemented operators (mpo/hamiltonian) act on MPS tensors with one physical leg. You can either fuse the two physical legs together, or use a nontrivial unit cell to split your MPS tensors (2 physical legs) into two 1-physical leg tensors. Do you have the dense MPO representation of the Lindbladian?

That said, there are serious problems with combining that publication with any infinite algorithms like DMRG1, DMRG2, VUMPS,... The problem is that your MPO tensor is not well behaved (it has an eigenvalue scaling as N^2 in system size). MPSKit's DMRG1/DMRG2 will fail to even initialize correctly in such a case, as we use the left/right dominant eigenvectors as initial start. You can probably work around that, but I kind of suspect that the numerical issues you're seeing with your own dmrg2 implementation are inherent to the kind of MPO tensor you constructed.

I don't remember the publication, but followup work proposed simply truncating the infinite double sum that occurs in your mpo. That's also a bit iffy, but infinite size algorithms should no longer struggle as much.

[4psireal2]

Yes I've tried fusing the two physical legs

γ = 1.0
Ω = 2.5
Id_mat = TensorMap(ComplexF64[1 0;0 1], ℂ^2 ← ℂ^2)
fuseIsometry = isomorphism(ℂ^4, ℂ^2 ⊗ ℂ^2) 

sigmax = 0.5 * TensorMap(ComplexF64[0 1; 1 0], ℂ^2 ← ℂ^2)
sigmax_r = fuseIsometry * (Id_mat ⊗ sigmax) * fuseIsometry'
sigmax_l = fuseIsometry * (sigmax ⊗ Id_mat) * fuseIsometry'

number_op = TensorMap(ComplexF64[0 0; 0 1], ℂ^2 ← ℂ^2)
number_op_r = fuseIsometry * (Id_mat ⊗ number_op) * fuseIsometry'
number_op_l = fuseIsometry * (number_op ⊗ Id_mat) * fuseIsometry'

annihilation_op = TensorMap(ComplexF64[0 1; 0 0], ℂ^2 ← ℂ^2)
creation_op = TensorMap(ComplexF64[0 0; 1 0], ℂ^2 ← ℂ^2)

on_site = γ * fuseIsometry * (annihilation_op ⊗ annihilation_op) * fuseIsometry' - (1/2) * number_op_r - (1/2) * number_op_l
on_site_dag = γ * fuseIsometry * (creation_op ⊗ creation_op) * fuseIsometry' - (1/2) * number_op_r - (1/2) * number_op_l

Lindbladian = Array{Any, 3}(missing, L, 6, 6)
Lindbladian[1, 1, :] = [1, sigmax_r, sigmax_l, number_op_r, number_op_l, on_site] 
for i = 2:L-1
    Lindbladian[i, 1, :] = [1, sigmax_r, sigmax_l, number_op_r, number_op_l, on_site]
    Lindbladian[i, 2, 6] = -1im * Ω * number_op_r
    Lindbladian[i, 3, 6] = 1im * Ω * number_op_l
    Lindbladian[i, 4, 6] = -1im * Ω * sigmax_r
    Lindbladian[i, 5, 6] = 1im * Ω * sigmax_l
    Lindbladian[i, 6, 6] = 1
end
Lindbladian[L, :, 6] = [on_site, -1im * Ω * number_op_r, 1im * Ω * number_op_l, -1im * Ω * sigmax_r, 1im * Ω * sigmax_l, 1]
Lindbladian = MPOHamiltonian(Lindbladian)

Lindbladian_dag = Array{Any, 3}(missing, L, 6, 6)
Lindbladian_dag[1, 1, :] = [1, sigmax_r, sigmax_l, number_op_r, number_op_l, on_site_dag] 
for i = 2:L-1
    Lindbladian_dag[i, 1, :] = [1, sigmax_r, sigmax_l, number_op_r, number_op_l, on_site_dag]
    Lindbladian_dag[i, 2, 6] = 1im * Ω * number_op_r
    Lindbladian_dag[i, 3, 6] = -1im * Ω * number_op_l
    Lindbladian_dag[i, 4, 6] = 1im * Ω * sigmax_r
    Lindbladian_dag[i, 5, 6] = -1im * Ω * sigmax_l
    Lindbladian_dag[i, 6, 6] = 1
end
Lindbladian_dag[L, :, 6] = [on_site_dag, 1im * Ω * number_op_r, -1im * Ω * number_op_l, 1im * Ω * sigmax_r, -1im * Ω * sigmax_l, 1]
Lindbladian_dag = MPOHamiltonian(Lindbladian_dag)

Lindbladian_hermit = Lindbladian_dag*Lindbladian