# Unexpected density waves for one dimensional spinless free fermions with periodic boundary condition **URL:** \ **Category:** ITensor Julia Questions\ **Tags:** julia, dmrg\ **Created:** [December 25, 2023, 9:39am UTC](https://itensor.discourse.group/t/unexpected-density-waves-for-one-dimensional-spinless-free-fermions-with-periodic-boundary-condition/1377 "2023-12-25T09:39:48Z")\ **Posts on this page:** 9\ **Page:** 1
**Author:** ![crliu](https://avatars.discourse-cdn.com/v4/letter/c/b3f665/32.png) [@crliu](https://itensor.discourse.group/u/crliu)\ **Post date:** [December 25, 2023, 9:39am UTC](https://itensor.discourse.group/t/unexpected-density-waves-for-one-dimensional-spinless-free-fermions-with-periodic-boundary-condition/1377/1 "2023-12-25T09:39:48Z")
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\_{\}(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, ```julia 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, ```Julia 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, ```Julia 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, ```Julia 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 ``` ---
**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\ **Post date:** [December 27, 2023, 5:30am UTC](https://itensor.discourse.group/t/unexpected-density-waves-for-one-dimensional-spinless-free-fermions-with-periodic-boundary-condition/1377/3 "2023-12-27T05:30:55Z")
Please try putting the opposite sign for the hopping term connecting site N to site 1. Does that address the issue? ---
**Author:** ![crliu](https://avatars.discourse-cdn.com/v4/letter/c/b3f665/32.png) [@crliu](https://itensor.discourse.group/u/crliu)\ **Post date:** [December 27, 2023, 6:05am UTC](https://itensor.discourse.group/t/unexpected-density-waves-for-one-dimensional-spinless-free-fermions-with-periodic-boundary-condition/1377/4 "2023-12-27T06:05:00Z")
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? ---
**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\ **Post date:** [December 27, 2023, 8:29pm UTC](https://itensor.discourse.group/t/unexpected-density-waves-for-one-dimensional-spinless-free-fermions-with-periodic-boundary-condition/1377/5 "2023-12-27T20:29:29Z")
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. 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. ---
**Author:** ![crliu](https://avatars.discourse-cdn.com/v4/letter/c/b3f665/32.png) [@crliu](https://itensor.discourse.group/u/crliu)\ **Post date:** [December 28, 2023, 9:25am UTC](https://itensor.discourse.group/t/unexpected-density-waves-for-one-dimensional-spinless-free-fermions-with-periodic-boundary-condition/1377/6 "2023-12-28T09:25:59Z")
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. ---
**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\ **Post date:** [January 4, 2024, 3:26am UTC](https://itensor.discourse.group/t/unexpected-density-waves-for-one-dimensional-spinless-free-fermions-with-periodic-boundary-condition/1377/7 "2024-01-04T03:26:17Z")
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. ---
**Author:** ![crliu](https://avatars.discourse-cdn.com/v4/letter/c/b3f665/32.png) [@crliu](https://itensor.discourse.group/u/crliu)\ **Post date:** [January 4, 2024, 4:46am UTC](https://itensor.discourse.group/t/unexpected-density-waves-for-one-dimensional-spinless-free-fermions-with-periodic-boundary-condition/1377/8 "2024-01-04T04:46:14Z")
Hi, Miles Thanks for the clarification. It’s really helpful for me. Best, Crliu ---
**Author:** ![Abraham](https://avatars.discourse-cdn.com/v4/letter/a/6a8cbe/32.png) [@Abraham](https://itensor.discourse.group/u/Abraham)\ **Post date:** [January 5, 2024, 5:50am UTC](https://itensor.discourse.group/t/unexpected-density-waves-for-one-dimensional-spinless-free-fermions-with-periodic-boundary-condition/1377/9 "2024-01-05T05:50:01Z")
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 ---
**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\ **Post date:** [February 3, 2024, 8:31pm UTC](https://itensor.discourse.group/t/unexpected-density-waves-for-one-dimensional-spinless-free-fermions-with-periodic-boundary-condition/1377/10 "2024-02-03T20:31:40Z")
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