\[ H = -t \sum_{< i,j >} \left( a_i^+ a_j + b_i^+ b_j + H.c. \right) + \frac{U}{2} \sum_i [a_i^+ a_i (a_i^+ a_i - 1) + b_i^+ b_i (b_i^+ b_i - 1)] + U_{AB} \sum_i a_i^+ a_i b_i^+ b_i, \]

where \(U\) and \(U_{AB}\) are the inter- and intrastate onsite interactions, respectively. I have implemented this as its own `CouplingMPOModel` class, while heeding some of the things I've learned from studying this topic on combining multiple charges on the same site.

Essentially I set up a superlattice with two sites in every unit cell: an _A_ site and a _B_ site. In the `init_lattice()` routine I use a `pairs` dictionary to get the geometry right:

Code: Select all

```
class TwoComponentBoseHubbardModel(CouplingMPOModel):
def init_sites(self, model_params):
n_max = model_params.get('n_max', 3)
filling = model_params.get('filling', 1)
conserve = model_params.get('conserve', 'N')
if conserve == 'best':
conserve = 'N'
if self.verbose >= 1.:
print(self.name + ": set conserve to", conserve)
siteA = BosonSite(Nmax=n_max, conserve=conserve, filling=filling)
siteB = BosonSite(Nmax=n_max, conserve=conserve, filling=filling)
multi_sites_combine_charges([siteA, siteB])
return [siteA, siteB]
def init_lattice(self, model_params):
pairs = {'onsite': [(0, 1, np.array([0]))], 'nearest_neighbors': [(0, 0, np.array([1])), (1, 1, np.array([1]))]}
L = model_params.get('L', 10)
sites = self.init_sites(model_params)
lat = lattice.Lattice([L], sites, pairs = pairs)
return lat
def init_terms(self, model_params):
t = model_params.get('t', 1.)
U = model_params.get('U', 10.)
UAB = model_params.get('UAB', 0.)
mu = model_params.get('mu', 0)
for u in range(len(self.lat.unit_cell)):
self.add_onsite(-mu - U / 2., u, 'N')
self.add_onsite(U / 2., u, 'NN')
for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
self.add_coupling(-t, u1, 'Bd', u2, 'B', dx, plus_hc=True)
for u1, u2, dx in self.lat.pairs['onsite']:
self.add_coupling(UAB, u1, 'N', u2, 'N', dx)
```

It does seem to me, though, that it runs a little inefficiently, so perhaps my choice of geometry is suboptimal. I need a bond dimension of about 300 to get a converged result for 20 sites. I must say I've never tried to use MPS algorithms on the Bose-Hubbard model before, so it may just be my inexperience (in which case, too bad!) On the other hand I wonder if the MPS representation is somehow more efficient if I use a ladder geometry, with the _A_ and _B_ sites being on either side of every rung.