Calculating observables from density matrices

Question

Hey everyone, I am rather new to MPS so please forgive the naive question. I want to replicate the results from the paper, https://arxiv.org/pdf/1702.04210.pdf, to learn MPS methods. The primary quantity of interests is the spin density,
$$ s(x, t) = \text{tr} (\rho (t) S^{z}_{x})$$, where x is the spin-site. In this work, they prepare the system in a "soft" domain wall density matrix.

My approach: The way I went about this is by thinking of the density matrix as an MPO and then evolve this in time using trotterisation methods. The first problem, I have is computing the $$\text{tr} (\rho (t) S^{z}_{x})$$. You can't just use the "*" method to contract the MPO's and then trace.

My first question is how to perform this operation correctly?

My solution attempt: I just went ahead and calculated real(inner(rho(t), S^z)) for each site which matches small scale exact diagonalization methods, but I am not sure if this is still correct. Here I just used a density matrix made from a Neel state

My second question is how to go about making this initial density matrix as an MPO?

My solution attempt: I first wrote a function in Julia that can create the matrix for this soft domain wall state. I then define the necessary site indices and convert to an MPO but it doesn't have the correct structure and also this won't work for large N because it explicitly makes the matrix.

Any feedback or ides is much appreciated

Comments

Hi, so I think your question has two parts: 1. computing local observables from a density matrix MPO and 2. preparing an MPO to be some specific initial density matrix, correct?

Regarding 1, your method of using inner(rho(t),S^z) is a correct procedure. The inner function for MPO's computes the quantity inner(A,B) = Tr[A^\dagger B] which is known as the Frobenius inner product. (You can also view the MPO's as having their physical indices all folded into the same direction to "vectorize" them and then the inner product is computed by treating them as MPS with two indices per site.)

If you are just measuring S^z_j on individual sites j, there is a more efficient way however than representing S^z_j as an MPO. Our documentation could benefit from a code example about this, but what you would do is trace all of the other sites of the MPO and contract the MPO tensor at site j with your operator S^z_j represented as a single ITensor with two indices.

Finally, regarding how to prepare an initial density matrix, the answer depends on which one you are trying to prepare. Could you provide more detail about that?

---

Note that I've started an issue related to this question here: https://github.com/ITensor/ITensors.jl/issues/730

I've made a simple example code in that issue for computing the expectation value of an MPO.

---

Hey Matt,

Thank you for posting that little code. I think in your section where you say 'used like:' you have an error. I think it should be, expect(M, "Sx", 3) instead of "op(M, "Sx", 3)". Perhaps I am not understanding correctly.

---

Yes that's right, thanks for pointing that out, I've fixed it.

---

Thanks, you guys for the insightful comments let me address Miles question. The initial density matrix I want to prepare is the following,
$$\rho \sim (1+\mu S_{z})^{\otimes N/2} \otimes (1-\mu S_{z})^{\otimes N/2}$$, where $$\mu \in [-1,1]$$. You can see that for $$\mu \sim \pi/1800$$ this is effectively the infinite temperature state on each site.

Also, how exactly can I control the maximum bond dimension with this kind of evolution? The reason I ask is that if I do apply(gates, rho; cutoff=cutoff, maxdim=maxdim, apply_dag=true) for evolution I get the same curve if I change maxdim. I would expect that for larger max bond dim I get a better approximation to the dynamics.

Chris

---

`cutoff` and `maxdim` are the correct ways to control the truncation/maximum bond dimension. Have you tried adjusting both the `cutoff` and `maxdim`, or just `maxdim`? We would need more details about your calculation to understand what trouble you might be having with the truncation/convergence with bond dimension (for example, I'm not sure what curve you are trying to compute).

Note that in general with real time evolution of tensor networks, the entanglement grows linearly with time, so generally you need exponentially large bond dimensions to reach longer and longer times, though the system you are looking at may be different.

---

Sorry I am rather new to tensor networks so I may not understand everything correctly yet. But what I want to calculate is the following:
$$C(x,t)=\langle \rho (t) S^{z}_{x} \rangle$$, where x indicated the site of the chain and $$\rho(t)$$ is the evolved MPO discussed above. Specifically, I want to fix the site position and plot this curve for different bond dimensions to see the convergence.

---

I see. So when you plot @@C(x, t)@@ for different bond dimensions, you are seeing that the curve does not change with bond dimension? Is that true for every time @@t@@? At later times, the state likely becomes more entangled, so you may not be running to long enough time scales to see the entanglement saturate the bond dimension you set. Also, it is helpful to look at the entanglement of @@\rho(t)@@ as well as the smallest singular value of @@\rho(t)@@ as a function of time to get an idea for how complex your state is to represent.

---

Hey Matt,

I based my TEBD code off the gate evolution example for an MPO, but when I try to use conserve_qns=true, I get a weird svd error which says, "In 'svd' the left or right indices are empty." Should I open a new issue for this?

thanks!

---

Hi, I’d say this is a good candidate for opening an issue. We may find that it’s not a bug, but rather something about the setup of your code or inputs, etc. but it would be good in this case for us to have a look. Thanks.

Please do see if you can find the simplest case or version of your code that reproduces this issue.

Miles

Answer 1

See the discussion above.

Answer 2

For the part of your question (in comments above) about preparing the initial state MPO, since it is a tensor product of local operators, you can make each operator separately as an ITensor and then "load" them into an MPO by assigning the tensors of the MPO one by one (in a for loop of course) to be the operators you want.

For example, just showing the first site:

N = 100
s = siteinds("S=1/2",N)

rho = MPO(N)

rho[1] = op("Id",s[1]) + mu*op("Sz",s[1])

and similar for sites 2,...N/2 in a for loop (use div(N,2) to get N/2 in Julia). Then another loop to go from N/2+1 to N with the other sign for the Sz term.

Some of our algorithms may complain if there are not virtual or link indices between the tensors of the MPO. If that happens, let us know and we can recommend a course of action. Thanks!

Miles

Comments

Hm I see this will come from more experience. I want to TEBD evolution, so I am just creating the two sites gates and applying across the chain and then doing the reverse. I then apply the dagger to get Heisenberg evolution of an MPO. But after each time-step is it necessary to orthogonalize the MPO?

---

Hi, so in TEBD is it very important to keep the MPS or MPO in an orthogonal form. The part about applying the gates will be fine and always numerically exact, but the part where you factorize and truncate the state afterward can incur uncontrolled errors if the other tensors to the left and right are not left- and right-orthogonal.

Note that when doing TEBD on an MPO, one is really treating it as an MPS with two physical indices on every MPS tensor, in terms of how the canonical forms are defined. There aren’t really so many useful canonical forms specific to MPOs besides this MPS-derived one.

---

I see, thank you for taking all this time to answer probably rather simple questions. Currently, after I evolve the MPO I take the expectation value with some spin operator on each site. What I did was to write a function that takes the (MPO, "Sz", sitesindex, site) where the MPO is the evolved state. I always first orthogonalize the MPO before computing the expectation value. It also seems the using number conserving for my initial density matrix does not work as well.

---

I see - sounds good. I had thought you were asking about keeping an orthogonal form during the time evolution itself, which is very important to do, as I think you saw already from my answer above.

---

Correct, for the evolution I just followed the example given on one of the pages for the quantum simulator on github.