Fermion site-type operator `c`: `inner()` indicates `c` on different sites commute

[christopherdavidwhit]

Hello,
I wanted to make sure c,cdag on Fermion sites worked in the way I expected, so I made a 3-site MPS, occupied on every site, and did fermion annihilation on sites 1 and 2, in both orders. (Code below.)

I can tell from the indices that apply(c2, ψ) (where c2 is annihilation on site 2 and ψ is my MPS) does something nonlocal. I assume that nonlocal something is applying the Jordan-Wigner string—there’re bits and pieces of documentation that hint that direction, and a quick poke at the code does too. (Maybe I’m reading those hints incorrectly!)

But then I do

order1 = apply(c2, apply(c1, ψ))
order2 = apply(c1, apply(c2, ψ))
  
@show inner( order1, order2)

and find that the inner product is 1, where I’d expect -1! (Again, full code and output below.)

My suspicion is that there’s an off-by-one error, and the J-W string with c_j & c^\dagger_j mistakenly includes site j, where it should only include everything to the right of (but not including) j. Part of why I think this is that when I start with a state

φ = MPS(ss, ["0", "1", "0"])

and do c^\dagger_1 or c^\dagger_3, one or the other should pick up a minus sign, but in fact neither picks up a minus sign—or, rather, in both cases the tensors pick up two minus signs. (I’m fairly confused about which way the JW string is going.) (Again, full code and output below.)

It’s tempting to say this is a bug, but I don’t understand the intended function of the code well enough to say! Maybe my expectations are wrong?

So—what’s the story? Is there something basic I’m missing?

Best,
Christopher

Full code:

  using ITensors
  using ITensorMPS
  using Pkg

  versioninfo()
  Pkg.status("ITensors")
  Pkg.status("ITensorMPS")

  ss = siteinds("Fermion", 3)
  ψ = MPS(ss, ["1", "1", "1"])


  @show ψ

  c1 = op("c", ss[1])
  c2 = op("c", ss[2])

  @show apply(c2, ψ)

  order1 = apply(c2, apply(c1, ψ))
  order2 = apply(c1, apply(c2, ψ))

  @show inner( order1, order2)

  println("=====")

  φ = MPS(ss, ["0", "1", "0"])
  @show φ[1] |> array
  @show φ[2] |> array
  @show φ[3] |> array

  println("-----")

  cdag1φ = apply(op("cdag", ss[1]), φ)
  @show cdag1φ[1] |> array
  @show cdag1φ[2] |> array
  @show cdag1φ[3] |> array
  
  println("-----")

  cdag3φ = apply(op("cdag", ss[3]), φ)
  @show cdag3φ[1] |> array
  @show cdag3φ[2] |> array
  @show cdag3φ[3] |> array

Output:

Julia Version 1.8.2
Commit 36034abf260 (2022-09-29 15:21 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
  Threads: 1 on 8 virtual cores
Status `~/.julia/environments/v1.8/Project.toml`
  [9136182c] ITensors v0.6.16
Status `~/.julia/environments/v1.8/Project.toml`
  [0d1a4710] ITensorMPS v0.2.4
ψ = MPS
[1] ((dim=2|id=578|"Fermion,Site,n=1"), (dim=1|id=379|"Link,l=1"))
[2] ((dim=1|id=379|"Link,l=1"), (dim=2|id=238|"Fermion,Site,n=2"), (dim=1|id=102|"Link,l=2"))
[3] ((dim=1|id=102|"Link,l=2"), (dim=2|id=11|"Fermion,Site,n=3"))

apply(c2, ψ) = MPS
[1] ((dim=2|id=578|"Fermion,Site,n=1"), (dim=1|id=283|"Link,l=1"))
[2] ((dim=2|id=238|"Fermion,Site,n=2"), (dim=1|id=283|"Link,l=1"), (dim=1|id=540|"Link,l=2"))
[3] ((dim=2|id=11|"Fermion,Site,n=3"), (dim=1|id=540|"Link,l=2"))

inner(order1, order2) = 1.0
=====
φ[1] |> array = [1.0; 0.0;;]
φ[2] |> array = [0.0 1.0;;;]
φ[3] |> array = [1.0 0.0]
-----
cdag1φ[1] |> array = [0.0; -1.0;;]
cdag1φ[2] |> array = [0.0; -1.0;;;]
cdag1φ[3] |> array = [1.0; 0.0;;]
-----
cdag3φ[1] |> array = [1.0; 0.0;;]
cdag3φ[2] |> array = [0.0; -1.0;;;]
cdag3φ[3] |> array = [0.0; -1.0;;]

[ryanlevy]

You should assume that odd numbers of fermion operators will need their strings handled manually (OpSum doesn’t support this yet).

Another important thing to note, once you write c1 = op("c", ss[1]), then c1 is an ITensor and has no knowledge of being fermionic as it has no QNs.
The strings are handled at an operator level, so OpSum and some of correlation_matrix can insert JW strings.

See the discussion here for some more: Proper way to construct fermion gates

[christopherdavidwhit]

Thanks for linking to https://itensor.discourse.group/t/proper-way-to-construct-fermion-gates/172: afaict it doesn’t solve my problem, but it is a help!

I’m a little confused—where does OpSum come into it? AFAICT the only ITensor functions I’ve used are MPS, op and apply, and it’s not obvious to me that any of those is going to call OpSum.

Also, you write

once [I] write c1 = op("c", ss[1]), then c1 is an ITensor and has no knowledge of being fermionic as it has no QNs. The strings are handled at an operator level, so OpSum and some of correlation_matrix can insert JW strings.

I’m a little confused by this. It doesn’t seem consistent with Miles’ comment in that other thread that apply automatically handles the JW string. It also doesn’t seem consistent with the observed behavior: in the code I posted, applying c or cdag

• changes the link indices at the site at which I apply it, and
• changes tensors at other sites by multiplying in a sign.
  (You can see both of these in the output in my original question.) Is this really the expected behavior for multiplying in an ITensor that doesn’t know anything about JW strings?

Again I’m probably missing something basic?

[christopherdavidwhit]

If this particular apply does an orthogonalize behind the scenes—IIRC applying an MPO ends up doing something like this—that could explain the observed behavior?

But the documentation (and Miles’ comment in the other thread) do suggest that apply should take care of JW strings for me.

[ryanlevy]

@christopherdavidwhit

I’m a little confused—where does OpSum come into it? AFAICT the only ITensor functions I’ve used are MPS, op and apply, and it’s not obvious to me that any of those is going to call OpSum.

Sorry for the confusion, you’re right. I should have added some context, your code snippets don’t use it of course. If you do decide to use OpSum anywhere though, it wouldn’t support it (and would error saying this).

I did a quick look and cannot see anywhere that apply handles JW strings, particularly odd strings with dense tensors. Sometimes things are passed to OpSum sneakily which may make it seem like strings are being taken care of.

@christopherdavidwhit

changes tensors at other sites by multiplying in a sign.
(You can see both of these in the output in my original question.) Is this really the expected behavior for multiplying in an ITensor that doesn’t know anything about JW strings?

You are also correct that there is an orthogonality shift when applying tensors, so the id of the linkdims can change via that. Because of the gauge degree of freedom, looking at the tensors can sometimes be unhelpful, particularly away from the orthogonality center where everything is a unitary. I expect the orthgonalize call may have introduced these negatives.

φ = MPS(ss, ["0", "1", "0"])
@show φ[1] |> array
@show φ[2] |> array
@show φ[3] |> array

println("-----")

orthoφ = orthogonalize(φ,1) 
@show orthoφ[1] |> array
@show orthoφ[2] |> array
@show orthoφ[3] |> array

output:

φ[1] |> array = [1.0; 0.0;;]
φ[2] |> array = [0.0 1.0;;;]
φ[3] |> array = [1.0 0.0]
-----
orthoφ[1] |> array = [-1.0; -0.0;;]
orthoφ[2] |> array = [0.0; -1.0;;;]
orthoφ[3] |> array = [1.0; 0.0;;]

---

For your first example, here’s how I would write it with a string

c1 = op("c", ss[1])
c2 = op("c", ss[2])
F1 = op("F",ss[1])

order1 = apply([c2*F1,c1],ψ) # note this is inefficient doing c1*F1 !!
order2 = apply([c1,c2*F1],ψ)

@show inner( order1, order2)

Which produces

inner(order1, order2) = -1.0

[christopherdavidwhit]

I guess this is it: thank you!