# How to calculate flux of a reduced density matrix?

**URL:** <https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243>\
**Category:** ITensor Julia Questions\
**Tags:** julia, dmrg, mps\
**Created:** [July 10, 2022, 4:21pm UTC](https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243 "2022-07-10T16:21:46Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![ZacheryZhou](https://avatars.discourse-cdn.com/v4/letter/z/cdc98d/32.png) [@ZacheryZhou](https://itensor.discourse.group/u/ZacheryZhou)\
**Post date:** [July 10, 2022, 4:21pm UTC](https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243/1 "2022-07-10T16:21:46Z")

</div>

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

> **[Bosonic Integer Quantum Hall Effect in an Interacting Lattice Model](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.115.116803)**
>
> We study a bosonic model with correlated hopping on a honeycomb lattice, and show that its ground state is a bosonic integer quantum Hall (BIQH) phase, a prominent example of a symmetry-protected topological (SPT) phase. By using the infinite density...

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 🙂

Best,  
Zach

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [July 10, 2022, 7:36pm UTC](https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243/2 "2022-07-10T19:36:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [July 10, 2022, 9:53pm UTC](https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243/3 "2022-07-10T21:53:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ZacheryZhou](https://avatars.discourse-cdn.com/v4/letter/z/cdc98d/32.png) [@ZacheryZhou](https://itensor.discourse.group/u/ZacheryZhou)\
**Post date:** [July 11, 2022, 2:41pm UTC](https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243/4 "2022-07-11T14:41:18Z")

</div>

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 😛

Best,

Zach

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [July 11, 2022, 3:39pm UTC](https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243/5 "2022-07-11T15:39:45Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![ZacheryZhou](https://avatars.discourse-cdn.com/v4/letter/z/cdc98d/32.png) [@ZacheryZhou](https://itensor.discourse.group/u/ZacheryZhou)\
**Post date:** [July 13, 2022, 2:39pm UTC](https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243/6 "2022-07-13T14:39:21Z")

</div>

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](http://itensor.org/support/3217/how-to-get-density-matrix-operator-from-mps-wave-function?show=3217#q3217)) and partially traced out half of the system. (for example, the discussion here [trace and partial trace of MPO - #3 by jisutich](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/3)) However, this method seems to be very slow when the bond dimension is large.

```auto
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](http://itensor.org/support/3895/von-neumann-entanglement-entropy-of-1d-hubbard-model-julia?show=3895#q3895)) 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](http://itensor.org/support/2885/converting-itensors-to-matrices-in-julia?show=2885#q2885))

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

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [July 14, 2022, 1:31pm UTC](https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243/7 "2022-07-14T13:31:42Z")

</div>

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](https://itensor.github.io/ITensors.jl/stable/QN.html#ITensors.val-Tuple%7BQN,%20Any%7D) 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`.