How to calculate flux of a reduced density matrix?

Hi Matt and Mile,

I am trying to reproduce the Laughlin charge pumping on a cynlinder with DMRG. Following Frank Pollmann’s paper Phys. Rev. Lett.115.116803

The charge can be obtained by

\langle Q(\Phi_y)\rangle =\sum_i \lambda_i(\Phi_y)Q_i(\Phi_y)

with \lambda_i(\Phi_y) are the eigenvalues of the reduced density matrix and Q_i(\Phi_y) is the corresonding U(1) quanutm number.

We know there is a method called flux() or totalqn() to calculate the flux of a MPS. So my naive question is that is there a method to calculate the flux of a reduced density matrix or get individual flux on every site of a MPS ?

Any comments are welcome :slight_smile:

Best,
Zach

Good question: so for any ITensor you can also call flux. So if you obtain a reduced density matrix (rdm) as an ITensor rho you can call flux(rho).

Also if I’m reasoning about it correctly, it should be the case that all rdm’s obtained from some pure state have a flux of zero. My reasoning is that for some full density matrix \rho = |\Psi\rangle \langle \Psi | if |\Psi\rangle has some flux q then \langle \Psi | has a flux of -q and the whole density matrix \rho has zero flux. Then tracing out any subset of sites to get an rdm will also give a tensor with zero flux.

So perhaps in the formula you cited above, the quantum number Q_i(\Phi_y) means the quantum number that flows “through” an rdm \rho in the sense that if we factorize it as U \Lambda U^\dagger then one can associate quantum numbers with the index connecting U to \Lambda and those might be the Q_i.

1 Like

(Just edited my answer above, adding the last paragraph, so just wanted to note that.)

Hi Miles,
Thanks for the comments and the last paragraph really helps a lot! I think your suggestions probably will lead to the observables in formula. I will have a try then and update the results here :stuck_out_tongue:

Best,

Zach

1 Like

Yes please let us know how it goes. I’m not totally sure I’m right about that formula so please check independently but it’s my best guess of what they are saying in the paper (without having read the whole thing).

1 Like

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

Hi Zach,
Thanks for the follow-up question here. So for the case of a pure state which is an MPS, there is a shortcut to getting the kind of rdm you want (where sites 1:k are traced over). The shortcut is to “gauge” the MPS so that the tensors on sites 1,2,…,k are all “left-orthogonal”. This can be done in ITensor by calling orthogonalize!(psi,k+1) on an MPS psi. Afterward, the “orthogonality center” will be at site k+1, meaning that the tensors 1,2,…,k will all be left-orthogonal (and tensors k+2,k+3,...,N will also be right-orthogonal).

After that step, you can show that tracing sites 1:k just amounts to leaving out those tensors (they cancel) and just contracting MPS tensor number k+1 with its primed conjugate over the virtual or link index on bond k. If we call the virtual index for bond k+1 b_{k+1} and the site index for site k+1 s_{k+1} then the resulting rdm will have indices s_{k+1},b_{k+1}, s'_{k+1},b'_{k+1}. Then you can input this rdm to our eigen function in ITensor and factorize it to get its eigenvalues.

I’d have to play around with the code to remember how to extract the QNs. Once you have the QN objects themselves, you can call the val function on them to get their numerical values. Say the name of the QN you want is "Sz" so like inside of a QN object like q-QN("Sz",-1,"Nf",1). Then you do val(q,"Sz") and it will return the number -1.

2 Likes