Coupling fermionic operators with bosonic operators

DMRG and Numerical Methods
,
1

Hi!

We are working with a model that couples fermionic operators with bosonic operators such that our Hamiltonian looks something like:

H = (a^{\dagger}+a)\sum_i c^{\dagger}_i c_{i+1} +h.c.

Every time one of our fermions hop, it releases a global photon that can talk to each of the sites. We can truncate our bosons to a two-level system in which we can write our bosonic operators in terms of the spin raising and lowering operators changing our Hamiltonian to:

H = (S^{+}+S^-)\sum_i c^{\dagger}_i c_{i+1} +h.c.

When we have attempted to write a script for this that couples the bosonic operator to the fermionic operator, we get an error

I = (dim=3|id=910|"l=2,Link") <Out>
  1: 2 QN()
  2: 1 QN({"Nf",0,-1},{"Sz",2})
Q = QN({"Nf",-1,-1},{"Sz",1})
From line 683, file index.cc

Index does not contain given QN block.

Index does not contain given QN block.
Abort trap: 6

Here is an example of the code we are using in C++

#include "itensor/all.h"
#include <vector>
#include <ostream>
#include <fstream>
#include <iostream>
#include <math.h> 
#define _USE_MATH_DEFINES

using namespace itensor;
using namespace std;

int main()
{
    int N = 21;
    auto sites = tJ(N,{"ConserveQNs=",true});

    float t2 = 1;
    float t1 = 0;
    float g = M_PI/2;
    std::complex<double> Opera(cos(g), sin(g));
    float mu0 = .0,
    mu1 = .5;

    auto ampo = AutoMPO(sites);
    auto state = InitState(sites);

    //Hopping
    for(int j = 2; j <= N-1; ++j)
        {
            
            // In unit cell hopping
            {
                ampo += -t1,"Cdagup",j+1,"Cup",j,"S+",1;
                ampo += -t1,"Cdagup",j,"Cup",j+1,"S+",1;
                ampo += -t1,"Cdagup",j+1,"Cup",j,"S-",1;
                ampo += -t1,"Cdagup",j,"Cup",j+1,"S-",1;
            }
            // Out of unit cell hopping
            else 
            {
                ampo += -t2,"Cdagup",j+1,"Cup",j,"S+",1;
                ampo += -t2,"Cdagup",j,"Cup",j+1,"S+",1;
                ampo += -t2,"Cdagup",j+1,"Cup",j,"S-",1;
                ampo += -t2,"Cdagup",j,"Cup",j+1,"S-",1;
            }
        }

    auto H = toMPO(ampo);

    for(auto j : range1(N-1))
    {
        auto k = j+1;
        state.set(k,(k%2==1 ? "Up" : "Emp"));
    }
    state.set(1,"Up");
    auto psi0 = MPS(state);


    for(int j = 1; j<= N; ++j)
            {
            psi0.position(j);

            Real Nt = elt(psi0(j)
                       * op(sites,"Ntot",j)
                       * dag(prime(psi0(j),"Site")));

            println("Ntot",j," = ",Nt);

            }

    auto sweeps = Sweeps(5);
    sweeps.maxdim() = 10,20,50,100;
    sweeps.noise() = 1E-7,1E-8,1E-10;
    sweeps.cutoff() = 1E-5,1E-6,1E-7,1E-8;
    sweeps.niter() = 2,5,10,15,25;

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

        // Finding the inner product of N
    for(int j = 1; j<= N; ++j)
            {
            psi.position(j);

            Real Nt = elt(psi(j)
                       * op(sites,"Ntot",j)
                       * dag(prime(psi(j),"Site")));


            println("N",j," = ",Nt);

            }
    
    return 0;
}

Eventually, we would would like to go beyond the two level system and I think we would need to use something like our own custom siteset or use MixedSiteSet but this is our current roadblock. Any help on how we can do this would be greatly appreciated!

Thanks!

Alex Vargas