How to add Langevin terms to the semiclassical Bose-Hubbard model? 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?
 A: 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!
