How to add Langevin terms to the semiclassical Bose-Hubbard model?

Question

I would like to add Langevin terms to the Hamilton equations of motion of the semiclassical Bose-Hubbard model.

Here's what I have:

I start with the standard example of Brownian motion, a particle in a potential. Its Hamilton function reads:

$H = \frac{1}{2m} p^2 + V\left(q\right)$,

The corresponding Hamilton equations of motion (EoM) read:

$\dot{p}=-\partial_{q}H=-\partial_{q}V\left(q\right)$

$\dot{q}=\partial_{p}V\left(q\right)=\frac{p}{m}$

One can convert these coupled differential equations of first order into a single differential equation of second order:

$\ddot{q}=-\frac{1}{m}\partial_{q}V\left(q\right)$

and rewrite them as

$m\ddot{q}=-\partial_{q}V\left(q\right)$.

In this form one can add Langevin terms (see the Wikipedia entry on Langevin dynamics) and one obtains:

$m\ddot{q}=-\partial_{q}V\left(q\right)-\gamma m \dot{q} + \sqrt{2\gamma m k_\text{B}T} \;\xi\left(t\right)$,

where $\gamma$ is the damping (free parameter), and $\xi\left(t\right)$ a delta-correlated stationary Gaussian process with zero-mean, satisfying:

$\left\langle \xi\left(t_{i}\right) \xi\left(t_{j}\right)\right\rangle=\delta\left(t_{i}-t_{j}\right)$.

In order to solve this numerically with an SDE solver (e.g., Heun scheme), we need to write this as a system of two first-order differential equations:

$\dot{p}=-\frac{1}{m}\partial_{q}V\left(q\right)-\gamma \dot{q} + \sqrt{2\gamma m k_\text{B}T} \;\xi\left(t\right)$

$\dot{q}=\frac{p}{m}$

We can achieve the same for the $xy$-model (which is very similar to the semiclassical Bose-Hubbard model, my target Hamiltonian). Its Hamilton function reads

$H^{\text{xy}}=-\sum_{\left\langle ij\right\rangle}J_{ij}\cos\left(\theta_{i}-\theta_{j}\right)+\frac{1}{2}U\delta n_{i}^2$,

where the canonical conjugate variables are the onsite phase $\theta_{i}$ and the density fluctuations $\delta n_{i}$.

The corresponding EoM are

$\dot{\theta}_{i}=U\delta n_{i}$

$\delta\dot{n}_{i}=-\sum_{j\left(i\right)}J_{ij}\sin\left(\theta_{i}-\theta_{j}\right)$.

Formally, the term $\frac{1}{2}U\delta n_{i}^2$ is like a kinetic energy, with $U$ playing the role of an inverse mass. Thus by analogy we can add Langevin terms to the second equation, like in the previous example:

$\delta\dot{n}_{i}=-\sum_{j\left(i\right)}J_{ij}\sin\left(\theta_{i}-\theta_{j}\right)-\gamma \delta\dot{n}_{i} + \sqrt{2\gamma k_{\text{B}}T/U}\;\xi\left(t\right)$.

What I would like to have is a similar expression for the semiclassical Bose-Hubbard model.

I start with the semiclassical (with complex numbers instead of field operators) Bose-Hubbard Hamiltonian in coherent state representation,

$H^{\text{BHM}}\left(\psi^{\ast}_{i},\psi_{i}\right) = -\sum_{\left\langle ij\right\rangle}t_{ij}\left( \psi^{\ast}_{i}\psi_{j}+\text{c.c.}\right) + \frac{1}{2}U n^{2}_{i}$ where

$n_{i} = \psi^{\ast}_{i} \psi_{i}$,

and then transform that to coordinate-momentum representation, using

$\psi_{i}=\frac{1}{\sqrt{2}}\left(q_{i}+\imath p_{i}\right)$

$\psi^{\ast}_{i}=\frac{1}{\sqrt{2}}\left(q_{i}-\imath p_{i}\right)$,

$H^{\text{BHM}}\left(q_{i},p_{i}\right) = -\sum_{\left\langle ij\right\rangle}t_{ij}\frac{1}{2}\left(q_{i} p_{j} + q_{j} p_{i}\right)+ \frac{1}{2}U \frac{1}{4}\left(q^{2}_{i}+p^{2}_{i}\right)^2$

with EoM:

$\dot{q}_{i}=-\sum_{j\left(i\right)}t_{ij}p_{j}+\frac{1}{2}U\left(q^2_{i}+p^{2}_{i}\right)p_{i}$

$\dot{p}_{i}=+\sum_{j\left(i\right)}t_{ij}q_{j}-\frac{1}{2}U\left(q^2_{i}+p^{2}_{i}\right)q_{i}$

How does one correctly add Langevin terms?

Update (after Ted Pudlik's comment):

Following Ted Pudlik's suggestion, I write the semiclassical Bose-Hubbard Hamiltonian in density-phase notation:

$H^{\text{BHM}} \left(n_i,\theta_i\right)=-\sum_{\left\langle ij\right\rangle}\left(t_{ij}\sqrt{n_i}\sqrt{n_j}\mathrm{e}^{\imath\left(\theta_j-\theta_i\right)}+\text{c.c.}\right)+\frac{1}{2}U\sum_{i}n_i^2$

The corresponding Hamilton equations of motion are:

$\dot{\theta}_i=\partial_{n_i}H^{\text{BHM}}\left(n_i,\theta_i\right)=U n_i-\sum_{j\left(i\right)}\left(t_{ij}\frac{\sqrt{n_j}}{2\sqrt{n_i}}\mathrm{e}^{\imath\left(\theta_j-\theta_i\right)}+\text{c.c.}\right)$

$\dot{n_i}=-\partial_{\theta_i}H^{\text{BHM}}\left(n_i,\theta_i\right)=\sum_{j\left(i\right)}\left(\imath t_{ij}\sqrt{n_i n_j}\mathrm{e}^{\imath\left(\theta_j-\theta_i\right)}+\text{c.c.}\right)$

As in the XY-model I add the following terms to the density (not phase, as previously stated---see comment below) derivative:

$-\gamma n_i +\sqrt{2\gamma k_{\mathrm{B}}T/U}\;\xi\left(t\right)$

where I choose $\gamma$ according to my needs (e.g., smaller than smallest eigenfrequency of the system or overdamped).

What I'm still worried about is the $\frac{1}{2\sqrt{n_i}}$ term in the phase derivative: When the onsite density is very low compared to its neighbor, this term diverges. That is for example the case in the thermal cloud of an ultracold gas in a harmonic trap.

Is there a way to transform the EoM including the Langevin terms into $\left(q,p\right)$-representation or $\left(Re,Im\right)$-representation in order to avoid this?

Answer 1

I dug around in the literature a bit and found that this formulation (semiclassical Bose-Hubbard plus Langevin-type dissipation) has been studied before. Here is the relevant reference: http://arxiv.org/abs/1304.5071. What you are trying to do is derive their equation (9). You probably missed it because they refer to their model as the discrete nonlinear Schroedinger equation, but this is a different name for the same Hamiltonian (as mentioned, e.g., in this paper).

In your notation, their result reads (assuming all the hoppings $t_{ij} = t$ for simplicity),

$$ \imath \frac{d\psi_i}{dt} = (1 + \imath\gamma)\left(U |\psi_i|^2\psi_i - t(\psi_{i+1} + \psi_{i-1})\right) + \imath\gamma\mu\psi_i + \sqrt{\gamma T} \xi_i(t) $$

Appendix A contains a derivation of this equation for an arbitrary Hamiltonian. I did not read it closely enough to guarantee that it's correct, but the result looks sensible.

Unlike the amplitude-phase formulation, this one should not give you any nasty divergences for nearly-empty sites.

EDIT: As we discussed in the comments, I'm not too sure about the minus signs in this equation---as it stands, the well populations will diverge for positive $\gamma$! I think the equation should maybe read,

$$ \imath \frac{d\psi_i}{dt} = (1 + \imath\gamma)\left(U |\psi_i|^2\psi_i - t(\psi_{i+1} + \psi_{i-1})\right) - \imath\gamma\mu\psi_i + \sqrt{\gamma T} \xi_i(t) $$

Or maybe the idea is that $\mu > 0$ means that the reservoir is at a higher chemical potential, so particles keep entering the system? I would email the authors asking about these sign issues, I can't quite wrap my mind around this. Sorry!