Thanks a lot. I’m using C++ version of iTensor. Now I’ve come with a function like that:

```
complex<double> Correlation(MPS &psi, string opname1, int i, string opname2, int j, const SiteSet &sites)
{
// if (j <= i)
// cout << "Error in Correlation: i=" << i << " j=" << j << endl, exit(0);
ITensor op_i, op_j;
for (int o = 1; o <= 2; o++)
{
ITensor &op_s = (o == 1) ? op_i : op_j;
string name_s = (o == 1) ? opname1 : opname2;
int s = (o == 1) ? i : j;
if (name_s == "P") // Pair annihilation operator
{
op_s = op(sites, "Cup", s);
op_s *= prime(op(sites, "Cdn", s));
op_s.mapPrime(2, 1);
}
else if (name_s == "Pdag") // Pair creation operator
{
op_s = op(sites, "Cdagdn", s);
op_s *= prime(op(sites, "Cdagup", s));
op_s.mapPrime(2, 1);
}
else //Other single-site operator defined in ITEnsor for the FermionSite local Hilbert space
{
op_s = op(sites, name_s, s);
}
}
// Are the two operators "fermionic" (i.e. anticommmute) ?
const bool IsFermionic1 = (opname1 == "Cdagup" || opname1 == "Cdagdn" || opname1 == "Cup" || opname1 == "Cdn");
const bool IsFermionic2 = (opname2 == "Cdagup" || opname2 == "Cdagdn" || opname2 == "Cup" || opname2 == "Cdn");
if (IsFermionic1 != IsFermionic2)
cout << "Error: mixing a fermionic and a non-fermion operator." << endl, exit(1);
psi.position(i); //'gauge' the MPS to site i
//Create the <psi| ('bra') from |psi> ('ket')
auto psidag = dag(psi);
//Prime the link indices to make them distinct from
//the original ket links
psidag.prime("Link");
//below we will assume j > i
//index linking i-1 to i:
auto li_1 = leftLinkIndex(psi, i);
auto C = prime(psi(i), li_1) * op_i;
C *= prime(psidag(i), "Site");
// if (i!=j){
if (IsFermionic1)
{
for (int k = i + 1; k < j; ++k)
{
C *= psi(k);
C *= op(sites, "F", k); // Jordan-Wigner fermion 'string' operator F =(−1)^{n}
C *= prime(psidag(k), "Site");
}
}
else
{
for (int k = i + 1; k < j; ++k)
{
C *= psi(k);
C *= psidag(k);
}
}
// }
//index linking j to j+1:
auto lj = rightLinkIndex(psi, j);
C *= prime(psi(j), lj) * op_j;
C *= prime(psidag(j), "Site");
auto result = eltC(C);
return result;
}
```

Did I get that right that for c^{\dag}_{i,\uparrow}c_{i,\uparrow} I should just calculate

```
psi.position(i);
auto psidag = dag(psi);
psidag.prime("Link");
auto li_1 = leftLinkIndex(psi, i);
auto C = prime(psi(i), li_1) * op_i;
C *= prime(psidag(i), "Site");
auto lj = rightLinkIndex(psi, j);
C *= prime(psi(j), lj) * op_j;
C *= prime(psidag(j), "Site");
auto result = eltC(C);
```

right? For `i = j`

it seems like there’s no need in Jordan-Wigner “strings”. The same thing should apply for “pair” creation/annihilation operators c^{\dag}_{i,\uparrow}c^{\dag}_{I,\downarrow} (pair creation) and c_{i,\uparrow}c_{I,\downarrow}(pair annihilation) as pairs are bosons.

If so, I’m getting an error:

`Wrong number of IndexVals passed to elt/eltC (expected 0, got 2)`