[Julia] computing correlator using efficient MPS techniques

Question

Hi all,

I am trying to understand the code example for computing 2pt correlator of spinless fermion, using efficient MPS techniques. (the built-in function correlation_matrix can not be applied to mixed types of sites which I'm currently working on)

The codes are quoted below:

orthogonalize!(psi,1)
C = zeros(N,N)
L = ITensor(1.)
for i=1:N
Li = L*psi[i]*op("Cdag",s,i)*dag(prime(psi[i]))
LiF = Li
for j=i+1:N
    lind = commonind(psi[j],LiF)
    LiF *= psi[j]
    cij = LiF*op("C",s,j)*dag(prime(prime(psi[j],"Site"),lind))
    C[i,j] = scalar(cij)
    C[j,i] = conj(C[i,j])
    LiF *= op("F",s,j)*dag(prime(psi[j]))
  end
  L *= psi[i]*dag(prime(psi[i],"Link"))
end    

The only line makes me confused is LiF*=op("F",s,j)*dag(prime(psi[j])). According to the definition given here, op("F") counts the parity of the fermion number at a given site: if there is a fermion at that site, it gives an additional minus sign.

It must somehow relate to the anti-commutation relation of fermions, but why we need it, can someone give me some insight on this?

Best,
Junsen

Comments

BTW, unexpectedly, I find this implementation is much slower than the naive one (i.e., using OpSum to make an MPO for each i,j)

Here is my code for calculate the fermion 2pt correlator using this MPS techniques:


```
begin
    orthogonalize!(psi,1)
    corr1 = zeros(div(N,3),div(N,3))
    L = ITensor(1.)
    for i=1:N
        if mod(i,3)==0
            Li = L*psi[i]*op("Cdag",sites,i)*dag(prime(psi[i]))
            LiF = Li
            for j=i+3:3:N
                LiF *= psi[j-2]*op("F",sites,j-2)*dag(prime(psi[j-2]))
                LiF *= psi[j-1]*dag(prime(psi[j-1],"Link"))
                lind = commonind(psi[j],LiF)
                LiF *= psi[j]
                cij = LiF*op("C",sites,j)*dag(prime(prime(psi[j],"Site"),lind))
                corr1[div(i,3),div(j,3)] = scalar(cij)
                corr1[div(j,3),div(i,3)] = conj(corr1[div(i,3),div(j,3)])
                LiF *= op("F",sites,j)*dag(prime(psi[j]))
            end
        end
        L *= psi[i]*dag(prime(psi[i],"Link"))
    end
               
    # fill in diagonal elements
    for i=1:div(N,3)
        orthogonalize!(psi,3*i)
        corr1[i,i] = scalar(psi[3i]*op("N",sites,3*i)*dag(prime(psi[3i],"Site")))
    end
end
```

Note the sites are defined as follows (basically a fermion ladder with dynamical Z2 gauge fields is deformed into a fermion-spin-fermion chain):

```
sites = [ isodd(mod1(n,3)) ? siteind("Fermion",n,conserve_qns=true) : Index([QN()=>2];tags="S=1/2,Site,n=$n") for n in 1:N]
```

Why this method is much slower than a naive MPO construction from OpSum?

---

Hi JW,
So I ran your code just now for an MPS of 30 sites and bond dimension 10, and I found that your code is much faster than using OpSum to make the MPO's. Did you use OpSum to get the entire correlation matrix for all the same i and j values as in your code? What size MPS did you input (both number of sites and bond dimension)?

Miles

Answer 1

Hi, Junsen,

I'm a C++ user, so I hope everything translates easily enough to the Julia side. I also won't go through the Julia code too carefully or test it since I don't use the language at this point.

Why we need that extra negative sign is related to the normal ordering of the Fermionic operators and the anti-commutation properties of A/Adag and C/Cdag (not the same thing) when not acting on the same site.

Take some time to really make sure you understand the documentation on this if working with Fermions.
For the basic labels of the operators:
http://itensor.org/docs.cgi?vers=cppv3&page=classes/fermion
(When the Julia version is selected it throws an error).
And for their properties:
http://itensor.org/docs.cgi?page=tutorials/fermions&vers=julia
Note that the end of this emphasizes that AutoMPO will do this for you, but when wanting to measure properties, one needs to be careful to include or exclude those terms as desired. One of the first things a new user can do with the Fermions is create your favorite Fermionic model and measure "good" and "bad" correlations including or excluding the Fermionic phases from the anti-commutation relations.

For your question on OpSum, I'm not sure that I follow. Looking at the Github documentation https://itensor.github.io/ITensors.jl/stable/OpSum.html#ITensors.OpSum
The operators used there are sums of local operators (namely, writing down ampo like a local Hamiltonian).
The correlation function constructed in the code example uses multiplication of operators. That is, when we use the Jordan-Wigner transform we go from a local basis of Fermions, which are difficult to use directly in a computational setting to a non-local set of spins or bosons. These strings of operators are non-local and so I don't see how OpSum applies.

There will be a substantial difference in calculation time between global and local measurements. There is a nice diagrammatic introduction to these two-point operators for spins (but only in C++, sorry!):
https://www.itensor.org/docs.cgi?page=tutorials/correlations&vers=cppv3
In the first picture, site "3" is really just all the sites between i and j (exclusive). Coming back to Fermions, this is the region over which we need the phase contribution for each basis state within the MPS from the anticommutation relations. We cannot reduce this into the average fermionic phase contribution and do a sum multiplied by the average value. We need the string for each basis state in the MPS over that region to calculate the phase, then include that when doing the summation over basis states.

Best wishes,

Jared

Comments

Dear Jared,

Sorry for this delayed reply, and thank you very much for the detailed answer!

Recently (after finishing my previous unrelated work) I've taken some times to check these details on Jordan-Wigner string and its implementation in ITensor. It turns out that the doc
http://itensor.org/docs.cgi?page=tutorials/fermions&vers=cppv3
you recommend is really helpful. Let me present my own understanding:

The series of operators `op("F",s,j)` is literally just a Jordan-Wigner string one must add by hand when computing correlation functions involving fermions. This is a direct consequence of the fact that the fermion operator defined by ITensor uses Jordan-Wigner's approach, but WITHOUT a Jordan-Wigner string (except in `OpSum`).

Best,

Junsen