Unexpected density waves for one dimensional spinless free fermions with periodic boundary condition

ITensor Julia Questions
,
1

Dear ITensor team,
I am a fresh man to using the ITensor library and want to test some physical systems. I make a DMRG calculation for the one dimensional half-filling spinless fermions which have a simplest Hamiltonian,

\hat{H} =J_f\sum_{<i,j>}(c_i^\dag c_j + h.c.).

I employ both the periodic boundary condition (PBC) and the open boundary condition (OBC) for a 32-site chain. When using the PBC, the density profile shows some unexpected density waves. While when applying the OBC, I get the desired results that each site has an averaged particle number 0.5. So, do I missed something? Is this just the boundary effect or something else? I known that the OBC is naturally for DMRG, but when I want to extract the central charge, PBC is better. I really appreciate if someone can help me fix this issue.

Best,

Crliu

The code is as following,

using ITensors

  include(get(ARGS, 1, "input.jl"))

  sweeps = Sweeps(nsweep, sweeps_args)


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

  os = OpSum()
  for b in 1:N-1
    os .+= -Jf, "Cdag", b, "C", b + 1
    os .+= -Jf, "Cdag", b + 1, "C", b
  end

  os .+= -Jf, "Cdag", 1, "C", N
  os .+= -Jf, "Cdag", N, "C",1

  H = MPO(os, sites)

  state = [isodd(n) ? "0" : "1" for n in 1:N]


  # Initialize wavefunction to be bond
  # dimension 10 random MPS with number
  # of particles the same as `state`
  psi0 = MPS(sites, state)

  # Check total number of particles:
  @show flux(psi0)

  # Start DMRG calculation:
  energy, psi = dmrg(H, psi0,sweeps)

  nmea = fill(0.0, N)
for j in 1:N
    orthogonalize!(psi, j)
    psidag_j = dag(prime(psi[j], "Site"))
    nmea[j] = scalar(psidag_j * op(sites, "N", j) * psi[j])
  end

  println("Density:")
  for j in 1:N
    println("$j $(nmea[j])")
  end
  println()

And the input file is,

N = 32
Jf = 1.0
nsweep = 10
sweeps_args = [
  "maxdim" "mindim" "cutoff" "noise"
  50 20 1e-12 1e-6
  100 20 1e-12 1e-7
  200 20 1e-12 1e-8
  800 20 1e-12 0
  800 20 1e-12 0
  800 20 1e-12 0
  800 20 1e-12 0
  1000 20 1e-12 0
  1000 20 1e-12 0
  1000 20 1e-12 0
  1000 20 1e-12 0
  2000 20 1e-12 0
  2000 20 1e-12 0
  2000 20 1e-12 0
  2000 20 1e-12 0
  2000 20 1e-12 0]

The output file with PBC is,

flux(psi0) = QN("Nf",16,-1)
After sweep 1 energy=-19.868252432167086  maxlinkdim=4 maxerr=1.08E-06 time=59.945
After sweep 2 energy=-20.008149046553754  maxlinkdim=16 maxerr=4.15E-08 time=0.688
After sweep 3 energy=-20.302023339629734  maxlinkdim=64 maxerr=6.89E-09 time=1.267
After sweep 4 energy=-20.306338337138502  maxlinkdim=249 maxerr=9.99E-13 time=14.975
After sweep 5 energy=-20.30634077359003  maxlinkdim=471 maxerr=9.99E-13 time=5.642
After sweep 6 energy=-20.306340774982772  maxlinkdim=506 maxerr=9.99E-13 time=7.296
After sweep 7 energy=-20.306340775016608  maxlinkdim=512 maxerr=9.99E-13 time=5.241
After sweep 8 energy=-20.30634077502393  maxlinkdim=513 maxerr=9.99E-13 time=5.692
After sweep 9 energy=-20.30634077502555  maxlinkdim=513 maxerr=1.00E-12 time=6.103
After sweep 10 energy=-20.30634077503053  maxlinkdim=513 maxerr=1.00E-12 time=6.827
Density:
1 0.5152587190593534
2 0.48474128093742286
3 0.5152587190477542
4 0.4847412809299222
5 0.5152587190423987
6 0.484741280930165
7 0.5152587190337196
8 0.4847412809547065
9 0.515258719018605
10 0.4847412809674652
11 0.5152587190184027
12 0.48474128096861385
13 0.5152587190222062
14 0.48474128097003827
15 0.5152587190276828
16 0.4847412809691997
17 0.5152587190297699
18 0.48474128097246433
19 0.5152587190266187
20 0.4847412809774208
21 0.5152587190301435
22 0.484741280978604
23 0.5152587190263962
24 0.48474128098039276
25 0.5152587190355917
26 0.4847412809718471
27 0.5152587190679927
28 0.4847412809648106
29 0.5152587190729093
30 0.4847412809555902
31 0.5152587190666889
32 0.4847412809451359

The output file with OBC is,

flux(psi0) = QN("Nf",16,-1)

After sweep 1 energy=-19.86890676875406 maxlinkdim=4 maxerr=4.17E-07 time=48.451

After sweep 2 energy=-20.01340215325526 maxlinkdim=16 maxerr=4.14E-08 time=0.374

After sweep 3 energy=-20.016383335926896 maxlinkdim=58 maxerr=4.17E-09 time=0.538

After sweep 4 energy=-20.01638789703774 maxlinkdim=76 maxerr=9.97E-13 time=10.015

After sweep 5 energy=-20.016387900367846 maxlinkdim=80 maxerr=9.96E-13 time=0.808

After sweep 6 energy=-20.01638790036378 maxlinkdim=80 maxerr=9.92E-13 time=0.599

After sweep 7 energy=-20.016387900365576 maxlinkdim=80 maxerr=9.95E-13 time=0.501

After sweep 8 energy=-20.01638790036632 maxlinkdim=80 maxerr=9.92E-13 time=0.510

After sweep 9 energy=-20.016387900366198 maxlinkdim=80 maxerr=9.92E-13 time=0.507

After sweep 10 energy=-20.016387900366222 maxlinkdim=80 maxerr=9.92E-13 time=0.565

Density:

1 0.49999999999686257
2 0.4999999999986735
3 0.49999999999702505
4 0.5000000000001048
5 0.4999999999989829
6 0.5000000000007757
7 0.4999999999958785
8 0.5000000000068969
9 0.49999999999552275
10 0.5000000000020647
11 0.4999999999967063
12 0.50000000000123
13 0.5000000000002627
14 0.5000000000039149
15 0.49999999999395034
16 0.5000000000026557
17 0.500000000004377
18 0.49999999999983447
19 0.49999999999801303
20 0.4999999999986424
21 0.49999999999657146
22 0.49999999999759576
23 0.5000000000047478
24 0.499999999992202
25 0.5000000000000361
26 0.49999999999460587
27 0.500000000002678
28 0.5000000000010292
29 0.5000000000047818
30 0.4999999999884952
31 0.5000000000208796
32 0.5000000000039868
3

Please try putting the opposite sign for the hopping term connecting site N to site 1. Does that address the issue?

4

Hi, Miles, thank you for your reply, it works. But I still have two questions: 1. Why this can address the issue? It seems that the Hamiltonian is changed after adding an opposite sign in the N-to-1 hopping term. 2. Should I aways need to add an opposite sign in the N-to-1 hopping term If I work for an arbitrary fermionic system?

5

Glad it worked.

  1. To understand this issue, itā€™s important to understand that the way the ITensor OpSum system currently handles fermions is by mapping them to long-range ā€œbosonicā€ (i.e. commuting) operators with ā€œstringsā€ (Jordan-Wigner string operators) attached to them. So an operator like c_N becomes F_1 F_2 \cdots F_{N-1} b_N. This mapping imposes a certain one-dimensional ordering on the operators.

    But I observe from the above that actually the way the mapping is done in ITensor is very consistent with the notation one would want ā€œon paperā€ whereas the way you coded it originally is the opposite. What I mean is that if you keep continuing the pattern of c^\dagger_b c_{b+1} that you started with, you would get c^\dagger_1 c_2 then c^\dagger_2 c_3 etc. until c^\dagger_{N-1} c_N then finally c^\dagger_N c_1, if you think of identifying (N+1) \leftrightarrow 1.

    So then you could write the code:

      os .+= -Jf, "Cdag", N, "C", 1
  os .+= -Jf, "Cdag", 1, "C",N

    which would continue to have a minus sign like your other terms, but where Iā€™ve flipped the order of the operators to the more consistent ordering above, and it should have the same effect as changing the sign in your original code since the operators anticommute.

  2. Iā€™m not sure what you are asking by asking if you ā€œalwaysā€ need to add an opposite sign. The answer to this is that the operators anticommute so that their ordering is important. If you write them in one ordering you will need a certain sign for it to be correct, and if you write it in the other ordering you will need the opposite sign. I think the real question is whether the pattern followed by the operators through the bulk of the system continues around the periodic boundary or not, and happily I think it does (otherwise it could be a surprise to users).

Please let me know if you have any further questions about this.

6

Hi, Miles, thank you for your patience. For question 1, if I understand correctly, what you mean is that the Fermionic operators anticommute. Numerically, we can pre-define an ordering like: c_1^\dag c_2,c_2^\dag c_3, \cdots,c_{N-1}^\dag c_N. In the case of PBC, an additional term c_N^\dag c_1 arises. This additional term does not obey the pre-defined ordering, and we can correct this term by rewriting it as c_1 c_N^\dag or by introducing an opposite sign like -c_N^\dag c_1. This correction does not change anything of the original Hamiltonian. Am I right? For question 2, you have answered all my confusions.

7

Thanks for pushing to clarify, because I found my answer above was too quick and I missed something important.

The real answer about whether the sign can be the one you wrote or needs to be different turns out to depend on the total number of fermions in the system. The reason for this can be understood by thinking of mapping the operator c^\dagger_1 c_N to a pair of bosonic operators connected by Jordan-Wigner string, and then noticing that the string will measure the total number of fermions in the system.

So then in some experiments Iā€™ve done, Iā€™ve confirmed that if the total number of sites is even, but is two times an odd number, then I get the flat density you are seeking. But if the number of sites is even and is two times an even number, then I get the staggered density.

I think this is the expected behavior, because there can be different types of boundary conditions with fermions which can be viewed as the presence or absence of a ā€˜fluxā€™ threaded through the middle of the loop.

Really the most important thing to compare to would be a free-fermion calculation done just by diagonalizing the single-particle Hamiltonian and reasoning about filling the N lowest energy states and computing the total density that way. What weā€™d want is for ITensorā€™s OpSum system to lead to the same results that one would get that way. So please let us know if thatā€™s not the case.

8

Hi, Miles
Thanks for the clarification. Itā€™s really helpful for me.
Best,
Crliu

9

Hi Miles,

The thing you mentioned about the total number of sites being even and twice an odd number to get a flat density, do you think the same applies to spinfull fermions in 2D lattices?
Because Iā€™ve got similar issues with the triangular lattice

Thanks,
Abraham

10

Hi Abraham,
Yes I would assume that it would apply in some form. This topic is a rather tricky and technical one, and is really about how to properly define global symmetry ā€˜sectorsā€™ of fermionic systems with periodic boundary conditions, more than it is about tensor networks as I think you might agree. I have been seeing theoretical talks recently where people are currently working on how to best think about these matters (e.g. they come up a lot when studying dualities, like mapping spinless fermions to Ising spins and back using Jordan-Wigner). The current thinking about it, as I understand, is that one really has to study multiple types of boundary conditions and multiple filling sectors (i.e. even and odd numbers of bulk fermions) to fully understand an entire model. That can seem weird and counterintuitive because with different boundary conditions it seems like they are separate models (different Hamiltonians). But itā€™s the right thing to do for various reasons. One reason is that in the dual language (e.g. spin language) the seemingly separate models all become the same model and you end up with a ground state degeneracy. Another more physical reason is that for fermions the twisted boundary conditions can be viewed as a magnetic flux threaded through the ring (in the periodic 1D example).

Again, I think your best bet might be to do some studies of the non-interacting limit where you can quickly try various boundary conditions and fillings, and also get an understanding of the structure of the states. (E.g. certain boundary conditions lead to the single-particle states being labeled by half-integer momenta rather than integer momenta.)

Best regards,
Miles

2 Likes