Integration of nonlinear PIDE via spectral methods

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:

$$g_{k,q_1,q_2}=g\, X^*(k)\, X^*(q_1+q_2-k)\, X(q_1)\, X(q_2)$$

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

$$\psi(k)=\int dr e^{-ikr}\tilde{\psi}(r)$$

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

• Generally, if you're working in one dimension and periodic BCs are acceptable, then Fourier pseudospectral methods are a good option. That would mean that the integrand is evaluated in physical space. – David Ketcheson May 30 '14 at 4:23
• I assume you are using $\star$ to denote the complex conjugate? Usually it's $*$. – David Ketcheson May 30 '14 at 4:24
• Periodic BCs are acceptable but unfortunately most of the time I am working in 2D. Are there any tricks to evaluate the integrand faster? – Andrei May 30 '14 at 11:10
• I think you can discretise the equation in $k$-space and solve the resulting set of coupled ODEs. I don't see what's wrong with that approach, which you mentioned yourself. What is the domain of $q_1,q_2$? You can evaluate the sum approximately, by summing over the discretised $q_{1,2}$'s. This is slow for sure, but this is a $(2+1)$-dimensional PIDE, after all. I can't think of a better way. – Kirill May 30 '14 at 16:05
• If you know that $\tilde{\psi}$ decays rapidly, then you could probably use something like fast multipole to rapidly evaluate the integral. But that's a major undertaking and will only be useful for a fairly large-scale discretization. – David Ketcheson May 30 '14 at 18:10