# Extracting total parity from MPS

+1 vote
edited

Hi Miles,

I'm trying to calculate the total parity of a MPS. For example, consider a state with even parity |even> = a|00>+b|11>. The naive calculation of (\sum < n_j > )% 2 gives (2*|b|^2) % 2, not 0 as expected. What is the correct way of doing this? Thanks!

Best,
Chengshu

selected by

Hi Chengshu,
Could you say more about how the parity you want to measure is defined? If I could see the form of the operator you want to measure, then I can tell you more about a good strategy for doing the measurement.

If it is defined as (\sum ) %2, then if it's not giving the answer you expected it could be one of any number of things. For example, maybe the code to do the measurement is incorrect (e.g. you didn't move the orthogonality center of the MPS by calling .position). Or maybe the QNs you are conserving don't include parity?

Best regards,
Miles

commented by (680 points)
edited by
Hi Miles,
Basically what I want to calculate is exp(i \pi \sum n_j) = \prod(1 - 2 n_j). So it is by definition a nonlocal operator.

In a spinless 2-site system, for example, |00>, |11> and all linear combinations of these two have parity 1, while |01>, |10> and their linear combinations have parity -1.
commented by (46.8k points)
Hi Chengshu, ok now I see what you are trying to measure, thanks. The reason measuring the density did not work is that the operator expinentiaion does not commute with taking an expectation value. I.e. the parity operator is sensitive to higher order correlations in the wavefunction not captured by the sum of the local densities.

But to measure an operator like this is relatively straightforward with an MPS. First write the operator as a product of local site operators, using the fact that the n_i operators commute. (I see you've already done this above.) Then sketch out the tensor diagram for the expectation value of the operator. It will look like the diagram for the overlap of the wavefunction with itself, except at every site one of the (1-2n_i) operators is inserted. Finally, compute this diagram by contracting tensor together from left to right, say.

auto P = psi.A(1) * op1 * dag(prime(psi.A(1)));
P *= psi.A(2);
P *= op2;
P *= dag(prime(psi.A(2)));

Etc.

Here by op1, op2, etc I mean the (1-2*n_i) operators which you can make either on the fly or ahead of time. You could make them from operators provided by the site set or by explicitly setting tensor elements.

At the end, P will be a scalar tensor if all goes well and you can call P.real() to extract its value.

Let me know if you have some more questions about it.
commented by (680 points)
Hi Miles, the code works! Thank you very much!
commented by (46.8k points)
Glad to know it works! If you have time and energy to clean up your code and send it to me, I can turn your measurement code into a "Code Formula" on the ITensor website. Other people might like to see it so they can do the same type of measurement in the future.