Table of Contents
The model
The model is the following discretised generalised coupled GPE:
$$ \frac{d\psi_\pm}{dt} = \left[\frac{1}{2}(Rn_\pm - \gamma) - \frac{i}{\hbar}\left(E_\pm + \alpha_\pm|\psi_\pm|^2 + \beta_\pm|\psi_\mp|^2 + K_\pm\left(n_\mp + \frac{\mathcal{P}_\mp}{W}\right) + G_\pm\left(n_\pm + \frac{\mathcal{P}_\pm}{W}\right)\right)\right]\psi_\pm - \frac{iE_{XY}}{2\hbar}\psi_\mp $$ $$ \frac{dn_\pm}{dt} = -(\Gamma + \Gamma_{s,\pm} + R|\psi_{\pm}|^2)n_\pm + \Gamma_{s,\mp}n_\mp + \mathcal{P}_{\pm} $$
This is somewhat different from the model use in a previous post in that opposite spin interactions between the polaritons are accounted for by the $\beta_\pm$ and $K_\pm$ terms. These are generally neglected, but for this specific system the behaviour we're looking for should be caused by precisely this interaction.
The procedure here is to initialize $\psi$ to small random values, run the simulation until a steady state solution is reached and find the difference in energy levels between $\psi_+$ and $\psi_-$. I do this by taking the Fourier transform of the product $\psi_+\psi_-^*$ and