How to calculate flux of a reduced density matrix?

Hi Miles,

Sorry for the late update. It seems i have some troubles for obtaining rdm and its factorization.

  1. For obtaining rdm from a MPS |\psi\rangle, a naive way is to call outer(psi',psi) (for example, How to get density matrix operator from MPS wave function? - ITensor Support Q&A) and partially traced out half of the system. (for example, the discussion here trace and partial trace of MPO - #3 by jisutich) However, this method seems to be very slow when the bond dimension is large.

    Therefore i was trying to make a conjugate copy \langle\psi| and partially trace out the indice from 1 to k (or from k+1 to length(psi)) with the help of delta(s[i],s[i]'). I tried

using ITensors
n = 10
s = siteinds("S=1/2",n;conserve_qns=true)
state = ["Dn" for i = 1:n]
psi = randomMPS(s,state)

k = 5
for i = 1 : k
   psi[i] = prime(dag(psi[i]),siteinds(psi))*delta(s[i],s[i]')*psi[i]
end

But it seems failed. ( sad face…)
So could you pleas give me more hints on finding rdm from a given |\psi\rangle ?

  1. Another naive question is for the factorizing the rdm. We know one can obtain the entanglement spectrum by simply perform svd decomposition. (von Neumann entanglement entropy of 1D hubbard model (Julia) - ITensor Support Q&A) But for a MPO \rho here, should I first convert it to a matrix and perform eig-like function and then convert it back to a MPS?
    (for example, see there Converting ITensors to matrices in Julia - ITensor Support Q&A)

  2. And for a MPS, one can check its quantum number by calling flux(), but how to extract it as a float number?

I guess the last two questions are due to my ignorance about julia.

Thanks,
Zach