Extend `correlator` to pair correlation

[ZongYongyue]

Hi Lukas,

I'm trying to extend the function of correlator to pair correlation like $$\langle c^\dagger_i c^\dagger_j c_k c_l\rangle$$. I'm doing two things for this purpose: one is adding new method for transfer_left:

function transfer_left(v::AbstractTensorMap{T, S, 4, 1}, A::MPSTensor{S}, Ab::MPSTensor{S}) where {T, S}
    check_unambiguous_braiding(space(v, 4))
    @plansor v[-1 -2 -3 -4; -5] := v[1 2 3 4; 5] * A[5 6; -5] * τ[4 7; 6 -4] * τ[3 8; 7 -3] * τ[2 9; 8 -2] * conj(Ab[1 9; -1])
end

function transfer_left(v::AbstractTensorMap{T, S, 3, 1}, A::MPSTensor{S}, Ab::MPSTensor{S}) where {T, S}
    check_unambiguous_braiding(space(v, 2))
    @plansor v[-1 -2 -3; -4] := v[1 2 3; 4] * A[4 5; -4] * τ[3 6; 5 -3] * τ[2 7; 6 -2] * conj(Ab[1 7; -1])
end

to deal with the transfer matrix like this:

        @plansor Vᵢ[-1 -2; -3] := state.AC[i][3 4; -3] * conj(U[1]) * I[1 2; 4 -2] *
                            conj(state.AC[i][3 2; -1])
        if k > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(k - 1)])
        end                     
        @plansor Vₖ[-1 -2 -3 -4; -5] := Vᵢ[1 2; 3] * state.AR[k][3 4; -5] * 
                            K[5 6; 4 -4] * τ[5 7; 6 -3] * τ[2 8; 7 -2] * conj(state.AR[k][1 8; -1])
        if j > (k + 1)
             Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(j - 1)])
        end

It works, but I am not sure if I write it in a right way and I don't know the means of this functioncheck_unambiguous_braiding here.

The another thing is to contract different diagrams for different orders of i, j ,k ,l , like i<j<k<l,i<j=k<l, i=k<l<j ... I have successfully done some of them, but I encounter a problem with this diagram:

      @plansor Vₖ[-1 -2 -3; -4] := state.AC[k][1 2; -4] * K[3 4; 2 -3] *
                          τ[3 5; 4 -2] * conj(state.AC[k][1 5; -1])
  #     if i > (k + 1)
  #         Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(i - 1)])
  #     end
      @plansor Vᵢ[-1 -2 -3 -4; -5] := Vₖ[1 2 3; 4] * state.AR[i][4 5; -5] * τ[3 6; 5 -4] * 
                          τ[2 7; 6 -3] * conj(U[8]) * I[8 9; 7 -2] * conj(state.AR[i][1 9; -1])

It returns ERROR: LoadError: ArgumentError: not a planar diagram expression: ((var"##state.AR[i]#1082"[4 5; -5] * (var"##τ#1091"[(3, 6); (5, -4)] * ((Vₖ[1 2 3; 4] * var"##state.AR[i]#1082"'[(-1,); (1, 9)]) * var"##τ#1092"[(2, 7); (6, -3)]))) * I[8 9; 7 -2]) * U'[(); (8,)]

I'm not sure where is wrong, and it is a little weird because I follow the same thought and contract other diagrams successfully like:

@plansor Vᵢ[-1 -2; -3] := state.AC[i][3 4; -3] * conj(U[1]) * I[1 2; 4 -2] *
                  conj(state.AC[i][3 2; -1])
if k > (i + 1)
  Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(k - 1)])
end                     
@plansor Vₖ[-1 -2 -3 -4; -5] := Vᵢ[1 2; 3] * state.AR[k][3 4; -5] * 
                  K[5 6; 4 -4] * τ[5 7; 6 -3] * τ[2 8; 7 -2] * conj(state.AR[k][1 8; -1])
if l > (k + 1)
  Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(l - 1)])
end
@plansor Vₗ[-1 -2 -3; -4] := Vₖ[1 2 3 4; 5] * state.AR[l][5 6; -4] * L[4 7; 6 10] * U[10] *
                  τ[3 8; 7 -3] * τ[2 9; 8 -2] * conj(state.AR[l][1 9; -1])
if j > (l + 1)
  Vₗ = Vₗ * TransferMatrix(state.AR[(l + 1):(j - 1)])
end
G = @plansor Vₗ[1 2 3; 4] * state.AR[j][4 5; 9] *
                  τ[3 6; 5 7] * J[2 8; 6 7] * conj(state.AR[j][1 8; 9])

[lkdvos]

Thanks for opening this, this would indeed be a nice feature to have.

The check_unambiguous_braiding function is there mostly for working with non-symmetric braiding rules (such as for anyonic systems).
The point is that in principle, MPSKit consists entirely of planar diagrams, and therefore supports these anyonic systems.
However, there are a few places where we do actually have crossings that show up, often because we use auxiliary trivial legs to uniformize the code.
For these crossings, the choice of over- or under-braiding can make a difference in the result, and given that we make an arbitrary choice here we put a check in place to make sure that we only do this whenever the result is unambiguous (for example when braiding trivial legs, over and undercrossings are the same).
Long story short: when we dont or cannot distinguish between using τ or adjoint(τ) here we throw an error to avoid silently returning wrong answers.

For the diagrams, there is a subtle interplay between contraction order and planarity of the diagram.
First note that when using @plansor with integers as labels, these are interpreted as indicating the order in which this network will be contracted.
In particular, that means that in the diagram that first label 1, then label 2 etc will be contracted.
Secondly, there is one exception to this: whenever a label is selected to be contracted over, all connecting legs between these two contractions will be contracted, independent of their label.
For example, in the diagram below first label 1 is selected, so C is contracted with A, but since they also share label 3 that is immediately contracted as well. The final order of contraction would therefore be: (1, 3), (2, 4) where I bundled together contractions that happen simultaneously.

@tensor C[1 2 3] * A[1 3 4] * B[2 4]

Given these two points, in the contraction you gave, the last contractions would be I with the remainder of the network, followed by U with that (note that 9 is contracted before 8 because it will happen when 7 is selected). This cannot happen in a planar way because of the open leg with label 8, which would become enclosed by contracting 7 and 9.
From what I can tell, simply swapping 7 and 8 should do the trick. (You might want to double check if the contraction orders are optimal though, see eg the costcheck=:warn keyword for the @tensor macro, which should also work here.
https://jutho.github.io/TensorOperations.jl/latest/man/indexnotation/#TensorOperations.@tensor

As a final note, it might be easier to simply not deal with the cases where ijkl are unsorted, and mapping them to sorted cases instead by permuting the operator. Taking c_i c_j c_k c_l as a 4,4 tensor map, you should be able to permute the tensor to sort the values first, and then apply the correlator implementation for that new operator. For the overlapping cases it requires a bit more care since you'd need to turn these into 3,3 or 2,2 operators first, but that might still be easier than figuring out all of these contractions with different connections going back and forth.

[ZongYongyue]

Thank you for your detailed explanation and very helpful teaching, it resolved my confusion!

[ZongYongyue]

I've more or less finished the pair correlation for correlator, and at least in my use cases, it works. I contracted different diagrams for different order of i, j, k, l instead of first contracting c_i c_j c_k c_l into a (4,4) tensor and then permuting, because I felt the latter adds an unnecessary contraction and decomposition step for the user. I did not checked whether the current contraction order for each diagram is optimal, because I'm in a hurry to use it now. I won't pull request because it is so rough. I share the code here in case it might be helpful to others or as a reminder for developers.

"""
    correlator(state::AbstractFiniteMPS, O₁::AbstractTensorMap, O₂::AbstractTensorMap, i, j, k, l)
"""
function correlator(state::AbstractFiniteMPS, O₁::AbstractTensorMap, O₂::AbstractTensorMap, i, j, k, l)
        I, J = decompose_localmpo(O₁)
        K, L = decompose_localmpo(O₂)
        U = ones(scalartype(state), _firstspace(O₁))
    if i < j < k < l
        @plansor Vᵢ[-1 -2; -3] := state.AC[i][3 4; -3] * conj(U[1]) * I[1 2; 4 -2] *
                            conj(state.AC[i][3 2; -1])
        if j > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(j - 1)])
        end
        @plansor Vⱼ[-1 -2; -3] := Vᵢ[1 2; 3] * state.AR[j][3 4; -3] * J[2 5; 4 -2] *
                            conj(state.AR[j][1 5; -1])
        if k > (j + 1)
            Vⱼ = Vⱼ * TransferMatrix(state.AR[(j + 1):(k - 1)])
        end
        @plansor Vₖ[-1 -2; -3] := Vⱼ[1 2; 3] * state.AR[k][3 4; -3] * K[2 5; 4 -2] *
                            conj(state.AR[k][1 5; -1])
        if l > (k + 1)
            Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(l - 1)])
        end
        G = @plansor Vₖ[2 3; 5] * state.AR[l][5 6; 7] * L[3 4; 6 1] * U[1] *
                        conj(state.AR[l][2 4; 7])

    elseif i < j == k < l
        @plansor Vᵢ[-1 -2; -3] := state.AC[i][3 4; -3] * conj(U[1]) * I[1 2; 4 -2] *
                            conj(state.AC[i][3 2; -1])
        if j > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(j - 1)])
        end
        @plansor Vⱼ[-1 -2; -3] := Vᵢ[1 2; 3] * state.AR[j][3 4; -3] * K[5 6; 4 -2] *
                            τ[7 8; 5 6] * J[2 9; 7 8] * conj(state.AR[j][1 9; -1])
        if l > (j + 1)
            Vⱼ = Vⱼ * TransferMatrix(state.AR[(j + 1):(l - 1)])
        end
        G = @plansor Vⱼ[2 3; 5] * state.AR[l][5 6; 7] * L[3 4; 6 1] * U[1] *
                            conj(state.AR[l][2 4; 7])

    elseif i < k < j < l
        @plansor Vᵢ[-1 -2; -3] := state.AC[i][3 4; -3] * conj(U[1]) * I[1 2; 4 -2] *
                            conj(state.AC[i][3 2; -1])
        if k > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(k - 1)])
        end                     
        @plansor Vₖ[-1 -2 -3 -4; -5] := Vᵢ[1 2; 3] * state.AR[k][3 4; -5] * 
                            K[5 6; 4 -4] * τ[5 7; 6 -3] * τ[2 8; 7 -2] * conj(state.AR[k][1 8; -1])
        if j > (k + 1)
             Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(j - 1)])
        end
        @plansor Vⱼ[-1, -2; -3] := Vₖ[1 2 3 4; 5] * state.AR[j][5 6; -3] * τ[4 7; 6 -2] * 
                            τ[3 8; 7 9] * J[2 10; 8 9] * conj(state.AR[j][1 10; -1])
        if l > (j + 1)
             Vⱼ = Vⱼ * TransferMatrix(state.AR[(j + 1):(l - 1)])
        end
        G = @plansor Vⱼ[2 3; 5] * state.AR[l][5 6; 7] * L[3 4; 6 1] * U[1] *
                            conj(state.AR[l][2 4; 7])
    
    elseif i < k < j == l
        @plansor Vᵢ[-1 -2; -3] := state.AC[i][3 4; -3] * conj(U[1]) * I[1 2; 4 -2] *
                            conj(state.AC[i][3 2; -1])
        if k > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(k - 1)])
        end                     
        @plansor Vₖ[-1 -2 -3 -4; -5] := Vᵢ[1 2; 3] * state.AR[k][3 4; -5] * 
                            K[5 6; 4 -4] * τ[5 7; 6 -3] * τ[2 8; 7 -2] * conj(state.AR[k][1 8; -1])
        if j > (k + 1)
             Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(j - 1)])
        end
        G = @plansor Vₖ[1 2 3 4; 5] * state.AR[j][5 6; 11] * L[4 7; 6 12] * U[12] *
                            τ[3 8; 7 9] * J[2 10; 8 9] * conj(state.AR[j][1 10; 11])

    elseif i == k < j < l
        @plansor Vᵢ[-1 -2 -3 -4; -5] := state.AC[i][3 7; -5] * K[5 6; 7 -4] * τ[5 4; 6 -3] *
                            conj(U[1]) * I[1 2; 4 -2] * conj(state.AC[i][3 2; -1])
        if j > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(j - 1)])
        end
        @plansor Vⱼ[-1 -2; -3] := Vᵢ[1 2 3 4; 5] * state.AR[j][5 6; -3] * τ[4 7; 6 -2] *
                            τ[3 8; 7 9] * J[2 10; 8 9] * conj(state.AR[j][1 10; -1])
        if l > (j + 1)
            Vⱼ = Vⱼ * TransferMatrix(state.AR[(j + 1):(l - 1)])
        end
        G = @plansor Vⱼ[2 3; 5] * state.AR[l][5 6; 7] * L[3 4; 6 1] * U[1] *
                            conj(state.AR[l][2 4; 7])
    
    elseif i == k < j == l
        @plansor Vᵢ[-1 -2 -3 -4; -5] := state.AC[i][3 7; -5] * K[5 6; 7 -4] * τ[5 4; 6 -3] *
                            conj(U[1]) * I[1 2; 4 -2] * conj(state.AC[i][3 2; -1])
        if j > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(j - 1)])
        end
        G = @plansor Vᵢ[1 2 3 4; 5] * state.AR[j][5 6; 11] * L[4 7; 6 12] * U[12] * 
                            τ[3 8; 7 9] * J[2 10; 8 9] * conj(state.AR[j][1 10; 11])
    
    elseif i < k < l < j
        @plansor Vᵢ[-1 -2; -3] := state.AC[i][3 4; -3] * conj(U[1]) * I[1 2; 4 -2] *
                            conj(state.AC[i][3 2; -1])
        if k > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(k - 1)])
        end                     
        @plansor Vₖ[-1 -2 -3 -4; -5] := Vᵢ[1 2; 3] * state.AR[k][3 4; -5] * 
                            K[5 6; 4 -4] * τ[5 7; 6 -3] * τ[2 8; 7 -2] * conj(state.AR[k][1 8; -1])
        if l > (k + 1)
            Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(l - 1)])
        end
        @plansor Vₗ[-1 -2 -3; -4] := Vₖ[1 2 3 4; 5] * state.AR[l][5 6; -4] * L[4 7; 6 10] * U[10] *
                            τ[3 8; 7 -3] * τ[2 9; 8 -2] * conj(state.AR[l][1 9; -1])
        if j > (l + 1)
            Vₗ = Vₗ * TransferMatrix(state.AR[(l + 1):(j - 1)])
        end
        G = @plansor Vₗ[1 2 3; 4] * state.AR[j][4 5; 9] *
                            τ[3 6; 5 7] * J[2 8; 6 7] * conj(state.AR[j][1 8; 9])

    elseif i == k < l < j
        @plansor Vᵢ[-1 -2 -3 -4; -5] := state.AC[i][3 7; -5] * K[5 6; 7 -4] * τ[5 4; 6 -3] *
                            conj(U[1]) * I[1 2; 4 -2] * conj(state.AC[i][3 2; -1])
        if l > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(l - 1)])
        end
        @plansor Vₗ[-1 -2 -3; -4] := Vᵢ[1 2 3 4; 5] * state.AR[l][5 6; -4] * L[4 7; 6 10] * U[10] *
                            τ[3 8; 7 -3] * τ[2 9; 8 -2] * conj(state.AR[l][1 9; -1])
        if j > (l + 1)
            Vₗ = Vₗ * TransferMatrix(state.AR[(l + 1):(j - 1)])
        end
        G = @plansor Vₗ[1 2 3; 4] * state.AR[j][4 5; 9] *
                            τ[3 6; 5 7] * J[2 8; 6 7] * conj(state.AR[j][1 8; 9])

    elseif k < i < j < l
        @plansor Vₖ[-1 -2 -3; -4] := state.AC[k][1 2; -4] * K[3 4; 2 -3] *
                            τ[3 5; 4 -2] * conj(state.AC[k][1 5; -1])
        if i > (k + 1)
            Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(i - 1)])
        end
        @plansor Vᵢ[-1 -2 -3 -4; -5] := Vₖ[1 2 3; 4] * state.AR[i][4 5; -5] * τ[3 6; 5 -4] * 
                            τ[2 8; 6 -3] * conj(U[7]) * I[7 9; 8 -2] * conj(state.AR[i][1 9; -1])

        if j > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(j - 1)])
        end
        @plansor Vⱼ[-1 -2; -3] := Vᵢ[1 2 3 4; 5] * state.AR[j][5 6; -3] * τ[4 7; 6 -2] * τ[3 8; 7 9] *
                            J[2 10; 8 9] * conj(state.AR[j][1 10; -1])
        if l > (j + 1)
            Vⱼ = Vⱼ * TransferMatrix(state.AR[(j + 1):(l - 1)])
        end
        G = @plansor Vⱼ[1 2; 3] * state.AR[l][3 4; 6] * L[2 5; 4 7] * U[7] * conj(state.AR[l][1 5; 6])

    elseif k < i < j == l
        @plansor Vₖ[-1 -2 -3; -4] := state.AC[k][1 2; -4] * K[3 4; 2 -3] *
                            τ[3 5; 4 -2] * conj(state.AC[k][1 5; -1])
        if i > (k + 1)
            Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(i - 1)])
        end
        @plansor Vᵢ[-1 -2 -3 -4; -5] := Vₖ[1 2 3; 4] * state.AR[i][4 5; -5] * τ[3 6; 5 -4] * 
                            τ[2 8; 6 -3] * conj(U[7]) * I[7 9; 8 -2] * conj(state.AR[i][1 9; -1])

        if j > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(j - 1)])
        end
        G = @plansor Vᵢ[1 2 3 4; 5] * state.AR[j][5 6; 12] * L[4 8; 6 7] * U[7] * τ[3 9; 8 10] *
                            J[2 11; 9 10] * conj(state.AR[j][1 11; 12])
        
    elseif k < i < l < j
        @plansor Vₖ[-1 -2 -3; -4] := state.AC[k][1 2; -4] * K[3 4; 2 -3] *
                            τ[3 5; 4 -2] * conj(state.AC[k][1 5; -1])
        if i > (k + 1)
            Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(i - 1)])
        end
        @plansor Vᵢ[-1 -2 -3 -4; -5] := Vₖ[1 2 3; 4] * state.AR[i][4 5; -5] * τ[3 6; 5 -4] * 
                            τ[2 8; 6 -3] * conj(U[7]) * I[7 9; 8 -2] * conj(state.AR[i][1 9; -1])
        if l > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(l - 1)])
        end
        @plansor Vₗ[-1 -2 -3; -4] := Vᵢ[1 2 3 4; 5] * state.AR[l][5 6; -4] * L[4 8; 6 7] * U[7] *
                            τ[3 9; 8 -3] * τ[2 10; 9 -2] * conj(state.AR[l][1 10; -1])
        if j > (l + 1)
            Vₗ = Vₗ * TransferMatrix(state.AR[(l + 1):(j - 1)])
        end
        G = @plansor Vₗ[1 2 3; 4] * state.AR[j][4 5; 9] * τ[3 6; 5 7] * J[2 8; 6 7] * 
                            conj(state.AR[j][1 8; 9])

    elseif k < i == l < j
        @plansor Vₖ[-1 -2 -3; -4] := state.AC[k][1 2; -4] * K[3 4; 2 -3] *
                            τ[3 5; 4 -2] * conj(state.AC[k][1 5; -1])
        if i > (k + 1)
            Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(i - 1)])
        end
        @plansor Vᵢ[-1 -2 -3; -4] := Vₖ[1 2 3; 4] * state.AR[i][4 5; -4] * L[3 7; 5 6] * U[6] *
                            τ[2 9; 7 -3] * I[8 10; 9 -2] * conj(U[8]) * conj(state.AR[i][1 10; -1])
        if j > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(j - 1)])
        end
        G = @plansor Vᵢ[1 2 3; 4] * state.AR[j][4 5; 6] *
                            τ[3 7; 5 8] * J[2 9; 7 8] * conj(state.AR[j][1 9; 6])

    elseif k < l < i < j
        @plansor Vₖ[-1 -2 -3; -4] := state.AC[k][1 2; -4] * K[3 4; 2 -3] *
                            τ[3 5; 4 -2] * conj(state.AC[k][1 5; -1])
        if l > (k + 1)
            Vₖ = Vₖ * TransferMatrix(state.AR[(k + 1):(l - 1)])
        end
        @plansor Vₗ[-1 -2; -3] := Vₖ[1 2 3; 4] * state.AR[l][4 6; -3] * L[3 7; 6 5] * 
                            U[5] * τ[2 8; 7 -2] * conj(state.AR[l][1 8; -1])
        if i > (l + 1)
                Vₗ = Vₗ * TransferMatrix(state.AR[(l + 1):(i - 1)])
        end
        @plansor Vᵢ[-1 -2 -3; -4] := Vₗ[1 2; 3] * state.AR[i][3 4; -4] * τ[2 6; 4 -3] * 
                            I[5 7; 6 -2] * conj(U[5]) * conj(state.AR[i][1 7; -1])
        if j > (i + 1)
            Vᵢ = Vᵢ * TransferMatrix(state.AR[(i + 1):(j - 1)])
        end
        G = @plansor Vᵢ[1 2 3; 4] * state.AR[j][4 5; 9] * τ[3 6; 5 7] * J[2 8; 6 7] * 
                            conj(state.AR[j][1 8; 9]) 
    else
        throw(ArgumentError("invalid input indices (i, j, k, l) for ($i, $j, $k, $l), only (i < j) && (k < l) is valid"))
    end
    return G
end

[lkdvos]

Thanks for sharing this! Any chance you can also share the transfer functions?

[ZongYongyue]

Sure!

function transfer_left(v::AbstractTensorMap{T, S, 4, 1}, A::MPSTensor{S}, Ab::MPSTensor{S}) where {T, S}
    @plansor v[-1 -2 -3 -4; -5] := v[1 2 3 4; 5] * A[5 6; -5] * τ[4 7; 6 -4] * τ[3 8; 7 -3] * τ[2 9; 8 -2] * conj(Ab[1 9; -1])
end

function transfer_left(v::AbstractTensorMap{T, S, 3, 1}, A::MPSTensor{S}, Ab::MPSTensor{S}) where {T, S}
    @plansor v[-1 -2 -3; -4] := v[1 2 3; 4] * A[4 5; -4] * τ[3 6; 5 -3] * τ[2 7; 6 -2] * conj(Ab[1 7; -1])
end

[lkdvos]

As a small remark here: This is not entirely the same as the correlation function, in the sense that specifically that function is meant to compute correlators for different values of j in a more efficient manner. For a single value of i and j, this should be completely equivalent to the expectation_value function.
In that sense, I would argue that this function is more a generalization of the expectation_value function to also deal with non-sorted inputs, I'm not sure if you agree with this?

For actually computing a pair correlation function for different values of i,j,k,l in an efficient way I'm not entirely sure how to tackle that, since my naive ideas to reuse part of the environments would involve having to store n^4 different Vl etc, which might not be feasible.
Ultimately I guess this isn't really ever a bottleneck, and you can just compute all of them separately (and in parallel?).
I'd be very happy to hear your ideas about this, any suggestions are always welcome.

[ZongYongyue]

I get used to call $\langle O_i \rangle$ expectation value of O and $\langle A_i B_j \rangle$ AB static correlation function, but I agree with you that there is no difference between them in the view of tensor contraction.

Actually, fix $c_ic_j$ and rung $c_kc_l$ in (ks, ls), crossover and overlap of i, j, k, l will occur only the two pairs are close to each other. For example the d wave pairing correlation $\Psi$(r) at r<3. Usually, I think that developers consider how to compute the order of i<j<k<l efficiently is enough. As to the other cases (the other orders), computing the expectation value of $\langle c_ic_jc_kc_l\rangle$ one by one is not too slow from my use cases.

Of course, this is just a consideration from a pragmatic point of view and may be somewhat one-sided — it's merely one perspective for reference.😄

[lkdvos]

I see, the naming might indeed be confusing.
In my head, the main difference was really about the implementation: we are trying to reuse part of the contractions in correlator, while expectation_value simply has independent evaluations.

I'll look into it to see if this can be added in some form that is useful to the main package too, I agree with you that this is a useful function too have.
Thank you for your perspective though, it is definitely useful to have input from users/other developers