At the mean-field level, the dynamics of a polariton condensate can be described by a type of nonlinear Schrodinger equation (Gross-Pitaevskii-type), for a classical (complex-number) wavefunction $\psi$. Its form in momentum space reads (in 2+1 dimensions):
\begin{multline} i \frac{d}{dt}\psi(k) =\left[\epsilon(k) -i\frac{\gamma(k)}{2}\right] \psi(k) +F_{p}(k)\,\, e^{-i\omega_{p}t} \\ + \sum_{q_1,q_2} g_{k,q_1,q_2}\, \psi^*(q_1+q_2-k) \, \psi(q_1)\, \psi(q_2). \end{multline}
The function $\epsilon(k)$ is the dispersion of the particles (polaritons). The polaritons are a non-equilibrium system, due to their finite lifetime (damping rate $\gamma$). Therefore, they need continuous pumping with amplitude $F_p$ at energy $\omega_p$. Finally, there exists a momentum-dependent nonlinear interaction $g_{k,q_1,q_2}$ that depends of the so-called Hopfield coefficients $X$ (simple functions of momentum) as:
\begin{equation} g_{k,q_1,q_2}=g\, X^*(k)\, X^*(q_1+q_2-k)\, X(q_1)\, X(q_2) \end{equation}
Equivalently, the real-space correspondent of this equation is:
\begin{multline} i\frac{d}{dt}\tilde{\psi}(r)=\left[\epsilon(i\partial_{r})-i\frac{\gamma(i\partial_{r})}{2}\right]\tilde{\psi}(r)+F_{p}e^{i(k_{p}r-\omega_{p}t)}\\ +\iiint dr_{123}\,\tilde{g}(r_1-r,r-r_2,r-r_3)\tilde{\psi}^*\!(r_1)\tilde{\psi}(r_2)\tilde{\psi}(r_3) \end{multline}
with
\begin{multline} \tilde{g}\left(r_{1}-r,r-r_{2},r-r_{3}\right)=g\int dq_{3}X(q_{3})e^{iq_{3}\left(r-r_{3}\right)}\int dq_{1}X^*(q_{1})e^{iq_{1}\left(r_{1}-r\right)}\\ \times \int dq_{2}\, X^*(q_{2}+q_{3}-q_{1})\, X(q_{2})\, e^{iq_{2}\left(r-r_{2}\right)} \end{multline}
where i have used that
\begin{equation} \psi(k)=\int dr e^{-ikr}\tilde{\psi}(r) \end{equation}
My question is how to best solve this equation numerically. Boundary conditions are not important, one can use PBC if that makes it easier.
In fact, I was thinking that it would maybe be easier to solve the equation in k-space, wouldn't that be equivalent to using a spectral method? I am not sure how to do it in practice though, and especially how to take care of the nonlinear interaction term containing the 2 sums over $q_1$ and $q_2$.
