ferromagnetism in 1D chain

I am trying to reproduce the result for ferromagnetism in 1d chain using DMRG particularly from this paper DMRG study of ferromagnetism in a one-dimensional Hubbard model | SpringerLink. I am not able to figure out the problem with my code. It is giving me spin_up number same as spin_down instead of giving me a ferromagnetic ground state

‘’’
using ITensors
using LinearAlgebra
using Plots
a=Float64
b=Float64
let
L=20
U=5
h=0
N = 40
sites = siteinds(“Electron”, N;conserve_nf=true)
@show sites[1]
t =-1
tt=0
mu=0
ampo = OpSum()
for j=1:N-1
ampo += t,“Cdagup”,j,“Cup”,j+1
ampo += t,“Cdagup”,j+1,“Cup”,j
ampo += t,“Cdagdn”,j,“Cdn”,j+1
ampo += t,“Cdagdn”,j+1,“Cdn”,j
ampo += U,“Nup”,j,“Ndn”,j
ampo += mu,“Ntot”,j
ampo += h,“Sz”,j
end
for j=1:N-2
ampo += tt,“Cdagup”,j,“Cup”,j+2
ampo += tt,“Cdagup”,j+2,“Cup”,j
ampo += tt,“Cdagdn”,j,“Cdn”,j+2
ampo += tt,“Cdagdn”,j+2,“Cdn”,j
end
ampo += h,“Sz”,N
ampo += mu,“Ntot”,N
ampo += U,“Nup”,N,“Ndn”,N
state1 = [isodd(n) ? “Emp” : “Up” for n=1:N]
H = MPO(ampo, sites);
sweeps = Sweeps(100);
maxdim!(sweeps, 64,64,100,100,128,128,160,160,200,200,256,256,300,300,400,400,400);
setcutoff!(sweeps,1E-15);
setnoise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
psiii = randomMPS(sites,state1,8);
energy, psi = dmrg(H, psiii, sweeps);
expNup=expect(psi,“Nup”);
expNdn=expect(psi,“Ndn”);
nt=expect(psi,“Ntot”)
println(“$nt”)
Szexp=expect(psi,“Sz”)
nbar=sum(expNup)/N;
ntot=sum(expNdn)/N;
println(“$nbar”)
println(“$ntot”)
println(“$Szexp”)
zz=correlation_matrix(psi,“Sz”,“Sz”);
z1=diag(zz)
sss=size(zz)
println(“$z1”)
‘’’
could you please help me figure out the mistake?

We’ve just discussed a very similar situation in another thread, here’s what I think is going on.

tl; dr: you may be finding the exact ground state, which never breaks symmetry in a finite system.

Vladimir’s guess is probably the right one, I agree.

Symmetry breaking only happens when there is an infinitesimally small magnetic field “picking out” which way the symmetry actually breaks.

To be more precise, a symmetry broken ferromagnetic state is defined as the result of two limits:

  1. first one takes the system size to infinity while having a non-zero magnetic field
  2. then one takes the limit of the field going to zero within the infinite system

A nice technique to observe symmetry breaking without putting on a bulk magnetic field is to “pin” the open edges of the system. You can read more about the pinning technique in this article.