Simulating the driven-dissipative Bose-Hubbard model

1. Mar 1, 2015

matteo137

I'm trying to simulate with Mathematica the driven-dissipative Bose-Hubbard model, and I would like to know what is the most convenient way of doing it.

The Hamiltonian of the system is the following:
$$H= \Delta \sum_j a^\dagger_j a_j + U \sum_j a^\dagger_j a^\dagger_j a_j a_j + J \sum_{\langle j,j^\prime\rangle} a^\dagger_j a_{j^\prime} + F_p \sum_j \left(a^\dagger_j + a_j \right)$$
where $\Delta, U, J, F_p$ are real constants, and the index $j$ goes from 1 to M.
$F_p$ is the amplitude of the driving field, and the dissipations are introduced through a master equation

$$\dfrac{d\rho}{\text{dt}} = - i \left[ H,\rho\right] + \dfrac{\kappa}{2} \sum_j \left(2 a_j \rho a^\dagger_j -a^\dagger_j a_j \rho - \rho a^\dagger_j a_j \right)$$
where $\kappa$ is the real constant.

Case 1: for the case of $M=3$ sites (trimer), I described each site in the Fock basis. i.e. with the vectors $(1,0,0,...)=\vert 0\rangle$, $(0,1,0,...)=\vert 1\rangle$, ... up to a cutoff $\vert n_{\text{max}}\rangle$. Starting from this basis, I expressed the operators $a$, $a^\dagger$ as matrices, and I built the Hamiltonian matrix, of size $(n_{\text{max}})^M$, using Kronecker products of the kind $a_1 = a\otimes 1\otimes 1$.
The matrix $\rho(t)$ is then obtained propagating a starting $\rho(0)$ with the Runge-Kutta algorithm.

This method does not seems very efficient and stable...
would it be better to work in the "occupation" basis? i.e. writing the operators in the basis built from the vectors $(0,0,0)=\vert 0,0,0\rangle$, $(1,0,0)=\vert 1,0,0\rangle$, ...
Or, for example, in the momentum basis, mean-field, ...?

Case 2: what about in two dimensions? ...is it better to work with a mean-field description?

2. Mar 6, 2015