[Julia] Spin gap behavior and boundary effects

Question

Hello ITensors Team ,

Ok I decide ask a new question and close the last one.
i'm trying to investigate the spin Gap behavior of Heisenberg spin-1/2 J1-J2(J1 and J2 >0, non-frustrated) Isotropic model with interations between first pair of spins with exchange paramenter J1 and the following pair of spins with interations mediated by J2 an so on. This simple model presents an topological phase transitions at J1=J2 between "Haldane-Even" and "Haldane-Odd" phases.
I wrote an simple code for this model considering open boundary conditions as follows :

using ITensors
using DelimitedFiles
int = range(0.95, step=0.5, stop=1.5)
Data = []
Data_TW = []
for i in int   
    let
        N = 120
        sites = siteinds("S=1/2", N; conserve_qns=true)
        J_1 = 1.0 
        J_2 = i
        ampo = AutoMPO()

        for j=1:2:N
            ampo += 0.5*J_1,"S+",j,"S-",j+1
            ampo += 0.5*J_1,"S-",j,"S+",j+1
            ampo += J_1,"Sz",j,"Sz",j+1
        end

        for j=2:2:N-2
            ampo += 0.5*J_2,"S+",j,"S-",j+1
            ampo += 0.5*J_2,"S-",j,"S+",j+1
            ampo += J_2,"Sz",j,"Sz",j+1
        end

        H = MPO(ampo,sites)
        state1 = []
        push!(state1, "Up", "Up")
        for n=2:N/2
            push!(state1, "Up", "Dn")
        end
        init_state = [isodd(n) ? 1 : 2 for n = 1:N]
        psi0 = randomMPS(sites,init_state,32)
        psi1 = randomMPS(sites,state1,32)

        @show flux(psi0)
        @show flux(psi1)
        sweeps0 = Sweeps(10)
        maxdim!(sweeps0, 160,320,600,800,1000)
        cutoff!(sweeps0, 1e-11)
        noise!(sweeps0, 1e-8,1e-9, 1e-10,0.0)
        energy0,psi0 = dmrg(H,psi0,sweeps0)
        H02 = inner(H,psi0,H,psi0)
        E0 = inner(psi0,H,psi0)
        var0 = H02 - E0^2
        @show var0
        sweeps1 = Sweeps(18)
        maxdim!(sweeps1, 320,600,800,1000,1200)
        cutoff!(sweeps1, 1e-11)
        noise!(sweeps1, 1e-8, 1e-9,1e-10,0.0)
        energy1,psi1 = dmrg(H,psi1,sweeps1)
        H12 = inner(H,psi1,H,psi1)
        E1 = inner(psi1,H,psi1)
        var1 = H12 - E1^2
        @show var1
        gap = energy1-energy0
        sgap = N*gap
        println("The Gap is =",gap)
        println("The Scaled gap is=", sgap)
        println("The J_2  value is :", J_2)
        push!(Data, ["$J_2" "$energy0" "$gap" "$sgap"])
    end
    open("Data(test).dat", "w") do f
        writedlm(f, Data, '\t')
    end
end

In my code the first loop concerning the J1 interacting spin terms i sum up to N for each 2 points and the second loop for the J2 interacting terms i sum up to N-2 terms(OBC).
The strange fact is, i cannot reproduce the literature results because my spin gap behaves like an Gapped-Gapless transition. The curious fact is, when i change the first loop summation(Up to N-2), i reproduce the correct behavior with an clear quantum phase transition at J = J2/J1 = 1.0 , but with the penalty of the missing interaction at last pair of spins(which shows up when i compute the local magnetization).
I already tried several initial different states , productMps vs randomMPS wavefunctions and different values of noise term .

Probably i'm missing something, i'm grateful with any help ! Thanks !

Answer 1

Hi, thanks for the question. The j=1:2:N loop makes me a bit nervous since j+1 might be greater than N. But I guess not as long as N is an even number, since the last value there will be j=N-1.

But the key issue here, I would bet, is that you are just not doing enough DMRG sweeps. Especially when obtaining excited states with DMRG, it’s important to check that enough sweeps have been done, because for excited states convergence can be rather slow.

Relatedly, I don’t recommend putting a sequence of DMRG calculations into one big loop (the loop over i here). It’s tempting, but almost every time I’ve seen users do it, the resulting problem is that they are studying different points in their phase diagram using the same DMRG settings, whereas near phase transitions or at various other harder points, they need to do a more careful study of DMRG convergence in both number of sweeps and other parameters such as maxdim and cutoff.

Please also check that at various points your results are converged in those parameters (maxdim and cutoff) too, since if you don’t have convergence in those, it could also strongly affect what gaps you find.

Finally, these kinds of phases (Haldane) phase can have subtle zero-energy edge states which can complicate how one calculates and defines the gap. Please check that the state you are treating as the first excited state is not actually a second, third, or even fourth degenerate ground state.

Hope that helps!

Miles

Comments

Hello miles, Thank you for the reply

Concerning the convergence, i tried to chose an number of parameters(Sweeps, maxdim and cutoff) with returns an variation V  = <H^2> - <H>^2 proportional to E-8 in mostly points(i'm using maxdim up to 1200, cutoff e-12 for small chains). So its possible the system still stuck in an local minima even using the noise term and sweeps numbers sufficiently large to given an variation V proportional to e-8?
Concerning the degenerated ground state, indeed it seems me code target wrongly the ground state of the system after the Quantum phase transition, where the system presents an degenerated non-trivial topological phase with virtual spins at edges. So, just giving as i did an state1 array with total number QN(sz = 1) could not represent the first excited state ? Maybe i had to search excited states at the same sector QN(sz=0) ? can you clarify this to me please ? Thank you !

---

Yes, good question. So if the physics after the phase transition is similar to that of the Haldane phase (the one that is in the same phase as the S=1 Heisenberg chain or the AKLT model), then with open boundary conditions there are *four* degenerate ground states. You can think of this as if the edges each had an S=1/2 edge state and then the four states are the four combinations of two spin 1/2's (singlet, and 3 triplets).

So if you only obtained the lowest-energy Sz=0 and Sz=+1 states, those would have nearly the same energy (up to exponentially small corrections in the system size). To obtain the bulk gap, you'll either need to obtain the lowest 5 eigenstates or else, as you mentioned **obtain the lowest excited state in the Sz=0 sector**.

So yes, I'd recommend looking just at Sz=0, but remember there is a second ground state (Sz=0 triplet) there, so you will need the second excited state.

---

Another idea would be to obtain the ground state and first excited state in the Sz=+1 (or +2 in "ITensor units") sector.

---

Great , thank you Miles for the clear explanation !