# Spatial parity inversion of MPS using SWAP gates

+1 vote

I have a model of two sets of energy levels coupled to an impurity (like a two channel Kondo model). The MPO is implemented such that the MPS has a L(eft channel) - imp(urity) - R(ight channel) geometry. I am using the c++ version of ITensor.

When the physical parameters are symmetric, this system has spatial inversion symmetry. To measure it, I want to implement a parity operator that swaps the L and R channels in the MPS. For 7 sites the action of this operator would transform the mps (1, 2, 3) - 4 - (5, 6, 7) into (5, 6, 7) - 4 - (1, 2, 3). Measuring the overlap of the new state with the original should give me parity.

I implemented the permutation of channels by using SWAP gates to permute neighbouring sites so that the channels get swapped. I tried using the BongGate class and also the svd trick mentioned in Miles' comment from Sep 13, 2019 here. I check that the permutation is correct by comparing occupation numbers (or even local density matrices) of sites in the old and new state. Then, I compute parity as inner(psi, newpsi).

There are two problems:

a) The approach does not seem to be correct. For symmetrical parameters, the parity should be 1 or -1, but the result is typically at least a few percent off (depends on physical parameters and symmetry sector chosen - I have observed that taking half-filling gives results closer to 1 (but still on the order of 0.98), but am not sure whether this is relevant). I also tried permuting the siteIndeces of the MPS using the delta tensor as indicated here, but it does not seem to do anything - the overlap always gives 1, even with asymmetric physical parameters.

b) The approach will not work for large system sizes, as permuting the levels with SWAP gates blows up the MPS in memory.

My question is two-fold.
What is wrong with my approach?
And what would be a better way to obtain parity?

commented by (70.1k points)
Hi, so it’s a good question but I don’t think the answer is too obvious right away. But I will think about it. Definitely starting with swaps on small systems is a good way to go, and maybe there is a more clever trick or technique that can scale to larger systems.

One question back to you: I’m surprised that the symmetry would actually change which sites are next to or connected to the impurity. In the original ordering sites 3 and 5 are the neighbors of the impurity but in the second ordering sites 7 and 1 become the neighbors. So it’s not just a reflection of the system but like a reconnection of the interactions altogether. Or is it? Just wanted to ask since I don’t know a lot about this kind of symmetry.

Thanks, Miles
commented by (130 points)
Yes, there is a disconnect between the symmetry of the hamiltonian, its states being gerade or ungerade orbitals, and the symmetry of the order in which tensors appear in the MPS.
In the above example, levels 1, 2 and 3 are coupled to level 4 by hopping v_L, while 5, 6 and 7 are coupled to 4 by hopping v_R. But they also have single particle energies, with pairs (1, 5), (2, 6) and (3, 7) having the same.  A spatial inversion is a transformation that maps v_L <—> v_R and the left channel onto the right and the opposite. This is expressed (I think) as the permutation of the MPS as stated above.

Do you think it would be better to have the MPS with tensors ordered in a symmetric way? This would be a relatively cheap extension of our code.

For more detail, the model is explicitly written and described on the second page of https://arxiv.org/pdf/2110.13881.pdf, with the MPO given in appendix B, page 17.
commented by (70.1k points)
Yes, if you could order the tensors in a symmetric way then I think I know a nice trick you could use. Here by symmetric I mean such that 7-6-5-4-3-2-1 would be swapped ordering and 1-2-3-4-5-6-7 would be the original ordering.

auto psi_swap = MPS(7);
and then write a loop that fills it up with the tensors of the original MPS but in reverse order, so like:
psi_swap.ref(1) = psi(7);
psi_swap.ref(2) = psi(6);
and so on either explicitly or in a loop.

Finally, you could use 'delta' tensors (special kind of ITensor for replacing site indices) to replace the site indices of psi_swap with the original ones in the original order. So like:

psi_swap.ref(1) = psi_swap(1) * delta(dag(sites(7)),sites(1));

and then such a code would very cheaply make an swapped state MPS for you.

Finally you could compute properties of the swapped state such as its overlap with the original state because they have compatible site indices (but underneath, the swapped state will have swapped properties).

Hope that idea works!
commented by (130 points)
Sorry for a late response.
I do not think this approach works either. The results are equivalent to what I got before - even tough all local observables are symmetric, the overlap with the swapped state does not give 1.

I also tried swapping all link indices in the same way as the site ones, so the two MPSs have exactly the same tensor structure, but unsurprisingly that does not affect the result.

By the way, applying delta( sites(i), sites(i) ) raises an error: mismatched index arrows. I would expect this to give identity (ie. do nothing to the indices), but it seems like iTensor identifies the two "out" indices as the ones to contract instead of one "out" and one "in", and thus fails.
Is that the expected behaviour?
commented by (70.1k points)
Hm, it sounds like maybe the reason you aren't getting 1 for your overlap might have a different origin than the procedure for swapping the sites. Could you please check other parts of the code to see if there could be bugs there? I'm sure you have done that quite a bit already of course.

Back to the swapping, though, one thing you could do to rule out it as the culprit of why you aren't getting 1 is to create an MPS 'by hand' which you know for sure has the right symmetry, then input it to the part of your code that swaps the sites. Then if you get 1 for this artificial input you will be more confident about whether that part of the code is correct or not.

Finally, about the arrows we don't technically allow two identical indices to be on an ITensor, even if they have the same arrows. Although it would be convenient, it leads to ambiguous cases with our tensor contraction system's logic.

So the reason for the error is that you will need to do delta(dag(sites(i)), sites(j)) to replace the sites(i) Index with the sites(j) Index. (I'm not sure why you would want to replace the sites(i) Index with itself.)

Hope that helps!

Miles
commented by (130 points)
Thank you for the comment about trying things on known states! Somehow I have not thought of that. I tried on small clusters and by inspection of amplitudes I found the problem.

The inversion of my state is not just swapping of tensors, as one has to also take care of the fermionic minus sign.
Representing psi as a string of creation operators acting on the vacuum state, the inversion corresponds to changing their site index i --> N-i, AND commuting the string back into canonical order. The fermionic minus appears when commuting i past j, but only when both i and j are singly occupied.

Now if we define O_i as a diagonal operator which gives 1 for single occupation and 0 elsewhere, the minus sign structure of a permutation of i by j is given by a two level operator F(i, j) = 1 - 2 O_i O_j.

Parity inverting my MPS thus means first applying a string of F(1, 2) F(1, 3) ... F(1, N) (to commute level 1 to the end), then F(2, 3) ... F(2, N) and so on untill F(N-1, N). This will effectively multiply the correct amplitudes by -1.
The next step is to actually swap the tensors and the site indices, as discussed previously.

This approach works, and is not terribly slow or costly of memory (I think).
But as I implemented it currently, each F(i, j) consists of applying two single site operators (which is fine) and then summing up two MPS to get 1 - 2*O_i O_j. I don't think it is possible to write the expression as a product of two single-site operators.
So my question is whether you know of a way to improve this procedure? Do you think using autoMPO to obtain an MPO representation of the string of F(i, j) would be quicker? I haven't tried writing it as an MPO by hand yet, but that might be possible too.