Fix expectation values for PEPS/PEPO with dual physical spaces

[Yue-Zhengyuan]

This PR wants to fix expectation values for PEPS/PEPO with dual physical spaces (relevant for fermionic networks). An additional test is included to ensure that the expectation value of identity operators is 1.

A side product is I generalized physicalspace(peps, r, c) etc to arbitrary integer values of r, c.

(Right now I only added the test, and haven't fix the twists yet. I want to finish this before #274 to make sure that I made the correct twists for the physical legs.)

[github-actions]

Your PR no longer requires formatting changes. Thank you for your contribution!

[Yue-Zhengyuan]

One particular thing I'm confused about is the normalization of reduced density matrices with str (supertrace) which can be negative. @lkdvos Can you explain the reason? @leburgel What do you think about this?

[Yue-Zhengyuan]

I vaguely remember this being related to the anti-periodic boundary condition in the imaginary time for fermions in the case of tr(exp(-βH)). But for other cases I'm still confused...

[lkdvos]

Looking at this, I'm actually completely unsure why we had a supertrace and why that was working in the first place. I would have to go through and see what we used to do, but I definitely think that when I implemented this I just tried some stuff until it worked.

I think the main point though is that the density matrix we are supplying to these functions actually contains precisely the minus sign that is canceled by the str, which is presumably because of how we are contracting the peps into rho to begin with. I think in general if you start with @tensor, you probably have to stick to that, although that is again just a hunch.

The main problem is that we at some point got a little "lazy" and just played around with the twists until things worked, but unfortunately these twists can move around a network (twist(t, 1) is equivalent to twisting all the other legs), so it's not clear to me that we have the most logical choice, nor that it works for any choice of arrows.

In general though, there is the additional question of what we mean with a dual physical space as well. One way we could define this is to start from a regular space and then flip it, which is a local basis transformation, but this has to be performed on the operators as well, and I think this transforms the identity into the exact operator to make the example you linked in the other PR work. An alternative view is to really take this as a space of itself, so then the identity is the really the identity, but I think that viewpoint is probably not consistent with the whole @tensor approach. I would have to think this through a bit more to give a full explanation of this, but maybe a different question is if we actually need this?

[Yue-Zhengyuan]

what we mean with a dual physical space

This is also kind of bothering me, which originated from our choice to purify a mixed state with dual ancilla spaces. (On the other hand, in add_physical_charge, we are fusing dual auxiliary spaces into the domain of the LocalOperator. This implies that we are fusing non-dual spaces to the physical legs of the state to be acted on.)

So strictly speaking, the resulting "purified state" is an operator. We discussed this earlier, but I still don't quite get it... It appears that when calculating expectation values of 2-layer PEPO, we should use $tr(Oρ^2)$ instead of $\langle \rho|O|\rho \rangle$. Current (physical) tests are all bosonic, and may not catch the errors.

[Yue-Zhengyuan]

I constructed a minimal example that may help understand what happened, but I can't precisely pin it down.

Consider a very simple state on 2 sites.

using Random, TensorKit, PEPSKit, LinearAlgebra
Random.seed!(2872346)
V1 = Vect[FermionParity](0 => 2, 1 => 2)
V2 = V1
V = V1 ⊗ V2
ψ = randn(ComplexF64, V)

Suppose we want to measure the value of an operator A1 acting on site1. The straightforward way to do it is:

A2 = TensorKit.id(V2);
nume = @plansor conj(ψ[b1 b2]) * (A1 ⊗ A2)[b1 b2; k1 k2] * ψ[k1 k2]
deno = @plansor conj(ψ[1 2]) * ψ[1 2]
val = nume / deno
display([nume, deno, val])
#= output
3-element Vector{ComplexF64}:
  5.471371195603439 + 4.618427854927436im
  6.927824745720233 + 0.0im
 0.7897675527925934 + 0.666649060050276im
=#

Here we use @plansor because we are evaluating the action of A1 on ψ, and then the inner product of ψ and (A1 ψ). @tensor only works when both V1, V2 are non-dual spaces.

Alternatively, we first construct the reduced density matrix on site 1: $\rho_1 = \mathrm{tr}_2{|\psi\rangle}{\langle\psi|}$. The correct way to calculate the partial trace is using @plansor:

ρ12 = ψ * ψ'
@plansor ρ1[k1; b1] := ρ12[k1 2; b1 2];
# equivalent construction
@plansor ρ1′[-1; -2] := conj(ψ[-2 b2]) * A2[b2; k2] * ψ[-1 k2]
@assert ρ1 ≈ ρ1′

Its eigenvalues are all positive as expected, and gives the correct expectation value of A1 using tr:

D, U = eigh(ρ1)
@assert minimum(D.data) > 0
@assert nume ≈ tr(A1 * ρ1)
@assert deno ≈ tr(ρ1)

What happened in PEPSKit might be something like the following: since contraction with the CTMRGEnv is done with @tensor, one might naively think that ρ1 is constructed with

@tensor ρ1_[-1; -2] := conj(ψ[-2 b2]) * A2[b2; k2] * ψ[-1 k2] # wrong

However, this introduced an unwanted twist that happens to be cancelled by str when evaluating expectation values. In addition, the eigenvalues of this wrong ρ1_ in the FermionParity(1) sector becomes negative.