Mixed-site (Boson,spin) model Expected value of <ni>

[Ashinxinl]

Hi, itensor team
Thanks for your support so far! Thank you for your patient and detailed response last time. I learned a lot about MPS and MPO techniques. However, I encountered some issues while trying to reproduce the results shown in Figure 1 (a) and (b) of the paper.(Phys. Rev. Lett. 121, 090402 (2018) - Strongly Correlated Bosons on a Dynamical Lattice) No matter how I calculate, I find that the expectation value of sigmaz is always negative when the bosonic density is at ( 2/3 ) or ( 0.88 ). This outcome seems to prevent me from replicating the results presented in the paper. I suspect that this might be due to the DMRG algorithm optimizing for the lowest energy state when spins and bosons are placed in different positions.(Mixed-site (Boson,spin) model)
Below are my results along with a brief snippet of my code for your reference:

using ITensors
let
N = 119
U=10
t=1
phonon_dim=2
alpha=0.5
delta=0.85
beta=0.02
mu=0.53 #ρ=2/3


maxdim = [200, 400]
mindim = [200]
nsweeps=35
cutoff = [1E-12]
noise=[1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,0,0,0,0,0,0,0]
sites = siteinds(n->isodd(n) ? "Boson" : "S=1/2",N)
psi0 = randomMPS(sites)
terms = OpSum()
for i=1:2:(N-2)
terms += -alpha, "Adag",i,"Z",i+1,"A",i+2
terms += -alpha, "A",i,"Z",i+1,"Adag",i+2
end
for i=1:2:(N-2) 
    terms += -t, "Adag",i,"A",i+2
    terms += -t, "A",i,"Adag",i+2
end
for i=1:2:N
  terms += (U/2), "N",i,"N",i
  terms += - (U/2), "N",i
end
for i=1:2:N
  terms += -(mu), "N",i
end
for i=1:2:(N-2) 
  terms += (delta/2), "Z",i+1
end
for i=1:2:(N-2) 
  terms += beta, "X",i+1
end
H = MPO(terms, sites)
energy,psi = dmrg(H,psi0;nsweeps,maxdim,noise,observer=obs,cutoff)
densities = [] 
sigmazs=[]
for i in 2:2:N
     opz = expect(psi, "Z"; sites=i)
     push!(sigmazs, opz)  
end
for i in 1:2:N
    density = expect(psi, "N"; sites=i)
    push!(densities, density)  
end
@show sigmazs
end

I would greatly appreciate any insights or suggestions you might have regarding this issue. Thank you once again for your help!
Best regards,

[miles]

It sounds like maybe your DMRG calculation is not converging very well? Here are some thoughts.

• I notice that while you are using the noise term, you are using a very small (nearly zero) value for the noise. Please try turning it up to a larger value, as large as 1E-3 even to see if that helps.
• It’s good that you are doing a lot of sweeps
• You should probably start with a lower maxdim on your initial sweeps, then only raise it to 200 or 400 toward the end. So something like maxdim = [10,10,10,10, 20,20,40,40,80,80, .... If you don’t want to manually repeat the values you can write maxdim = [repeat([10],4)..., repeat([20],2)..., .... The repeat function is a built-in Julia function
• I would not set the mindim value so high. I would recommend actually just not defining it.

Please try those things and see if they might help

[Ashinxinl]

HI miles
Thank you for your timely and patient response. I followed your suggestions and adjusted parameters such as maxdim and noise , but after many attempts, I found that there were no significant changes in the results. When the boson density is 2/3, the expectation values of each σz site are still all negative, while the article mentions that there are both positive and negative values. Could this be related to the way the MPO is constructed? Since my DMRG calculations are converging, what do you think might be the issue? Here are the results I obtained after changing the parameters. (Zoperation is expectation values of each σz site)Thank you for your help!

Here is the image from the article.FIG1(a)(b)(Phys. Rev. Lett. 121, 090402 (2018) - Strongly Correlated Bosons on a Dynamical Lattice)

[miles]

As to why your results aren’t matching the ones in the article, does the calculation in the article use “hard core” bosons restricted to have a maximum occupancy of 1, or regular bosons which can (in principle) have any number of bosons on a site?

I ask because your code looks like it is just using our “Boson” site type without any “dim” parameter, in which case the maximum occupancy defaults to just 2. If you are trying to approach the limit of bosons with higher occupancies, it’s important to raise the maximum dimension of the boson sites until you can confirm that your results are converged with this dimension value.

You can make more general sites this way:

boson_dim = 4
sites = [isodd(j) ? siteind("Boson";dim=boson_dim) : siteind("S=1/2") for j=1:N]

If you do @show sites you will see that now all of your Boson sites will have a dimension of 4, and generally whatever maximum dimension you set.

[Ashinxinl]

Thank you, Miles! I understand your point now. I am currently studying physics, and when I set the boson_dim to 3, I successfully replicated the results shown in the image. I really appreciate your detailed assistance; it has boosted my confidence and interest in learning physics. I wish you all the best in your work!