Changing Bose-Hubbard model

[mariocat]

Hello,
I would like to change the Bose-Hubbard Model such that in the first sum in can add \(\exp{(j*\pi*n_{i})}\), where ni is the number operator acting on the i site. How can I do this?
I was thinking about working on this part of the classe:

Python: Select all

def init_terms(self, model_params):
        # 0) Read and set parameters.
        t = model_params.get('t', 1., 'real_or_array')
        U = model_params.get('U', 0., 'real_or_array')
        V = model_params.get('V', 0., 'real_or_array')
        mu = model_params.get('mu', 0, 'real_or_array')
        phi_ext = model_params.get('phi_ext', None, 'real')
        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']:
            if phi_ext is None:
                hop =self.add_coupling(-t,dx,[0,expm(1j*np.pi*'N')])
            else:
                hop = self.coupling_strength_add_ext_flux(-t, dx, [0, 2 * np.pi * phi_ext])
            self.add_coupling(hop, u1, 'Bd', u2, 'B', dx, plus_hc=True)
            self.add_coupling(V, u1, 'N', u2, 'N', dx)

[Jakob]

The term you want is called parity or 'P' in the BosonSite, with definition \(P = 1 - 2 (n \text{ mod } 2) = (-1)^n = \mathrm{exp}(\mathrm{i} \pi n)\).
So for H = \(\sum_i \gamma_i * \mathrm{exp}(\mathrm{i} \pi n_i) + H_\text{other}\) you could do something like

Python: Select all

def init_terms(self, model_params):
    ...  # other terms
    gamma = model_params.get('gamma', 0, 'real_or_array')
    for u in range(len(self.lat.unit_cell)):
        self.add_onsite(gamma, u, 'P')

[Jakob]

If you want to extend the BH model without rewriting everything, you can do sth like

Python: Select all

class MyExtendedBHModel(tenpy.BoseHubbardModel):
    def init_terms(self, model_params):
        super().init_terms(model_params)
        gamma = model_params.get('gamma', 0, 'real_or_array')
        for u in range(len(self.lat.unit_cell)):
            self.add_onsite(gamma, u, 'P')

which has all the Hamiltonian terms of the BH model, plus this new parity term

[mariocat]

Hello, thanks for the answer. Actually I do not want to add a new term in the Hamiltonian but simply modify the first one by putting the parity operator just after the sum symbol and before b_i+bj +h.c.
I dont know if I've been clear!

[Jakob]

Ah I see.
In that case you can replace the hopping

Python: Select all

for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
    self.add_coupling(-t, u1, 'Bd', u2, 'B', dx, plus_hc=True)

that generates \(\sum_{<i,j>} (b^\dagger_i b_j + hc) \)
with

Python: Select all

for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
    self.add_coupling(-t, u1, 'P Bd', u2, 'B', dx, plus_hc=True)

to get a coupling \(\sum_{<i,j>} (P_i b^\dagger_i b_j + hc) = \sum_{<i,j>} (\mathrm{e}^{\mathrm{i}\pi n_i} b^\dagger_i b_j + hc)\)

I have pushed some more documentation to main in https://github.com/tenpy/tenpy/commit/1 ... 538b74027e
to make this easier to see in the future.

[mariocat]

Thanks again! I do appreciate your help! Just one last thing: how would you implement something like $$-\tau \sum_{\langle i,j \rangle i<j} e^{i \alpha\, n_i}\left( b_i^\dagger b_j + b_j^\dagger b_i \right)$$ where alpha is just a generic number. So basically this accounts to taking the BH Hamiltonian and change the first term

[Jakob]

The operator you want is not one of the pre-defined operators of the BosonSite, so you would need to add it explicitly.
You probably want to do something like this:

Python: Select all

class MyExtendedBoseHubbarModel(tenpy.BoseHubbardModel):

    def init_sites(self, model_params):
        site = super().init_sites(model_params)
        alpha = model_params.get('alpha', np.pi, 'real')
        N_diag = np.arange(site.Nmax + 1)
        site.add_op('P_alpha', np.diag(np.exp(1.j * alpha * N_diag)), permute_dense=True)
        return site

    def init_terms(self, model_params):
        ...  # other terms
        tau = model_params.get('tau', 1, 'real_or_array')
        for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
            self.add_coupling(-tau, u1, 'P_alpha Bd', u2, 'B', dx)
            self.add_coupling(-tau, u1, 'P_alpha B', u2, 'Bd', dx)

Note that i did not verify carefully if this is correct, and you should do that, see
https://tenpy.readthedocs.io/en/latest/ ... ing-models

[mariocat]

Thanks, it works!

[mariocat]

Hey, still sorry for bothering!
Actually there's something I am not really getting:
1)

Python: Select all

 alpha = model_params.get('alpha', np.pi, 'real')
        N_diag = np.arange(site.Nmax + 1)
        site.add_op('P_alpha', np.diag(np.exp(1.j * alpha * N_diag)), permute_dense=True)

Could you please explain me what these 3 lines do?

2) If I use any value of alpha different from $$0$$ or $$2\pi$$ I get this error: WARNING:tenpy.linalg.krylov_based:poorly conditioned H matrix in KrylovBased! |psi_0| = 0.192107
Any clue about what it could mean?

[Jakob]

1)
- line before the three: initialize a BosonSite, with its opertors etc. next we modify it by adding the new operator
- extract alpha from the model_params, this is just a convenient way to tell the model the value for alpha
- N_diag is a helper for defining the operator. It is the diagonal of the number operator, i.e. N_diag = [0, 1, 2, ...]
- Form the operator we want \(\mathrm{exp}(\mathrm{i}\alpha n)\), where we now that since n is diagonal, we can just take np.exp on the diagonal. permute_dense=True ensures that you do not need to think about details of the basis order of the site. We then add that operator to the site and give it a name, st we can later reference it by that name.

2) \(\alpha=\pi\) should also make that go away, no? The Hamiltonian is not Hermitian for other values of alpha. tenpy does not generally support non-hermitian hamiltonians at the moment, in particular the lanczos eigensolver used in DMRG assumes a hermitian Hamiltonian.

Did you want to simulate a non-hermitian Hamiltonian?

If you wanted to do a term that is always hermitian
\(-\tau \sum_{\langle i,j \rangle}\left(e^{i \alpha\, n_i} b_i^\dagger b_j + e^{-i \alpha\, n_i} b_j^\dagger b_i \right) = -\tau \sum_{\langle i,j \rangle} \left(P_i b_i^\dagger b_j + hc \right)\)

you can do that instead by

Python: Select all

for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
    self.add_coupling(-t, u1, 'P_alpha Bd', u2, 'B', dx, plus_hc=True)