How to construct the square term in Hamiltonian

Question

Hello,
I am calculating a Bose-Hubbard model with a square term in its Hamiltonian:@@\widehat{H}=\sum_j{\left( \widehat{a}_{j}^{\dagger 2}+\widehat{a}_{j}^{2} \right)}@@ where @@\widehat{a}_{j}@@ is the annihilation operator of boson.
When I used code

ampo += 1, "Adag", i, "Adag", i;
ampo += 1, "A", i, "A", i;

with lattice class 'Boson' and run the code, the terminal gave feedback telling

  I = (dim=5|id=964|"l=1,Link") <Out>
  1: 3 QN()
  2: 1 QN({"Nb",-1})
  3: 1 QN({"Nb",1})
Q = QN({"Nb",-2})
From line 683, file index.cc
Index does not contain given QN block.
Aborted (core dumped)

which I also met when I used code like

ampo += 1, "S+", i, "S+", i;
ampo += 1, "S-", i, "S-", i;

in the spin system. I think I have to prime the operator to make them contract? But I have no idea how to do it in autoMPO class.
Could you please tell me how can I construct a Hamiltonian containing a term like this? Thank you.

Comments

I haven't messed with Bosons in ITensor yet, but you could make an alternate site set with the operators you'd like by adapting from the boson.h file (github version here: https://github.com/ITensor/ITensor/blob/v3/itensor/mps/sites/boson.h)

This can be done by copying the contents of the boson.h file and adding something like the following to the list of operators in your siteset file.

else
    if(opname == "Asq" || opname == "Asq")
            {
            for(auto n : range1(maxOcc-1))
                {
                Op.set(s=2+n,sP=n,std::sqrt(n*(n+1)));
                }
            }
else
    if(opname == "ADagsq" || opname == "ADagsq")
            {
            for(auto n : range1(maxOcc-1))
                {
                Op.set(s=n,sP=n+2,std::sqrt(n*(n+1)) );
                }
            }

The maxOcc-1 is because the onsite space can hold maxOcc+1 values for the number of Bosons (0, 1, 2, ..., maxOcc). Adjusting the upper end of that for loop should allow you to avoid if-thens and the ITensor code will assume all other inputs to give a resulting value of 0. You should double check all of the quantities above -  I'm writing this in haste.

You'll need to rename some things to avoid double-naming issues.

This is 100% untested by me and 100% *not* guaranteed, but something similar to this should work if there's no quick way to have pair creation/annihilation from the A/ADag operators that have already been written. Also, once you've made your own siteset (and debugged the naming errors!), it's far less intimidating to do it a second time for another model.

Cheers,
Jared

---

Dear Jared,
Thank you so much for your answer.
I tried your way carefully but it still replied that "Index does not contain given QN block.". I though it may came from the maximum occupancy because the default maximum occupancy is 1, so I set the maximum occupancy more than two, but it have the same errors...
The terminal told me that QN index only have Nb = +1/-1 while my QN index with this Hamiltonian term have Nb=-2. Maybe I need to change my QN index to fixed it?

My code is as follow
#include "itensor/all.h"

using namespace itensor;
using std::vector;

int main()
{
    int L = 8;
    auto sites = Boson(L,{"MaxOcc=",4,"ConserveQNs",true});
    auto sweeps = Sweeps(5);
    sweeps.maxdim() = 40, 80, 400, 400, 800;
    sweeps.cutoff() = 1E-16;

    Real miu = 1.0;
    Real U = 1.0;
    Real J = 1.0;
    Real V = 1.0;
    Real em = 1.0;

    auto ampo = AutoMPO(sites);
    for (int i = 1; i < L; i += 1)
    {
        ampo += -miu-U/2, "N", i;
        ampo += U/2, "N", i, "N", i;
        ampo += -J, "A", i, "Adag", i+1;
        ampo += -J, "Adag", i, "A", i+1;
        ampo += V, "N", i, "N", i+1;
        ampo += -em/2, "ADagsq", i;
        ampo += -em/2, "Asq", i;
    }
    ampo += -J, "A", 1, "Adag", L;
    ampo += -J, "Adag", 1, "A", L;
    ampo += V, "N", 1, "N", L;
    auto H = toMPO(ampo);

    auto state = InitState(sites);
    for (int i : range1(L))
    {
        if (i%2 == 1) state.set(i,"2");
        else state.set(i,"2");
    }
    auto psi = randomMPS(state);
    
    auto [energy,psi0] = dmrg(H,psi,sweeps,"Quiet");
    printfln("-----Energy = %.2f-----", energy/L);

return 0;
}

And the error report is as follow
I = (dim=5|id=69|"l=1,Link") <Out>
  1: 3 QN()
  2: 1 QN({"Nb",-1})
  3: 1 QN({"Nb",1})
Q = QN({"Nb",-2})
From line 683, file index.cc

Index does not contain given QN block.

Index does not contain given QN block.
Aborted (core dumped)

---

You're welcome, Zhiyao.

If you look near the top of the boson.h (or your custom siteset if you copied over from boson.h). you can see this:
        auto conserveQNs = args.getBool("ConserveQNs",true);
        auto conserveNb = args.getBool("ConserveNb",conserveQNs);

In file that you compile and run for this Hamiltonian, you tell the code to conserve quantum numbers but the code above (from the siteset with your pair-exchange operators) says that the default value for whether boson number should be conserved is the same as the value for conserveQNs. So.. Your Hamiltonian explicitly breaks number conservation via pair exchange, but you've told ITensor that you expect your Hamiltonian should conserve bosonic number. When it attempts to map outside of the fixed number space  ITensor produces an error then aborts.

In your main code, change
    auto sites = Boson(L,{"MaxOcc=",4,"ConserveQNs",true});
to
    auto sites = Boson(L,{"MaxOcc=",4,"ConserveQNs",true,"ConserveNb",false});

Then you can continue debugging.

Jared

---

Dear Jared,
Thank you  so much for your quick reply! It help me a lot!
I understant my fault now. Thank you for your guidance.

---

You're welcome!

Jared

Answer 1

Answered in the comments above. Adapting from the boson.h file (github version here: https://github.com/ITensor/ITensor/blob/v3/itensor/mps/sites/boson.h)
allows the desired operators to be defined in a custom siteset.

Copying the contents of the boson.h file and adding something like the following to the list of operators into the custom siteset file:

else
    if(opname == "Asq" || opname == "Asq")
            {
            for(auto n : range1(maxOcc-1))
                {
                Op.set(s=2+n,sP=n,std::sqrt(n*(n+1)));
                }
            }
else
    if(opname == "ADagsq" || opname == "ADagsq")
            {
            for(auto n : range1(maxOcc-1))
                {
                Op.set(s=n,sP=n+2,std::sqrt(n*(n+1)) );
                }
            }

The maxOcc-1 is because the onsite space can hold maxOcc+1 values for the number of Bosons (0, 1, 2, ..., maxOcc). Adjusting the upper end of that for loop should allow you to avoid if-thens and the ITensor code will assume all other inputs to give a resulting value of 0. I still haven't double-checked the quantities like the square roots, should double check all of that before using it.

Also note that for this model, number conservation is being broken, so that either needs to be removed from the custom siteset file or set to false in the declaration of 'sites' to avoid an error being thrown.

Cheers,
Jared