# whether itensor can realize cavity induced model?

I want to know whether dmrg can realize hamiltonian contains:
$$(b^{\dag}+b)\sum{i,\sigma}n{i,\sigma}$$
b is the boson(cavity) and n is the number of electron at site i. b is interacted with all electron site. I don't know whether dmrg can calculate this .

Hi Weijie,
Yes, DMRG and ITensor can be used to study a Hamiltonian of this type (and any Hamiltonian which is a sum of local terms on a lattice).

But special care has to be taken with bosons, as you know, since one has to make their Hilbert space finite in some way. So you could use the Boson site set in ITensor and adjust the MaxOcc parameter to make sure your results are converged with the maximum boson occupancy per site.

commented by (250 points)
Hi miles
I find one example.The hamiltonian is H=He+Hph+He-ph
$$H_{e}=-t\sum_{i,\sigma}(c^{\dag}_{i\sigma}+h.c.)+U\sum_{i}n_{i\uparrow}n_{i\downarrow} \\ Hph=\Omega b^{dag}b \\ He-ph=-\frac{g}{\sqrt(N)}(b^{\dag}+b)\sum_{i,\sigma}(-1)^{i}n_{i\sigma}$$
And I write a code. When I execute the code, the energy of every iteration is not always descending. It will be fluctuating, as -10,-11,-10,-9.... I try different parameters. However the result is same.
#include "itensor/all.h"
using namespace itensor;

using Holstein = MixedSiteSet<ElectronSite,BosonSite>;

int
main()
{
auto N = 20;

auto U = 0.5;
auto J = 1.0;
auto g = 1.0;
auto sigma = 1.0;

auto sites = Holstein(N,{"ConserveNf=",true,
"ConserveNb=",false,
"MaxOcc=",50});

auto ampo = AutoMPO(sites);
for(int j=1;j  <= N-2; j+=2)
{
ampo += -J,"Cdagup",j+2,"Cup",j;
ampo += -J,"Cdagup",j,"Cup",j+2;
ampo += -J,"Cdagdn",j+2,"Cdn",j;
ampo += -J,"Cdagdn",j,"Cdn",j+2;
}
ampo += -J,"Cdagup",1,"Cup",19;
ampo += -J,"Cdagup",19,"Cup",1;
ampo += -J,"Cdagdn",1,"Cdn",19;
ampo += -J,"Cdagdn",19,"Cdn",1;
for(int j=1;j < N; j += 2)
{
ampo += U,"Nup",j,"Ndn",j;
}
for(int j=1;j<N;j +=4)
{
}
for(int j=3;j<N;j +=4)
{
}
auto H = toMPO(ampo);

auto sweeps = Sweeps(20);
//very important to use noise for this model
sweeps.noise() = 1E-6,1E-6,1E-8,1E-12;
sweeps.maxdim() = 10,20,100,100,200;
sweeps.cutoff() = 1E-10;

auto state = InitState(sites);
for(auto j : range1(N))
{
if(j%4==1) state.set(j,"Up");
else if(j%4==3) state.set(j,"Dn");
else state.set(j,"0");
}
auto psi0 = MPS(state);

auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet=",true});

printfln("Ground State Energy = %.12f",energy);
commented by (44.9k points)
Hi, for cases like this I think it is best to compute local measurements (like of the density of particles) and print them out to a file, then plot the results to visualize what is going on. It could be that your calculation is getting stuck due to the details of the physics.

Miles