How to calculate 2-particle density matrix for the Hubbard Model when two indexes are equal?

Question

I'm studying the Hubbard Model and I want to calculate the following expected value:

$$\left< c_{i,\sigma_{i}}^{\dagger}c_{j,\sigma_{j}}^{\dagger}c_{k,\sigma_{k}}c_{l,\sigma_{l}}\right>,$$

which represents the two-particle density matrix by setting $\left(i,\sigma_{i}\ j,\sigma_{j}\right) \rightarrow \text{row; and }\left(k,\sigma_{k}\ l,\sigma_{l}\right) \rightarrow \text{column}$. I'm doing in the most optimized way, following -Correlator- and transforming the spinfull fermion operators into bosonic ones -Jordan Wigner-. I rewrite the Jordan Wigner to compress the different spin cases:

$$c_{i,\sigma_{i}} = F_{1}\cdots F_{i-1}\left(F_{i}^{\sigma_{i}}a_{i,\sigma_{i}}\right),$$

with $\sigma_{i} = 0\text{ corresponds } \uparrow \text{ and } \sigma_{i} = 1\text{ corresponds } \downarrow$. When $i\neq j\neq k\neq l$ the algorithm works flawlessly, but when two or more indexes are equal It doesn't.

Here's an example when (i==k) < j < l. I have to sort all operators in an increasing way. Using anti-commutator of fermions:

$$\left< c_{i,\sigma_{i}}^{\dagger}c_{j,\sigma_{j}}^{\dagger}c_{k,\sigma_{k}}c_{l,\sigma_{l}}\right> = -\left< c_{i,\sigma_{i}}^{\dagger}c_{k,\sigma_{k}}c_{j,\sigma_{j}}^{\dagger}c_{l,\sigma_{l}}\right>.$$

Substituting the transformation of Jordan Wigner we have:

$$\begin{align} =-\left( a_{i,\sigma_{i}}^{\dagger}F_{i}^{\sigma_{i}} \right) \underbrace{\left( F_{i-1}\cdots F_{1} \right)\left( F_{1}\cdots F_{k-1} \right)}_{=1} \left( F_{k}^{\sigma_{k}}a_{k,\sigma_{k}} \right) \left( a_{j,\sigma_{j}}^{\dagger}F_{j}^{\sigma_{j}} \right) \underbrace{\left( F_{j-1}\cdots F_{1} \right) \left( F_{1}\cdots F_{l-1} \right)}_{\left(F_{j}\cdots F_{l-1}\right)} \left( F_{l}^{\sigma_{l}}a_{l,\sigma_{l}} \right) \end{align}$$

$$=-\left(a_{i,\sigma_{i}}^{\dagger}F_{i}^{\sigma_{i}+\sigma_{k}}a_{k,\sigma_{k}}\right)\left(a_{j,\sigma_{j}}F_{j}^{\sigma_{j}+1}\right)\left(F_{j+1}\cdots F_{l-1}\right)\left(F_{l}^{\sigma_{l}}a_{l,\sigma_{l}}\right).$$

But I'm not getting the correct results. Here's a sample code:

IQIndex il = commonIndex(psi.A(i),psi.A(i+1),Link);
IQTensor C = psi.A(i);

C  *= (sigmai) ? sites.op("Adagdn",i) : sites.op("Adagup",i);
if ((sigmai+sigmak+1) % 2) {
    C.noprime(Site);
    C *= sites.op("F",i);
}
C.noprime(Site);
C *= (sigmak) ? sites.op("Adn",i) : sites.op("Aup",i);
C *= dag(prime(prime(psi.A(i),Site),il));

for(int m = i+1; m < j; ++m) {
    C *= psi.A(m);
    C *= dag(prime(psi.A(m),Link));
}

C *= psi.A(j);
C *= (sigmaj) ? sites.op("Adagdn",j) : sites.op("Adagup",j);
if (!sigmaj) {
    C.noprime(Site);
    C *= sites.op("F",j);
}
C *= dag(prime(prime(psi.A(j),Site),Link));

for(int m = j+1; m < l; ++m) {
    C *= psi.A(m);
    C *= sites.op("F",m);
    C *= dag(prime(prime(psi.A(m),Site),Link));
}

IQIndex kl = commonIndex(psi.A(l),psi.A(l-1),Link);
C *= psi.A(l);
if (sigmal) {
    C *= sites.op("F",l);
    C.noprime(Site);
}
C *= (sigmal) ? sites.op("Adn",l) : sites.op("Aup",l);
C *= dag(prime(prime(psi.A(l),Site),kl));

double result =  -C.real();

I still do not know where I'm going wrong. Whether it is in the Jordan Wigner transformation for spinfull fermions, or if in some part of the code. Maybe you can help me with some reference to this calculation. I thank you all!

Answer 1

Well, after posting this problem I ended up solving it. I decided to keep the post here if other people do the same silly thing as me. Another reason to leave the post here is that people who want to make a similar calculation can see and even exchange ideas with me.

The calculation is right, but the order of the operators' application is wrong in the code. Simple error, but since the correlator is calculated starting from the lowest to the highest site index (usually), I did not realize that in the same site the application of operators is obviously in the anti lexicographic order. The order of application is like this:

$$\begin{align} =-\underbrace{\left(a_{i,\sigma_{i}}^{\dagger}F_{i}^{\sigma_{i}+\sigma_{k}}a_{k,\sigma_{k}}\right)}_{\leftarrow}\underbrace{\left(a_{j,\sigma_{j}}F_{j}^{\sigma_{j}+1}\right)}_{\leftarrow}\underbrace{\left(F_{j+1}\cdots F_{l-1}\right)}_{\rightarrow}\underbrace{\left(F_{l}^{\sigma_{l}}a_{l,\sigma_{l}}\right)}_{\leftarrow}. \end{align}$$

Comments

Glad to hear you solved it, and thanks for posting the solution. These things are tricky to get right. We are hoping to add some features to ITensor later this year to make working with fermions easier