¿Cómo agregar términos de Langevin al modelo semiclásico de Bose-Hubbard?

Me gustaría agregar términos de Langevin a las ecuaciones de movimiento de Hamilton del modelo semiclásico de Bose-Hubbard.

Esto es lo que tengo:

Comienzo con el ejemplo estándar del movimiento browniano, una partícula en un potencial. Su función de Hamilton dice:

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

Las correspondientes ecuaciones de movimiento de Hamilton (EoM) dicen:

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

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

Se pueden convertir estas ecuaciones diferenciales acopladas de primer orden en una sola ecuación diferencial de segundo orden:

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

y reescribirlos como

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

De esta forma, se pueden agregar términos de Langevin (ver la entrada de Wikipedia sobre la dinámica de Langevin ) y se obtiene:

$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)$,

dónde$\gamma$es el amortiguamiento (parámetro libre), y$\xi\left(t\right)$un proceso gaussiano estacionario correlacionado delta con media cero, que satisface:

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

Para resolver esto numéricamente con un solucionador SDE (por ejemplo, el esquema de Heun), necesitamos escribir esto como un sistema de dos ecuaciones diferenciales de primer orden:

$\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}$

Podemos lograr lo mismo para el$xy$-modelo (que es muy similar al modelo semiclásico de Bose-Hubbard, mi hamiltoniano objetivo). Su función de Hamilton dice

$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$,

donde las variables canónicas conjugadas son la fase in situ$\theta_{i}$y las fluctuaciones de densidad$\delta n_{i}$.

Las EoM correspondientes son

$\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)$.

Formalmente, el término$\frac{1}{2}U\delta n_{i}^2$es como una energía cinética, con$U$jugando el papel de una masa inversa. Así por analogía podemos añadir términos de Langevin a la segunda ecuación, como en el ejemplo anterior:

$\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)$.

Lo que me gustaría tener es una expresión similar para el modelo semiclásico de Bose-Hubbard.

Empiezo con el hamiltoniano de Bose-Hubbard semiclásico (con números complejos en lugar de operadores de campo) en representación de estado coherente,

$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}$dónde

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

y luego transformar eso en una representación de momento coordinado, usando

$\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$

con MOE:

$\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}$

¿Cómo se agregan correctamente los términos de Langevin?

Actualización (después del comentario de Ted Pudlik):

Siguiendo la sugerencia de Ted Pudlik, escribo el hamiltoniano semiclásico de Bose-Hubbard en notación de fase de densidad:

$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$

Las ecuaciones de movimiento de Hamilton correspondientes son:

$\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)$

Al igual que en el modelo XY, agrego los siguientes términos a la derivada de densidad (no de fase, como se indicó anteriormente, vea el comentario a continuación):

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

donde elijo$\gamma$de acuerdo con mis necesidades (por ejemplo, menor que la frecuencia propia más pequeña del sistema o sobreamortiguado).

Lo que todavía me preocupa es el$\frac{1}{2\sqrt{n_i}}$término en la derivada de fase: cuando la densidad en el sitio es muy baja en comparación con su vecino, este término diverge. Ese es, por ejemplo, el caso de la nube térmica de un gas ultrafrío en una trampa armónica.

¿Hay alguna forma de transformar la EoM incluyendo los términos de Langevin en$\left(q,p\right)$-representación o$\left(Re,Im\right)$-representación para evitar esto?