spinless bosons

Question

Hello!

How one can construct spinless boson basis?
Here there is only hard-core bosons http://itensor.org/docs.cgi?page=classes/spinless&vers=cppv2

And here there is only spinfull bosons http://itensor.org/docs.cgi?vers=cppv2&page=classes/hubbard

I want to consider 1D bose-hubbard model.

thank you!

Answer 1

Hi,

Spinless boson is tricky, but there are many ways to deal with it. See section V of U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).

The easiest thing one can do is just cutting the local Hilbert space to some max occupancy, e. g. double occupancy, as long as one is only interested in low energy physics. You just need to define your own siteset similar to the one already defined in the ITensor library, with the index changed to be e. g.

s = IQIndex{
                nameint("site=",n),
                Index(nameint("Emp",n),1,Site),QN("Nb=",0),
                Index(nameint("N1",n),1,Site), QN("Nb=",1),
                Index(nameint("N2",n),1,Site), QN("Nb=",2)
                };

and changes the definition of the site operator accordingly.

Comments

thx. could u please write more in detail "...changes the definition of the site operator accordingly."

---

You just need to write a file similar to e. g. spinless.h, which is already defined in the ITensor library.
What you need to change:
1. Those "s = IQIndex{ ... }" lines.
2. The function "IQIndexVal state ...": now you have three cases in the if...else..., i.e. state == "0", state == "1", state == "2".
3. The function "IQTensor op...": since you are doing bosonic system, you no longer need those "C", "Cdag", "F"... and you need to change all other definitions of the bosonic operators. For example, "A" is changed to be:
        if(opname == "A")
            {
            Op.set(N1,EmpP, 1);
            Op.set(N2,N1P, sqrt(2));
            }
with EmpP=sP(1), N1=s(2), N1P=sP(2), N2=s(3).

---

Could you please see my modification of spinless.h file?

Answer 2

Could someone see my modification of spinless.h file to simulate spinless bosons in 1D Bose-Hubbard model?

#ifndef __ITENSOR_SPINLESS_H
#define __ITENSOR_SPINLESS_H
#include "itensor/mps/siteset.h"

namespace itensor {

class SpinlessBosons;

using SpinlessBos = BasicSiteSet<SpinlessBosons>;

class SpinlessBosons
    {
    IQIndex s;
    public:

    SpinlessBosons() { }
    SpinlessBosons(IQIndex I) : s(I) { }
    SpinlessBosons(int n, Args const& args = Args::global())    
        {
        s = IQIndex{
            nameint("SpinlessBos",n),
            Index(nameint("Emp",n),1,Site),QN("Nb=",0),
            Index(nameint("N1",n),1,Site), QN("Nb=",1),
            Index(nameint("N2",n),1,Site), QN("Nb=",2)
            };

        }

    IQIndex
    index() const { return s; }

    IQIndexVal
    state(std::string const& state)
        {

        if(state == "Emp" || state == "0") 
            {
            return s(1);
            }   
        else 
        if(state == "N1" || state == "1") 
            {
            return s(2);
            }
        else 
        if(state == "N2" || state == "2") 
            {
            return s(3);
            }
        else
            {
            Error("State " + state + " not recognized");
            }
        return IQIndexVal{};
        }

    IQTensor
    op(std::string const& opname,
       Args const& args) const
        {
        auto sP = prime(s);
        auto Emp  = s(1);
        auto EmpP = sP(1);
        auto N1  = s(2);
        auto N1P = sP(2);
        auto N2 = s(3);
        auto N2P = sP(3);


        auto Op = IQTensor(dag(s),sP);

        if(opname == "N" || opname == "n")
            {
            Op.set(N1,N1P,1);
            Op.set(N2,N2P,1);
            }
        else
        if(opname =="A")
            {
            Op.set(N1, EmpP,1); 
            Op.set(N2,N1P,1);
            }
        else
        if(opname == "Adag")
            {
            Op.set(Emp, N1P,1); 
            Op.set(N1,N2P,1);   
            }

        else
            {
            Error("Operator \"" + opname + "\" name not recognized");
            }

        return Op;
        }
    };

} //namespace itensor

#endif

the same code with comments (i.e. the code in spinless.h is commented out ) for comparison with the file spinless.h

#ifndef __ITENSOR_SPINLESS_H
#define __ITENSOR_SPINLESS_H
#include "itensor/mps/siteset.h"

namespace itensor {

//class SpinlessSite;
class SpinlessBosons;
//using Spinless = BasicSiteSet<SpinlessSite>;
using SpinlessBos = BasicSiteSet<SpinlessBosons>;

//class SpinlessSite
class SpinlessBosons
    {
    IQIndex s;
    public:

//    SpinlessSite() { }
    SpinlessBosons() { }
//    SpinlessSite(IQIndex I) : s(I) { }
    SpinlessBosons(IQIndex I) : s(I) { }
//    SpinlessSite(int n, Args const& args = Args::global())
    SpinlessBosons(int n, Args const& args = Args::global())    
        {
        s = IQIndex{
            nameint("SpinlessBos",n),
            Index(nameint("Emp",n),1,Site),QN("Nb=",0),
            Index(nameint("N1",n),1,Site), QN("Nb=",1),
            Index(nameint("N2",n),1,Site), QN("Nb=",2)
            };
        /*
        auto conserve_Nf = args.getBool("ConserveNf",true);
        auto oddevenupdown = args.getBool("OddEvenUpDown",false);

        if(!oddevenupdown) //usual case
            {
            auto q_occ = QN("Nf=",1);
            if(not conserve_Nf) q_occ = QN("Pf=",1);


            s = IQIndex{nameint("Spinless ",n),
                Index(nameint("Emp ",n),1,Site),QN(),
                Index(nameint("Occ ",n),1,Site),q_occ};

            }
        else
            {
            QN q_occ;
            if(n%2==1) q_occ = QN("Sz",+1,"Nf=",1);
            else       q_occ = QN("Sz",-1,"Nf=",1);
            s = IQIndex{nameint("Spinless ",n),
                Index(nameint("Emp ",n),1,Site),QN(),
                Index(nameint("Occ ",n),1,Site),q_occ};


            }*/ 
        }

    IQIndex
    index() const { return s; }

    IQIndexVal
    state(std::string const& state)
        {

        if(state == "Emp" || state == "0") 
            {
            return s(1);
            }   
        else 
        if(state == "N1" || state == "1") 
            {
            return s(2);
            }
        else 
        if(state == "N2" || state == "2") 
            {
            return s(3);
            }
        else
            {
            Error("State " + state + " not recognized");
            }
        /*
        if(state == "Emp" || state == "0") 
            {
            return s(1);
            }   
        else 
        if(state == "Occ" || state == "1") 
            {
            return s(2);
            }
        else
            {
            Error("State " + state + " not recognized");
            }*/
        return IQIndexVal{};
        }

    IQTensor
    op(std::string const& opname,
       Args const& args) const
        {
        auto sP = prime(s);
        auto Emp  = s(1);
        auto EmpP = sP(1);
        auto N1  = s(2);
        auto N1P = sP(2);
        auto N2 = s(3);
        auto N2P = sP(3);

        //auto Occ  = s(2);
        //auto OccP = sP(2);

        auto Op = IQTensor(dag(s),sP);

        /*
        if(opname == "N" || opname == "n")
            {
            Op.set(Occ,OccP,1);
            }
        else
        if(opname == "C")
            {
            Op.set(Occ,EmpP,1);
            }
        else
        if(opname == "Cdag")
            {
            Op.set(Emp,OccP,1);
            }
        else
        */  
        if(opname == "N" || opname == "n")
            {
            Op.set(N1,N1P,1);
            Op.set(N2,N2P,1);
            }
        else
        if(opname =="A")
            {
            Op.set(N1, EmpP,1); 
            Op.set(N2,N1P,1);
            }
        else
        if(opname == "Adag")
            {
            Op.set(Emp, N1P,1); 
            Op.set(N1,N2P,1);   
            }
        /*
        if(opname == "A")
            {
            Op.set(Occ,EmpP,1);
            }   
        else
        if(opname == "Adag")
            {
            Op.set(Emp,OccP,1);
            }
        else
        if(opname == "F" || opname == "FermiPhase")
            {
            Op.set(Emp,EmpP,1);
            Op.set(Occ,OccP,-1);
            }
        else
        if(opname == "projEmp")
            {
            Op.set(Emp,EmpP,1);
            }
        else
        if(opname == "projOcc")
            {
            Op.set(Occ,OccP,1); 
            }
        */  
        else
            {
            Error("Operator \"" + opname + "\" name not recognized");
            }

        return Op;
        }
    };

} //namespace itensor

#endif

Comments

1) I wonder about

auto conserve_Nb = args.getBool("ConserveNb",true);

should I add something like that for conserving the number of particles?
Will the number of particles be conserved in the code written above?

2) I don't really understand what is

 auto oddevenupdown = args.getBool("OddEvenUpDown",false);

and should I insert something similar in my file

3) Is it Ok to write
Op.set(N1,N2P,1);   
after    if(opname == "Adag")?
the action of Adag on the sate with 1 particle should give the state with 2 particles multiplied by sqrt(2).

---

Hi,
You can just ignore 1) and 2). For 3), You are right. I have corrected this typo in my answer.
Your first code looks good after changing "Op.set(N1,N2P,1)" to "Op.set(N1,N2P,std::sqrt(2))" and the corresponding one in opname ==  "A" (the correct is "Op.set(N2,N1P,std::sqrt(2))") and also in opname == "N" (the correct is "Op.set(N2,N2P,2)").
As Miles suggested below, you can test with max occupancy defined to be 2, 3, 4... to see if that makes any difference in your result.

Answer 3

Hi AbraDabra,
This is a site set type we should offer. We're working on it!

But in the meantime, yes, you can create your own kind of site set for this case by copying one of the existing ones (such as "Spinless") and creating a new file and a new site set (perhaps boson.h and "class BosonSite" and "class Boson").

Though Mingru is correct that sometimes just having a maximum site occupancy of 1 is enough, in general that's not the case, very much depending on the particular Hamiltonian.

But for many cases, you can truncate the maximum occupancy to, say, 3 or 4 without incurring much error. You should raise and lower this maximum occupancy and check that your results are converged with respect to it.

For the operators you'll need to define, here is an example for maximum occupancy of 3:

ITensor
op(std::string const& opname,
   Args const& args) const
    {
    auto sP = prime(s); //here "s" is the site index stored in this BosonSite object

    IQIndexVal Emp(s(1)),
               EmpP(sP(1)),
               One(s(2)),
               OneP(sP(2)),
               Two(s(3)),
               TwoP(sP(3)),
               Three(s(4)),
               ThreeP(sP(4));

    if(opname == "N" || opname == "n")
        {
        Op(One,OneP) = 1;
        Op(Two,TwoP) = 2;
        Op(Three,ThreeP) = 3;
        }
    else
    if(opname == "A")
        {
        Op(One,EmpP) = 1;
        Op(Two,OneP) = std::sqrt(2);
        Op(Three,TwoP) = std::sqrt(3);
        }
    else
    if(opname == "Adag")
        {
        Op(Emp,OneP) = 1;
        Op(One,TwoP) = std::sqrt(2);
        Op(Two,ThreeP) = std::sqrt(3);
        }

If you'd like to post your code when you're done with a first draft or have some more specific questions please just post a comment or a new question -

Best regards,
Miles