How to calculate the superfluid order parameter (<a>) for Bose-Hubbard model?

DMRG and Numerical Methods
1

Hi, I’m studying the 1D Bose-Hubbard model and trying to reproduce the Mott insulator-superfluid phase transition. Here’s my Hamiltonian:

\hat H=\sum_j [J(\hat a^\dagger_{j+1}\hat a_j+\hat a^\dagger_{j-1}\hat a_j)+\frac{U}{2}\hat n_j(\hat n_j-1)-\mu \hat n_j].

I would like to calculate the expectation value of the annihilation operator as the superfluid order parameter, i.e. \langle\hat a_j\rangle. As far as I know, DMRG doesn’t spontaneously break U(1) symmetry so I would need to add a pinning field for a few sweeps to force symmetry breaking. According to some reference I’ve found (arxiv:2409.15424) this can be done by adding a symmetry breaking term to the edge of the Hamiltonian, i.e.

\hat H_\mathrm{SB}=\sum_j [J(\hat a^\dagger_{j+1}\hat a_j+\hat a^\dagger_{j-1}\hat a_j)+\frac{U}{2}\hat n_j(\hat n_j-1)-\mu \hat n_j]+\varepsilon(a^\dagger_1+a_1),

do a few sweeps with H_\mathrm{SB}, then use the calculated MPS to do DMRG with the original Hamiltonian. Here’s my result:

1272×816 32.5 KB

I would expect the order parameter to be zero for \mu<\mu_{\mathrm c} (Mott insulator), and smooth for \mu>\mu_{\mathrm c} (superfluid). However, here the order parameter is not strictly zero for the Mott insulator phase, and shows fluctuation for the superfluid phase. The phase transition point seems accurate (according to PhysRevB.58.R14741). Is this behavior expected? Where did I do wrong? Thanks!

Here’s my code:

using ITensors, ITensorMPS, Statistics, Plots, HDF5, Printf
function hamiltonian(os::OpSum, J::Number, U::Number, mu::Number, Nsite::Int64)
  for i=2:(Nsite-1)
    os += J,"a†",(i+1),"a",i
    os += J,"a†",(i-1),"a",i
    os += U/2,"n",i,"n",i
    os += (-U/2-mu),"n",i
  end
  os += J,"a†",2,"a",1
  os += J,"a†",(Nsite-1),"a",Nsite
  for i=(1,Nsite)
      os += U/2,"n",i,"n",i
      os += (-U/2-mu),"n",i
  end
  return os
end
function hamiltonianSB(os::OpSum, J::Number, U::Number, mu::Number, Nsite::Int64)
  for i=2:(Nsite-1)
    os += J,"a†",(i+1),"a",i
    os += J,"a†",(i-1),"a",i
    os += U/2,"n",i,"n",i
    os += (-U/2-mu),"n",i
  end
  os += J,"a†",2,"a",1
  os += J,"a†",(Nsite-1),"a",Nsite
  for i=(1,Nsite)
      os += U/2,"n",i,"n",i
      os += (-U/2-mu),"n",i
  end
  os += (1E-2),"a†",1
  os += (1E-2),"a",1
  return os
end

mutable struct EnergyDiffObserver <: AbstractObserver
    energy_tol::Float64
    last_energy::Float64
    this_diff::Float64
 
    EnergyDiffObserver(energy_tol=0.0) = new(energy_tol,1000.0)
end

function ITensors.checkdone!(o::EnergyDiffObserver;kwargs...)
    sw = kwargs[:sweep]
    energy = kwargs[:energy]
    o.this_diff = abs(energy-o.last_energy)/abs(energy)
    if abs(energy-o.last_energy)/abs(energy) < o.energy_tol
      println("Stopping DMRG after sweep $sw")
      return true
    end
    # Otherwise, update last_energy and keep going
    o.last_energy = energy
    return false
end

n = 100
etol = 1E-6
d = 10
sites = siteinds("Boson", n, dim = d)
os = OpSum()
J = 0.1
U = 1
mu = 0.6
os = OpSum()
H_SB = hamiltonianSB(os, J, U, mu, n)
HMPO_SB = MPO(H_SB, sites)
psi0 = random_mps(sites; linkdims=10)
nsweeps = 10
maxdim = 10:20:220
cutoff = [1E-8]
obs = EnergyDiffObserver(etol)
@printf("calculating SB\n")
energy,psi_SB = dmrg(HMPO_SB, psi0;
nsweeps, maxdim, cutoff, observer=obs, outputlevel=1)
H = hamiltonian(os, J, U, mu, n)
HMPO = MPO(H, sites)
psi1 = psi_SB
nsweeps = 50
maxdim = 10:20:220
cutoff = [1E-8]
obs = EnergyDiffObserver(etol)
@printf("calculating no SB\n")
energy,psi = dmrg(HMPO,psi1;
nsweeps, maxdim, cutoff, observer=obs, outputlevel=1)
abs_a_expect = abs.(expect(psi, "a"; sites=1:n))
2

It’s difficult for us or anyone to say what you did “wrong” just based on your Hamiltonian and code. I say “wrong” in quotes because it’s possible that your code is working perfectly and these are the correct results. It’s important when using numerical methods like this to also consider the following things, which I’m sure you are aware of too:

  • what is the role of finite-size effects (e.g. what if you set n to smaller and larger values)?
  • what is the role of boundary effects? have you plotted the expected value \langle a_j \rangle and looked at it? does the shape look plausible and converged, symmetric?
  • are your results converged as a function of number of sweeps, cutoff, and max bond dimension (maxdim)? do they change if you raise or lower these?