Kondo problem of spin-electron mixed chain

Dear ITensors team,

I really appreciate your help for solving my previous problem of finding the ground state of a system with c-electrons and f-electrons. It helped me a lot in my research and now I am trying to study the system in the kondo limit. To save configuration space, I replaced the f-electron with spins.

My system is now an AFM Heisenberg ring (N sites) with external magnetic field and a noninteracting Hubbard ring (N sites). They have inter-ring on-site magnetic coupling:

H=\overbrace{\sum_i^N(\mathbf{S}_i\cdot\mathbf{S}_{i+1}+B_iS_i^z)}^{H_S}\overbrace{-t\sum_{i\sigma}^N(c^\dagger_{i,\sigma}c_{i+1,\sigma}+h.c.)}^{H_c}+\overbrace{g\sum_i^N\mathbf{S}_i\cdot \boldsymbol{\xi}_{i}}^{H_I}

where \boldsymbol{\xi}_{i} is the electron spin.
Following the discussion in the forum, I mapped the system to a mixed chain with spin-\frac{1}{2} on odd number sites and electron on even number sites. To avoid confusion, I will use i as index of the real ring, and n as the index of the mixed chain.

N2 = N*2
sites = Vector{Index}(undef,N2)
for n=1:N2
    if !isodd(n)
        sites[n] = siteind("Electron"; addtags="n=$n",conserve_qns=true,conserve_sz=false)
    else
        sites[n] = Index([QN() => 2],"S=1/2,Site,n=$n")
    end
end

I use N=8 to compare with Lanczos exact diagonalization. The electron number is conserved and set to be the same as the site number N_e=N. My aim is to reach N\sim20 with ITensors DMRG, Julia version.

What I do:
I add B field only for half of the spin sites, so B_{i,odd}=1,B_{i,even}=0.

  1. Starting with t=0.1,g=-1, I managed to find the 3 lowest eigenstates with
nsweeps = 50
maxdim = [32, 50, 50, 100, 200]
noise = [1E-4, 1E-5, 1E-6, 1E-7, 1E-8, 1E-10, 0.0]
cutoff = [1E-10]
obs = DMRGObserver(; energy_tol=1e-8)

state = ["Up", "Up", "Dn", "Dn", "Up", "Up", "Dn", "Dn", "Up", "Up", "Dn", "Dn", "Up", "Up", "Dn", "Dn"]
psi0_init = randomMPS(sites, state, 32)
energy, psi0 = dmrg([HS,HC,HI], psi0_init; nsweeps, maxdim, cutoff, noise, observer=obs, eigsolve_krylovdim=10)
    
psi1_init = randomMPS(sites, state, 32)
energy1,psi1 = dmrg(HS+HC+HI,[psi0],psi1_init; nsweeps, maxdim, cutoff, noise, observer=obs, eigsolve_krylovdim=10)

psi2_init = randomMPS(sites, state, 32)
energy2,psi2 = dmrg(HS+HC+HI,[psi0,psi1],psi2_init; nsweeps, maxdim, cutoff, noise, observer=obs, eigsolve_krylovdim=10)

the outputs are

inner(psi1, psi0) = -4.016258051996359e-9 - 1.04126681321512e-11im
inner(psi2, psi0) = -3.142391111259102e-9 - 2.8309474111014046e-11im
inner(psi2, psi1) = 2.062048445204612e-6 + 5.429848970465785e-12im

Ground State Energy = -8.102545324713711
1st Ex State Energy = -7.581251489870086
2nd Ex State Energy = -7.514436934018416

And from ED, I got

E_0=-8.1037312704267279,\\ E_1=-7.5812903736466550,\\ E_2=-7.5149187832551849.

I didn’t use large maxdim for DMRG but the eigenenergies from DMRG are still close to ED values (\delta E_{j=0,1,2}\sim10^{-3}). I also used the same parameters for a smaller N=4 system, and the eigenenergies from dmrg are very close to ED (\delta E_{j=0,1,2}\sim10^{-8}).

  1. Then I used t=0.5,g=-0.1. With the same maxdim, there was no convergence within 50 sweeps. It also happened for some runs that after 50 sweeps, E_1 given is lower than E_0. For those ‘wrong order’ runs, inner(psi1, psi0) is \sim0.01. After changing maxdim to
maxdim = [32, 50, 50, 100, 200, 400, 800]

I got the E_i :s in right order:

inner(psi1, psi0) = 1.993598762675486e-8 - 1.3074162636901027e-9im
inner(psi2, psi0) = 3.221085119119989e-7 - 9.54181523919444e-11im
inner(psi2, psi1) = -0.00014749734504719873 - 1.504935873834117e-9im

Ground State Energy = -9.75847701396777
1st Ex State Energy = -9.71823684874759
2nd Ex State Energy = -9.717512477902524

While the ED results are

E_0=-9.7607916289105034,\\ E_1=-9.7191608132422580,\\ E_2=-9.7182895952743280.

So with larger maxdim, I got same order of accuracy as in 1: \delta E\sim10^{-3}.

My questions:

  1. Increasing maximum link dimension seems to be necessary for large |t/g| and for large N. However, each sweep becomes slow: with maxdim=800 and N=20, a sweep is like
After sweep 50 energy=-24.856623993624577  maxlinkdim=800 maxerr=9.25E-07 time=937.064

I have tried starting with small t and increasing it in small steps using the output MPS as input, but that didn’t change the things much. For such spin-electron mixed chain, how can I get faster convergence?

  1. When dmrg function was not terminated by the observer within 50 sweeps, it happened that E_0>E_1, where the last 2 sweeps in finding E_0 have energy difference \sim10^{-6}. If I need to increase maxdim from 200 to 800 when I increase N from 4 to 8, in order to avoid E_0>E_1, should I also use larger maxdim when I use N=20 ?

  2. With N=8, there is a \delta E\sim10^{-3} difference between ED and DMRG eigenenergies. Is there anyway to improve that (it can be that my ED code is wrong, but such \delta E is very small for N=4 )?

  3. Remembering what Miles said for my previous question:

would it help if I add next-nearest-neighbor spin coupling for the current system? Or would it help if I go to higher dimension?

Many thanks and regards,
Zhen