The convergence of entropy in spinless p-wave superconductor chain

Hi, everyone
I am trying to perform some DMRG calculations for spinless p-wave Kitaev chain. Specifically, I want to calculate the entanglement entropy of this chain, but I can’t always get the convergent results through changing the the number of sweeps, maximal bond dimensional, cutoff and noise paramerters. A typical code is given as follows

using ITensors

let
N = 40
t=1.0
mu = 1.0
delta = 1.0
f = open(“Kitaev-test.dat”,“w”)
s = siteinds(“Electron”,N)
b=20

for nn = 1200 : 100 : 2000
@show nn

ampo = OpSum()
for j = 1 : N-1
	ampo += -t ,"Cdagup",j,"Cup",j+1
	ampo += -t ,"Cdagup",j+1,"Cup",j
	ampo += delta ,"Cup",j,"Cup",j+1
	ampo += delta ,"Cdagup",j+1,"Cdagup",j			
end

for j = 1 : N
    ampo += -mu ,"Cdagup",j,"Cup",j
end

H = MPO(ampo, s)

sweeps = Sweeps(30)
setmaxdim!(sweeps, 100,200,400,800,1000,nn)
setcutoff!(sweeps, 1E-5,1E-6,1E-7,1E-8,1E-8)
setnoise!(sweeps, 1E-5,1E-6,1E-7,1E-8,1E-10)

#@show sweeps

Initialize wavefunction to be bond

dimension 10 random MPS with number

of particles the same as state

psi0 = randomMPS(s, 40)

Start DMRG calculation:

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

function entropyvonneumann(psi::MPS, b::Int)
    s = siteinds(psi)
    orthogonalize!(psi, b)
    _,S = svd(psi[b], (linkind(psi, b-1), s[b]))
	SvN=0.0
    for n in 1: dim(S, 1)
        p = S[n,n]^2
        SvN -= p * log(p)
    end

return SvN
end

SvN = entropyvonneumann(psi, b)
println("$nn $energy $(SvN)")
write(f,"$nn $energy $(SvN)\n")
	
end
close(f)

end

I want to get the entanglement entropy as a function of maximal bond dimension. But my results always show oscillating behavior just like following picture and can’t obtain a convergent stable value.

I also test other parameters including the number of sweep, noise etc. but the results are similar. Who can help me to fix it or explain why ?

1 Like

Is the system gapless for the Hamiltonian you are using? Just curious. Also I would recommend starting with even shorter systems like N=16 just to get your parameter choices tested.

The main reason I think you might be seeing these inconsistent results is because you are doing a lot of sweeps at the largest bond dimension, whereas it’s usually a better plan to do lots of initial sweeps at smaller bond dimensions then only a few toward the end at larger bond dimensions.

So try e.g. doing like 40 sweeps at bond dimension 10, then do 10 more at bond dimensions 20 and 40, then ramp up to larger bond dimensions like in the hundreds, say, for the last 20 sweeps. Just what I would try but I can’t guarantee that will work. Please give it a try.

The Hamiltonian is a topological superconductor phase with superconductor gap for used parameters, e.g. mu=1.0, t=1.0, delta=1.0. When change the chemical potential mu, it become other phases. Actually for other phases, I also got a similar oscillating behavior for entanglement entropy.
I tryied the ramping method of bond dimensions you suggested, I still cannot obtain the convergent stable entanglement entropy for a system with N=20 sites. A typical setting as follows:

setmaxdim!(sweeps, 10,10,10,10,10,10,10,20,20,20,20,20,20,20,20,
40,40,40,40,40,40,40,40,80,80,80,100,100,100,100,200,200,200,200,
300,300,300,400,400,400,400,nn) #nn maximal bond dimension varible
setcutoff!(sweeps, 1E-4,1E-4,1E-4,1E-4,1E-4,1E-4,1E-4,1E-4,1E-4,1E-4,1E-5,1E-5,
1E-5,1E-5,1E-5,1E-5,1E-6,1E-7,1E-8)

I wonder if the ITensor is suitable for superconducting systems ? Because the Hamiltonian contains some non-conserved terms of particles, e.g. cupcup, cdagupcdagup.

Taking a closer look at your code, I believe I see the issue: your Hamiltonian only involves terms that connect up-spin electrons. There are no terms involving down-spin electrons. While this might possibly work in the case of quantum number conservation, you are also not conserving quantum numbers, so I think this combination of choices could lead to some rather strange behavior.

Indeed when I try to simulate this Hamiltonian, not only is the entropy behaving quasi-randomly but even the particle density is basically acting random in space each run. I think it’s because basically it’s an ill-defined Hamiltonian.

Please revisit the definition of the Hamiltonian you’re trying to simulate. Should there also be terms acting on down-spin electrons (involving operators like “Cdn” and “Cdagdn”)? Or were you intending to simulate spinless fermions? If that’s the case then you can use the “Fermion” site type which does not carry any spin.

Yes. after modifing the code according to your suggestion, it can get convergent results now. Thank you very much, Miles.

I change “Electron” to “Fermion” and “Cdagup” (“Cdag”) to “Cdag” (“C”), then everything is fine.

Glad to hear it!