Fermionic statistics between different flavors

Question

Hello Miles,

After fixing the single-flavor fermion chain following your advice at http://itensor.org/support/3237/strange-inconsistency-in-free-fermion-benchmark, I went on to add in a second flavor of free fermion. I would expect the GS energy would simply double if the number of fermions is the same for the two flavors. However, this is not what I got. Is this because of the fermionic statistics between the two flavors? For single-flavor, I know that it has been taken care of by "hasfermionstring", but for two flavors I'm not sure whether it has been included properly or not. If not, I would like to know how to incorporate the statistics.

Here is a minimal version of the code for your reference. Thanks a lot. - Mason

using ITensors

function ITensors.space(
  ::SiteType"CFermion";
  conserve_qns=false,
  conserve_Q1=conserve_qns,
  conserve_nfparity=conserve_qns,
  qnname_Q1="Q1",
  qnname_nfparity="NfParity",
 ) 
  if conserve_nfparity && conserve_Q1
    return [
      QN((qnname_nfparity, 0, -2), (qnname_Q1, 0)) => 1,
      QN((qnname_nfparity, 1, -2), (qnname_Q1, +1)) => 1,
      QN((qnname_nfparity, 1, -2), (qnname_Q1, +2)) => 1,
      QN((qnname_nfparity, 0, -2), (qnname_Q1, +3)) => 1
    ]
  end
  return 4
end

ITensors.state(::SiteType"CFermion",::StateName"Emp" ) = 1
ITensors.state(::SiteType"CFermion",::StateName"f1" ) = 2
ITensors.state(::SiteType"CFermion",::StateName"f2" ) = 3
ITensors.state(::SiteType"CFermion",::StateName"f1f2" ) = 4

function ITensors.op!(Op::ITensor, ::OpName"C1", ::SiteType"CFermion", s::Index)
  Op[s' => 1, s => 2] = 1.0
  return Op[s' => 3, s => 4] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"C1d", ::SiteType"CFermion", s::Index)
  Op[s' => 2, s => 1] = 1.0
  return Op[s' => 4, s => 3] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"C2", ::SiteType"CFermion", s::Index)
  Op[s' => 1, s => 3] = 1.0
  return Op[s' => 2, s => 4] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"C2d", ::SiteType"CFermion", s::Index)
  Op[s' => 3, s => 1] = 1.0
  return Op[s' => 4, s => 2] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"F", ::SiteType"CFermion", s::Index)
  Op[s'=>1,s=>1] = +1.0
  Op[s'=>2,s=>2] = -1.0
  Op[s'=>3,s=>3] = -1.0
  Op[s'=>4,s=>4] = +1.0
end

ITensors.has_fermion_string(::OpName"C1", ::SiteType"CFermion") = true
ITensors.has_fermion_string(::OpName"C1d", ::SiteType"CFermion") = true
ITensors.has_fermion_string(::OpName"C2", ::SiteType"CFermion") = true
ITensors.has_fermion_string(::OpName"C2d", ::SiteType"CFermion") = true

let
N = 4   
sweeps = Sweeps(10)
maxdim!(sweeps,10,20,50)
cutoff!(sweeps,1e-10)

sites = siteinds("CFermion",N; conserve_qns = true)

states = ["f1f2","Emp","Emp","Emp"]

ampo = AutoMPO()

for i = 1:(N-1)
    ampo += -1, "C1d", i, "C1", i+1
    ampo += -1, "C2d", i, "C2", i+1
    ampo += -1, "C1d", i+1, "C1", i
    ampo += -1, "C2d", i+1, "C2", i
end

# PBC
ampo += -1, "C1d", N, "C1", 1
ampo += -1, "C2d", N, "C2", 1
ampo += -1, "C1d", 1, "C1", N
ampo += -1, "C2d", 1, "C2", N

H = MPO(ampo,sites)
psi0 = productMPS(sites,states)

energy0,psi = dmrg(H,psi0,sweeps)

return
end

Comments

I'm not a Julia user, but as far as I know on the C++ side, you need to define the C operators with the correct phase, while the A operators are left without additional phases. I *think* the autoMPO is supposed to take care of this for building Hamiltonians, but I'm not sure the extent to which this feature has been fleshed out, but I'm not sure if this is extended to customized sitesets as I haven't read the autompo documentation thoroughly. https://www.itensor.org/docs.cgi?page=tutorials/fermions&vers=julia

In my custom sitesets for Fermions (both two-species and four-species sitesets), I've included those phases and been able to reproduce results from exact diagonalization and other numerical works. I use the A operators and insert Fermionic phases by hand using the F operator if measuring, say, a single particle Green function. Verifying against known cases to ensure that the base code is correct whenever I add a new feature to my test bed code.

The syntax for Julia is different enough that I'm not going to learn enough (right now) to double-check the include and Op[] portions; for instance, I don't know if you need an ITensors.QN or if just QN will suffice, I use the include statements at the top for the "lazy" method Miles mentioned in the other question. Hopefully one of the Julia users can check that part over.

Best wishes,
Jared

---

Jared, thanks for your comment. Following the documentation at the link you provided, I derived the anti-commutation between the two flavors based on the parity operator defined in my code. So, in principle, the fermionic statistics between flavors should already be included by AutoMPO. The "ITensors." issue has been fixed in their updated version.

---

Hi, Mason,

Right, so if it's not that, check your syntax over and over, try variations, etc. and save plenty of backup files.

Something that might be the cause - your qnname_Q1, this is the Fermionic number, right? As far as I know, you still need the -1 to instruct ITensor that it's a Fermionic symmetry. The 1 says that it's a Z-symmetry, while the minus sign indicates that it's Fermionic. This may or may not be handled automatically on the backend, I don't know.  But, for my C++ code, I would have something analogous to
QN((qnname_nfparity, 0, -2), (qnname_Q1, 0,-1)) => 1
instead of your
QN((qnname_nfparity, 0, -2), (qnname_Q1, 0)) => 1.

Again ,the Julia syntax and all is different from the C++, so it's difficult for me to debug in a foreign language. These are all just half-educated guesses. You can attempt to measure things like particle number etc to ensure that the quantities you want conserved are conserved.

Best,
Jared

---

Jared, thanks again, I think I figured it out, please refer to my answer.

Answer 1

I should have realized that ITensor already has a perfect example for two flavors of fermions: the "electron" site type at https://github.com/ITensor/ITensors.jl/blob/main/src/physics/site_types/electron.jl where the spin up and down are the two flavors.

By carefully comparing my code with the "electron" site type, we can see two issues:
1) the parity operator for each flavor is missing (i.e. "F1" and "F2"), even though I ahve the overall parity operator "F".
2) when defining the second fermion operator "C2" and "C2d", we should have Op[s' => 2, s => 4] = -1.0 and Op[s' => 4, s => 2] = -1.0, where the negative signs are crucial for the fermionic statistics.

After fixing the above two issues, everything will be consistent.