Simulating the driven-dissipative Bose-Hubbard model

[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?

 

[eljose79]

Thank you for reaching out and sharing your question with us. As a scientist familiar with using Mathematica for simulations, I can offer some suggestions for simulating the driven-dissipative Bose-Hubbard model.

Case 1: Working in a specific basis, such as the Fock or occupation basis, can be a good starting point for simulations. However, as you have observed, it may not always be the most efficient or stable approach. My suggestion would be to try using different bases and see which one yields the best results for your specific system. For example, the momentum basis may be more suited for studying certain properties of the system, while a mean-field approach may be better for others. It may also be helpful to consult literature or speak with colleagues who have worked on similar systems to get a better idea of which basis to use.

In addition, you may also want to experiment with different numerical methods for propagating the density matrix, such as the Runge-Kutta algorithm you mentioned or other techniques like the time-dependent variational principle. Each method has its own advantages and limitations, so it would be beneficial to explore and compare their performances for your specific system.

Case 2: For simulations in two dimensions, it may be necessary to adopt a mean-field description due to the increased complexity of the system. However, as with the single dimension case, it would be beneficial to try out different approaches and see which one works best for your specific system.

In summary, there is no one "most convenient" way of simulating the driven-dissipative Bose-Hubbard model in Mathematica. It ultimately depends on the specific properties and dynamics of your system. I hope these suggestions can help guide you in finding the most suitable approach for your simulations. Best of luck with your research!