Proper way to construct fermion gates

[nGn]

I am trying to construct the following state of fermions:

\begin{align} |\psi\rangle &= \frac{1}{2} \Big( 1 + c_3^{\dagger} c_1^{\dagger} \Big) \Big( 1 + c_2^{\dagger} c_4^{\dagger} \Big) \, |0_1, 0_2, 0_3, 0_4 \rangle \\ &= \frac{1}{2} \Big( 1 - c_1^{\dagger} c_3^{\dagger} + c_2^{\dagger} c_4^{\dagger} + c_1^{\dagger} c_2^{\dagger} c_3^{\dagger} c_4^{\dagger} \Big) \, |0_1, 0_2, 0_3, 0_4 \rangle, \end{align}

which has a correlation matrix

\langle \psi | \begin{pmatrix} c_1^\dagger c_1^\dagger & c_1^\dagger c_2^\dagger & c_1^\dagger c_3^\dagger & c_1^\dagger c_4^\dagger \\ c_2^\dagger c_1^\dagger & c_2^\dagger c_2^\dagger & c_2^\dagger c_3^\dagger & c_2^\dagger c_4^\dagger \\ c_3^\dagger c_1^\dagger & c_3^\dagger c_2^\dagger & c_3^\dagger c_3^\dagger & c_3^\dagger c_4^\dagger \\ c_4^\dagger c_1^\dagger & c_4^\dagger c_2^\dagger & c_4^\dagger c_3^\dagger & c_4^\dagger c_4^\dagger \end{pmatrix} | \psi \rangle = \begin{pmatrix} 0 & 0 & -1/2 & 0 \\ 0 & 0 & 0 & 1/2 \\ 1/2 & 0 & 0 & 0 \\ 0 & 0 & -1/2 & 0 \end{pmatrix}

using ITensors

s = siteinds("Fermion", 4);
ψ0 = productMPS(s, "Emp");

ampo24 = AutoMPO();
ampo24 += (1/sqrt(2), "Id", 1); 
ampo24 += (1/sqrt(2), "c†", 2, "c†", 4);

ampo31 = AutoMPO();
ampo31 += (1/sqrt(2), "Id", 1);
ampo31 += (1/sqrt(2), "c†", 3, "c†", 1);

o24 = MPO(ampo24, s);
o31 = MPO(ampo31, s);

ψ = product(o31, product(o24, ψ0));
correlation_matrix(ψ, "c†", "c†")

This gives the correct result.
At the same time, I had expected (based on Matt’s response in [julia] TEBD gates involving fermions - ITensor Support Q&A) that it would be possible to do something like

o24_naive = op("Id", s[2])*op("Id", s[4]) + op("c†", s[2])*op("c†", s[4]);
o31_naive = op("Id", s[1])*op("Id", s[3]) + op("c†", s[3])*op("c†", s[1]);
ψ_naive = product(o31_naive, product(o24_naive, ψ0));
correlation_matrix(ψ_naive, "c†", "c†")

However, this returns

4×4 Matrix{Float64}:
 0.0           0.0          -2.77556e-16  0.0
 0.0           0.0           0.0          5.55112e-17
 2.77556e-16   0.0           0.0          0.0
 0.0          -5.55112e-17   0.0          0.0

Am I missing something crucial with this naive approach? eg, is it necessary to use enable_auto_fermion()? Or is the best practice to always use OpSum/AutoMPO?

[miles]

Great question. So the situation with fermions is indeed a bit confusing in ITensor right now. It would be simpler (though harder to use) if you had to do everything by hand, but the current setup is that some things are done automatically for you and others aren’t. The handling of Jordan-Wigner string in particular as well as fermionic swaps is only done automatically for you in three places in the code:

• the AutoMPO / OpSum system
• the apply function
• the correlation_matrix function

In any other place you have to account for minus signs and Jordan-Wigner string yourself. So in your code you need two changes:

1. you need to add J-W string to the second terms in your o24_naive and o31_naive operators and also possibly flip some minus signs there
2. you need to use the apply function when applying these operators (may not strictly be necessary for this case but better to be safe)

So some code that almost gets this right is this:

  o24_naive = op("Id", s[2])*op("Id",s[3])*op("Id", s[4]) + op("c†", s[2])*op("F",s[3])*op("c†", s[4]);
  o31_naive = op("Id", s[1])*op("Id",s[2])*op("Id", s[3]) + op("c†", s[3])*op("F",s[2])*op("c†", s[1])

ψ_naive = apply([o31_naive, o24_naive], ψ0)

I say almost because I’m getting some sign flips from the correct correlation matrix, so I think there may be an extra minus sign needed in one or both definitions of o24_naive and o31_naive, when you go through the Jordan-Wigner mapping correctly.

Finally, the auto fermion system isn’t ready to use yet, and is only in the code right now as a developer-level feature for us to keep testing it out. I wish it was ready sooner but there are some conceptual issues I’m still working through!

[nGn]

Thanks for the very clear answer! When inserting the JW strings by hand, is it important then to order the operators by hand according to their index? It looks like this will make your suggested code give the right signs:

o24_naive = op("Id", s[2])*op("Id",s[3])*op("Id", s[4]) + op("c†", s[2])*op("F",s[3])*op("c†", s[4]);
o31_naive = op("Id", s[1])*op("Id",s[2])*op("Id", s[3]) - op("c†", s[1])*op("F",s[2])*op("c†", s[3])

ψ_naive = apply([o31_naive, o24_naive], ψ0)

Oh, and just to clarify: the order in which operators are applied in OpSum is from right to left, while the order used by apply is sequential in the array index, correct? That is,

\verb| product(MPO(OpSum() + ("A", 1, "B", 2), s), ψ)| \quad \longrightarrow \quad A_1B_2|\psi\rangle \\ {} \\ \verb| apply([op("A", s[1]), op("B", s[2])], ψ) | \quad \longrightarrow \quad B_2A_1 |\psi\rangle

It almost goes without saying, but thank you for the beautiful work you and Matt have done with creating/maintaining/supporting ITensor!

[miles]

Yes, you figured out the solution. Great. The reason for the minus sign is indeed related to the ordering of the operators, because the Jordan-Wigner string mapping is defined for the case of the operators ordered left to right (so c^\dagger_1 coming before c^\dagger_3). In that ordering, the operator you are trying to make has a minus sign in front of the c^\dagger_1 c^\dagger_3 term so that same minus sign is neeed in the construction using the op function that you did. A way to understand this is that after the mapping involving the string, the operators now commute, so the minus sign is like a “memory” of what order they were in before they were changed into commuting operators.

Yes your understanding about how things are ordered in OpSum versus apply is definitely correct.

Appreciate the kind words about ITensor! We really like it too and use it a lot ourselves and are working hard to make it better. Please let us know any suggestions you may have.

[hmargasak]

Sorry to reopen this thread if it’s not appropriate, but this isn’t exactly a separate question. I’m curious if the expect function falls under the apply function in your list of functions that don’t require manual JW strings?

[miles]

Well, expect is mainly for a single, local operator like on a single site. So then JW strings wouldn’t play a role i.e. there is no handling of them needed anyway. So a short answer would be “yes” that expect will work automatically for fermions. But were you thinking of a case where that wouldn’t be the case? What operator were you planning to compute the expected value of?

[hmargasak]

Sorry, I should’ve been more clear. I was worried that, since bosonic operators with a JW string will anticommute properly on different sites but not commute properly on the same site (this is at least the way I understand it, I could absolutely be wrong), I thought that single-site operators like Sz might not work properly for Electron systems without modification. Thanks for clearing that up!