I have issues, calculation correlation functions of 2-site operators in different ways. Before posting a lot of specific code, let me pls. ask in a more general context:

- Let \(|\psi\rangle\) be an MPS of finite even length L
- Let \(A_i B_{i+1}\) denote a nearest neighbor 2-site operator, i.e. some
`npc.Array shape=(2, 2, 2, 2) labels=['p0', 'p0*', 'p1', 'p1*']`

, which is composed from two 1-site operators \(\), and \(\) where - to simplify matter - both, \(A\) and \(B\) do also have string labels \(Astr\) and \(Bstr\) (Eg. one could just have Astr=Bstr='Sz' on some spin model). - To not have any questions about canonical forms later, let the 2-site op be unitary. (I don't know if this is relevant, but anyway)
- Eval \( C_a = \langle\psi|A_{L/2} B_{L/2+1} A_i B_{i+1}|\psi\rangle \) i.e.
`Ca = psi.correlation_function([Astr,Bstr],[Astr,Bstr], [L/2])`

- Now, set \(|\phi\rangle = A_{L/2} B_{L/2+1} |\psi\rangle\), i.e.
`phi = copy.deepcopy(psi); phi.apply_local_op(L/2,A B)`

- Define an MPSEnvironment like
`bk = MPSEnvironment(phi,psi)`

and eval the expectation value \(C_b = \langle\phi|A_i B_{i+1}|\psi\rangle\) with`Cb = bk.expectation_value(AB)`

So what is it that I do not understand?

And, how to properly eval correlation functions using MPSEnvironment?

PS.: If I play the same game with only a single site operator, then \(C_a = C_b\) as I expected.