How to get density matrix operator from MPS wave function?

Question

Hello all,

I'd like to ask a follow-up on the question posted here

http://itensor.org/support/2641/construction-of-matrix-product-density-operator?show=2641#q2641

since I do not fully understand the reply to that post.

Suppose I have a set of states obtained with a procedure like the one exemplified here:

https://itensor.org/docs.cgi?vers=cppv3&page=formulas/excited_dmrg

so I would have a number of excited states. Now, how would I go about writing a thermal density matrix where these states are added together with a thermal weight?

It seems to me we have a set of wave functions and that the density matrices for each of them should be written as a tensor product and then summed over, but it's not clear to me how that would be done with the MPS set resulting from the excited state search. Is there a straightforward way to create an MPO from an MPS? I checked the pages on manipulation of MPOs but could not find an answer.

Bonus question: assuming I have such an MPO describing a mixed state, how to get the expected value of an operator?

Best wishes,

Rafael

Comments

Hi Rafael,
Good to see Matt's excellent and detailed answer below. In case you're interested to discuss more,  I was wondering why you wanted to sum up a set of eigenstates into a single MPO? Is it because you want to perform some operation on that MPO or use it as input for another code that needs an MPO? Or do you want to compute something like the thermal average of a local observable? Just curious –

Miles

---

Hi Miles,

>> Or do you want to compute something like the thermal average of a local observable?

This is what I want to do (although for a global observable). I know there are other ways to do this (such as the ancilla method) but I want to calculate the density matrix for each state, as Matt wrote below. I'm sure there's a simple way to calculate the outer products in the c++ version, but I've seen this question come up a few times around here and the answer is not clear to me. Any help is appreciated.

Best wishes,

Rafael

---

Hi Rafael, I see thanks. So if you just want a thermal expectation value, then you don't actually need to make the density matrix first as a separate step, which could be very expensive. Instead you can take your pure (eigen) states and get the expected value of your operator within each one of them, then weight these values by exp(-beta*H). Please see these notes to show the derivation of this:

https://itensor.org/miles/EigenFiniteT.pdf

In the last expression in the notes, you can see that the expected value of the operator A (regardless of what kind of operator it is) can be written as a sum of expected values of A with respect to each one of the eigenstates. If you 'chop' the sum over eigenstates to only sum over a finite number of them, you would want to correspondingly compute Z with the same truncated sum so that taking A to be the identity operator gives 1.0 as it should.

---

Hi Miles,

Thanks again for the reply. That is for sure simpler, and I don't see any reason why it would not work too. Let me give it a try and see how it goes.

Really appreciate the help!

Rafael

Answer 1

Say that you have MPS psi1, psi2, psi3 as excited states. You can use the outer function (https://itensor.github.io/ITensors.jl/stable/MPSandMPO.html#ITensors.NDTensors.outer-Tuple{MPS,%20MPS} ):

rho1 = outer(psi1, psi1)
rho2 = outer(psi2, psi2)
rho3 = outer(psi3, psi3)

then you can scale them by constants c1, c2, c3 and add them up with:

rho = +(c1 * rho1, c2 * rho2, c2 * rho3; cutoff=1e-12)

Finally, I will leave it as an exercise to compute properties of the MPO by generalizing the case for MPS (http://itensor.org/docs.cgi?vers=julia&page=formulas/measure_mps ). I would recommend drawing the diagrams to start out. A hint is that in order to compute a local expectation value of an MPO, you will need to do a partial trace of the degrees of freedom of the MPO that are in the complement of the set of sites where the local operators you are interested in are located. One way to perform the partial trace of an MPO is by contracting the pairs of site indices on each site by a delta tensor.

Comments

Hello Matt,

Thanks for the fast reply!

I am assuming that the outer() function is valid for the C++ version too, which I am using. I will try your recommended method for calculating observables with an MPO. Since I am interested in calculating the expected value of a global operator, I am thinking it might also work to multiply the MPOs with nmultMPO() and then simply applying the trace() function.

Best,

Rafael

---

I'm actually not sure if there is an `outer` function in the C++ version, but it is pretty simple to write it (for example, prime one of the MPS, contract the tensors pairwise at each site, and then compress the result).

For a global operator, using AutoMPO should be a good way to go. A more efficient alternative to `nmultMPO` would be to use the version of `trace` that accepts two MPOs `A` and `B` and directly computes `tr(AB)`: https://github.com/ITensor/ITensor/blob/684bb4207da536dce59f1f2589b60ec23d1d61cb/itensor/mps/mpo.cc#L566-L605

---

Hello Matt,

Thanks again. About this:

>> prime one of the MPS, contract the tensors pairwise at each site, and then compress the result

could you provide an example in code form or the reference page where these operations are defined?

Best,

Rafael

---

Please take a look through the example codes here: http://itensor.org/docs.cgi?page=formulas&vers=cppv3 to get an idea for some of the operations and let us know where you get stuck.