Local parity operator expansion for electron site parity operator

Continuing the discussion from Issues with Nf parity conservation:

Dear Miles,

I have a comment and a question regarding @vernek 's question on parity conservation.
Am I correct that in Edson’s case the expectation value of the total parity operator e.g.
\displaystyle (-1)^{N}=\prod_{m,\sigma}(1-2n_{m,\sigma}) should capture the appropriate parity eigenvalue ?
If so what is the simplest way of implementing it?
I spent a couple of hours with a sanity check regarding parties in a related problem e.g. I also have superconductivity so parity is the only conserved quantity.
In the long run, I want to calculate parity content associated with a certain part of my system, so a nice operator built up of local parts would be nice.

I have came up, based on the documentation, the following code snippet:

s = siteinds(PSI)
PSI0=PSI
for i in 1:Length
    # this is (1-2n_up)(1-2n_dn)=(1-2n_tot-4Nupnd)
    O=op("Id", s[i])-2*op("Ntot",s[i])+4*op("Nupdn",s[i])
    psi=O*PSI0[i]
    noprime!(psi)
    PSI0[i]=psi
end

inner(PSI,PSI0)

it seems however that the last inner product is always 1. independent of the actual parity sector my original wavefunction PSI0 is located.

Thanks for the question. Some simple ways of implementing a product operator of the form you wrote, so some operator of the form W = \prod_j O_j, would be:

  1. making a copy of your original MPS, then looping over each site of the copied MPS and applying each O_j onto the j’th MPS tensor, then using inner to overlap this modified MPS with your original MPS

  2. putting each O_j into an MPO, i.e. making each MPO tensor equal O_j, then using the inner(psi,W,psi) to compute its expected value

Also, for what it’s worth, you can be sure that if you are using conserved quantum numbers which either explicitly conserve parity, or something that includes parity like total particle number, then if you call flux(psi) on your MPS psi, then it will accurately return the total quantum number flux (total QN) of your MPS. If that flux is in the parity sector you want, that the overall MPS is also guaranteed to be in that parity sector.

Thanks, Miles for the prompt answer!
I am aware that I can check the total parity with flux(psi), but my goal is to get the parity of two disjoint subsystems.
This might not be well defined in general, but I have good reason to believe that for my special case, it is going to deliver a proper parity eigenvalue, e.g. \pm 1.
Regarding your suggestion in your first answer, I believe that my code included in the original post does just that, yet it does not give me the appropriate answer. I prepared a more complete example below using a random state with a fixed parity

using ITensors


N=10 # make this even just for the sake of concreteness
@time sites = siteinds("Electron", N;
                         conserve_nfparity = true,
                         conserve_sz = false,
                         conserve_nf = false);


state = ["Up" for i in 1:N]
state[1]="Emp" # odd state if this line is commented even state can be checked
psi0 = randomMPS(sites, state, 5)
@show flux(psi0)   # this shows that the flux is  odd                    

PSI=psi0

s = siteinds(PSI)

for i in 1:N
    # this is (1-2n_up)(1-2n_dn)=(1-2n_tot-4Nupnd)
    O=op("Id", s[i])-2*op("Ntot",s[i])+4*op("Nupdn",s[i])
    psi=O*PSI[i]
    noprime!(psi)
    PSI[i]=psi
end

inner(PSI,psi0) # this should be -1 for odd and +1 for even

I found my bug. Instead of PSI=psi0 I needed PSI=copy(psi0). I guess this solves my Issue! Thanks for the suggestions Miles!

Glad you found the bug and the code is working now.