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!