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

• I think the coordinate-momentum representation is not the way to go: do an amplitude-phase decomposition instead ($\psi_i = \sqrt{n_i} \exp(\phi_i)$, with $n_i$, $\phi_i$ real). You then get a Hamiltonian just like the xy model, and can add Langevin terms in the same way. The amplitude-phase decomposition is also canonical (SE question). – Ted Pudlik Oct 19 '15 at 18:50
• Thank you very much for your suggestion. It has been very helpful. I have edited my question because in density-phase representation there's a divergent term when writing down the equations of motion. – Robert Nov 2 '15 at 18:32
• You add the Langevin term to the $\dot{n_i}$, not to the phase derivative $\dot{\theta_i}$, right? – Ted Pudlik Nov 2 '15 at 18:44
• The divergence of the phase derivative as $n_i \to 0$ reflects the physical fact that the phase of an empty mode is not well-defined. Also, note that the semiclassical approximation is only valid in the limit of large occupation numbers: if the onsite density is very low, it's not the right model to use! – Ted Pudlik Nov 2 '15 at 18:45
• Thanks, that's a point I hadn't thought about. In a typical optical lattice setup that involves a harmonic trap, densities in the center — the condensate part — seem to be high enough to justify the use of a semi-classical model, but outside in the thermal cloud they are not that high anymore. If thermal fluctuations are high enough though, one can hope they destroy all quantum effects that are not captured by the model. To avoid numerical troubles it would be nice not to have a term that could suddenly become very large. – Robert Nov 3 '15 at 0:14

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!

• It didn't even occur to me to look for the NLSE. I somehow associate it with photonics and continuous systems (Gross-Pitaevskii equation). The above equation seems to work (integration is stable) after changing the sign in $\left(1+\imath\gamma\right)$ to $\left(1-\imath\gamma\right)$. I think that's because in the paper they don't take the convention that the tunneling term has a negative sign in front of it. The derivation in the appendix is also very insightful. Thanks a lot for your help! – Robert Nov 3 '15 at 21:52
• The sign of the $\imath \gamma\mu\psi_i$ term should perhaps also be negative. In general, the sign convention used in this paper is a little confusing to me. (Shouldn't the equation of motion be $\imath \dot{z_i} = \partial H/\partial z_i^*$? They have an extra minus sign.) But it's clearly a relevant reference! I recommend you send an email to the authors asking if your translation of their equation into your notation is correct, or if there are sign issues. – Ted Pudlik Nov 3 '15 at 22:27