[Feature Request/Question] LocalOperator interface for correlation functions

[zxm403089989]

It seems that the updated correlators module in the recent version of PEPSKit supports the calculation of long-range correlation functions correlators.jl . I have two questions regarding this:

1. Does the current implementation fully support fermions? Specifically, if the tensors utilize fZ2 symmetry, can the correlation functions be calculated correctly?
2. Is it possible to add support for calculating correlation functions using LocalOperator, similar to the interface/style of expectation_value?
   Currently, correlator_horizontal and correlator_vertical seem to automatically convert input operators into MPOs. I believe adding direct support for LocalOperator (without explicit MPO conversion) would be more convenient and intuitive for calculating correlation functions in iPEPS.

[Yue-Zhengyuan]

1. We have support for fermions.

2. You can calculate correlation functions using expectation_value (check out our tests on correlator), but it is not contracting in the optimal way. I guess you are trying to calculate $c^\dagger c$ (Green's function). The reason we need to convert it to MPO is because each single c does not preserve fermion parity, which is not allowed by TensorKit (but $c^\dagger c$ preserves fermion parity). The virtual leg of the MPO may have a small bond dimension that will not increase the cost of evaluating the correlator.

[zxm403089989]

Thank you for your response and comments. I noticed that in the example you mentioned, it seems expectation_value is not restricted to calculating expectation values only within a single unit cell? From my reading of the code, I had understood that the results from expectation_value are actually expectation values within a single unit cell, because when calculating the reduced density matrix, coordinates beyond the unit cell range are taken mod to wrap back into the unit cell. Therefore, the result is not a true long-range correlation function, and correlator is what computes the true one. Or probably I misunderstanding on this part? Additionally, I'm trying to use ⊗ to define a local two-body operator for use with FiniteMPO, since correlator requires the MPO to act on two sites. My code is as follows and the returns are close to expectation_value of N. Is this approach reasonable?

let i = CartesianIndex(1, 1),js = [CartesianIndex(1, 3)]
    N=e_number(ComplexF64, U1Irrep, Trivial)
    pspace = space(N, 1)
    unit = TensorKit.id(pspace)
    MPON = (N ⊗ unit)
    PEPSKit.correlator_horizontal(peps,MPON,i,js,peps,env)
end

Here I assume that variables including peps, env are already evaluated and function e_number comes from MPSKitModels.jl

[lkdvos]

The expectation_value function is indeed not restricted to a single unit cell.
The mod coordinates are in fact correct, since the ansatz repeats the tensors, but if you inspect it carefully the patch that will be taken into account is different (larger) for reduced density matrices that extend beyond one unit cell.

In fact, the implementation for correlator, especially in the way that you are suggesting here, with a single i and j value should return the exact same result as the expectation_value. The difference is that expectation_value will take the entire patch into account and contract it, which might fail for large separations because at some point the contraction sequence finder will start taking prohibitively long.
On the other hand, the correlator hardcodes a contraction along a horizontal or vertical, being more restricted, but we can always find quasi-optimal contraction orders for that.

The biggest difference however is that correlator has the functionality of evaluating multiple expectation values in one go along a single line, which can greatly speed up the result since it doesn't have to recompute the partial contractions.

The approach you suggest seems reasonable, although I would recommend to use the correlator function instead of the correlator_horizontal. The result will be the same, but correlator is an exported and therefore public part of our API, while correlator_horizontal isn't, meaning we could at any point decide to change this implementation without it constituting a breaking change, which would break your code. You are of course free to use internal functions, but this is "use at own risk" 😉

[zxm403089989]

@Yue-Zhengyuan

1. We have support for fermions.

2. You can calculate correlation functions using expectation_value (check out our tests on correlator), but it is not contracting in the optimal way. I guess you are trying to calculate

   c
   †
   

   c
   (Green's function). The reason we need to convert it to MPO is because each single c does not preserve fermion parity, which is not allowed by TensorKit (but

   c
   †
   

   c
   preserves fermion parity). The virtual leg of the MPO may have a small bond dimension that will not increase the cost of evaluating the correlator.

Initially, I thought applying fermionic creation/annihilation operators to a PEPS only required naively contracting the auxiliary indices with those of a unit operator. However, after seeing issue 160, I realized that handling a single fermionic operator is not so straightforward—even in 1D, one must account for the Jordan-Wigner string. How can this type of problem be generalized to higher dimensions? Do you have any suggestions?

[Yue-Zhengyuan]

Honestly, I don't know about a general answer. Even for finite PEPS, there will be many different ways to drag the auxiliary leg to the boundary due to the 2D geometry. If the goal is to evaluate the expectation value of two c's, I still suggest trying to avoid dealing with a single c, unless you really need the excited state c+|psi>.

[lkdvos]

Note that this ambiguity can never really show up in "physical results", since at the end you always want to extract a scalar quantity. Since scalars transform trivially, they must always either have trivial charge (even number of fermionic operators) or be equal to 0.

The more difficult part is what to do with intermediary results, where you may have temporarily split off a part of the calculation which contains an odd fermionic charge, but will then connect it back again with another odd charged intermediate result. For these cases it is important to keep everything consistent, but that mostly means to just couple it back to the total calculation with an even charge and making sure to not handle one of the parts differently from the other. (Intuitively, you just have to make sure to keep the auxiliary string connected correctly)

[zxm403089989]

@lkdvos
Thank for reply. I’m not entirely sure if I have correctly understood your comment. Regarding the latter part, are you suggesting that as long as the auxiliary leg is always connected to another charged auxiliary leg, the fermionic even parity can be ensured and thus yield the correct result? If this is the case, then perhaps it is not strictly necessary, as described in the above issue, to pull the auxiliary index all locating the boundary.

[lkdvos]

I am indeed suggesting that the auxiliary leg is always connected to another auxiliary leg, but I'm also suggesting that this is a consequence of the physics, not a restriction of the implementation.
The point is that any scalar you can compute has to be even (or is forced to be zero), which is more or less the definition of a symmetry.
Therefore the only charged objects you can obtain are only an intermediate result, which will necessarily be combined with another charged object in the end to obtain a physical result.
As a concrete example, you could think of $$\braket{\psi |c^\dagger c | \psi}$$ or something similar, where the total is even, the $$\ket{\psi}$$ are even, and the $$c^\dagger c$$ are even, but we can choose to compute it as $$(\bra{\psi} c^\dagger) (c \ket{psi})$$, where now the expressions in the parentheses are odd.
In order to do so, we had to split the $$c^\dagger c$$ using some kind of factorization, which still has a string connecting the operators.
It is then perfectly fine to compute the odd subexpressions separately, but these will each have a dangling leg, and we must take care to connect these back together correctly.

For MPSKit, to handle this consistently, we choose the convention that we always migrate that leg to the boundary, mostly for ease of implementation reasons, as this works well with our datastructures.
In principle, there is nothing that prevents us from just keeping the dangling leg on the local MPS tensor and working with it like that though, but this would have required a bit more work to reimplement some of the contractions to account for that extra leg.
So indeed, it is absolutely true that it is not necessary to pull the auxiliary index to the boundary.

In the context of PEPSKit, since we are dealing with infinite PEPS, this wouldn't be possible to achieve, so the options are more limited.
The naive way would be to just have an extra leg on the PEPS tensor itself, which we could then either keep separate, or fuse into the physical leg, depending on what is more convenient for the code.
However, for expectation values and correlation functions, it turns out to be more convenient to not create the charged states at all, and compute the entire (even) expectation value instead, which is precisely what we is achieved by using this MPO -- connecting the string between the factorized $$c^\dagger c$$.

I hope this is slightly more clear, but feel free to ask more questions if anything is still hard to follow, some of it is quite subtle and I am not sure if I explain everything in the best way possible.

[zxm403089989]

@lkdvos

Thank you very much for your patient and detailed explanation. I believe I now basically understand your points. As you mentioned, the physical result must be charge 0, which follows from the definition of symmetry, just as described in the TensorKit tutorial:
https://quantumkithub.github.io/TensorKit.jl/latest/appendix/symmetric_tutorial/#Note-47d09bd0c202ac3d
https://quantumkithub.github.io/TensorKit.jl/latest/appendix/symmetric_tutorial/#Creation-and-annihilation-operators-as-symmetric-tensors
In both $Z_2$ and $U(1)$ symmetric systems, when defining a local operator, the symmetrically varying operator must introduce a nontrivial charged auxiliary index.
Unfortunately, I have not yet found that PEPSKit provides iPEPS computation tools capable of handling auxiliary indices. I mean modules that can be used directly for expectation values, propagators, or time evolution, in the same way as for ordinary iPEPS without auxiliary indices. Therefore, I speculate that if one really needs to compute iPEPS with auxiliary indices without adding new functions, a temporary approach could be to merge the auxiliary index with one of the physical indices according to certain conventions, thereby defining a charged iPEPS.

In my view, contracting a creation/annihilation operator with a “identity operator” carrying an auxiliary index at a neighboring site is essentially equivalent to merging the auxiliary index with the physical index of a spatial neighbor. And if we only need to compute the local fermion density, i.e., $(\left\langle \psi\right|c^{\dagger})(c\left|\psi\right\rangle )$, it might even be possible to fuse the outer index of the creation/annihilation operator with the auxiliary index (assuming the annihilation operator is defined as TensorMap(..., A ⊗ V ← V)), by adopting some method to fuse the codomain A and V.

However, for more complex fermionic correlator (suppose it is necessary to evaluate charged objects separately, for instance in the calculation of unequal-time propagators), for the convenience of computation, it seems necessary to ensure that when contracting two charged objects the auxiliary indices are located at the same lattice site, which would require using a string of parity-flip operators to move the auxiliary index (like a Jordan-Wigner string). The above reflects my personal understanding and could be mistaken.

As a minor question, does TensorKit or PEPSKit offer a function to fuse two indices of a TensorMap, or perhaps something analogous to dropdims in Base?

[Yue-Zhengyuan]

a temporary approach could be to merge the auxiliary index with one of the physical indices according to certain conventions, thereby defining a charged iPEPS.

You can take a look at the add_physical_charge function in PEPSKit/MPSKit which does exactly this.

does TensorKit or PEPSKit offer a function to fuse two indices of a TensorMap, or perhaps something analogous to dropdims in Base?

To fuse two spaces V1, V2, we construct an isomorphism(fuse(V1, V2) ← V1 ⊗ V2), and contract (multiply) it with the tensors whose legs you want to fuse.