Introduction

The starting point here is this 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 + 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} $$ We want to add a term, or modify existing terms, so that the Zeeman splitting goes back up a bit between $1.5P_{th}$ and $~1.8P_{th}$ as in figs. 4 (a) and (c) in Sawicki et al. 2024. This may be due to a Feshbach resonance where two polaritons form a biexciton. This effect should depend anti-symmetrically on the magnetic field.

The Feshbach resonance

In order to obtain a feshbach resonance term I think there are two methods available. One: go down to the quantum model of the system including the terms relating to the biexciton. A Feshbach resonance should be apparent in that model, and one should then be able to apply a mean field approximation.

Another option comes from the estimation that a term accounting for feshbach resonance should grow significantly when some detuning parameter $\beta$ reaches the resonance value $\beta_0$. A Feshbach resonance may be magnetically or optically tunable (maybe also with some other methods, idk). In a magnetically tunable system where the bound state doesn't decay, the scattering length follows the formula $$ a = a_0 \left(1 - \frac{\Delta}{B-B_0}\right), $$ i.e. eq. 21 in Chin et al. 2010. Here $\Delta$ is the resonance width $\Gamma_0/\delta\mu$, $\Gamma_0$ is the resonance strength and $\delta\mu = \mu_{atoms} - \mu_c$ is the difference between the magnetic moment of the separated atoms and the bound state (in the polaritonic case this would be the difference between the magnetic moments of the polaritons and the biexciton, or maybe the difference between the magnetic moment of the exciton and the biexciton?). Accounting for decay would simply amount to adding an imaginary term $i\gamma/2$ in the denominator.

In eq. 23 in the same text, the scattering length for an optically tuned system with losses is presented: $$ a = a_0\left(1 + \frac{\Gamma_{0}(I)}{h(\nu - \nu_{c} - \delta\nu(I)) + \frac{i\gamma}{2}}\right). $$ Here $I$ is the intensity of the laser that is tuning the system, $\nu$ is its frequency, $\nu_c$ is "the frequency of the unshifted optical transition between the excited bound state and the collisional state of the two atom at $E=0$".

Applying to the GPE model

So basically, if we have a detuning parameter $\beta$, then we want to somehow add a term $$ \frac{\Delta}{\beta - \beta_0 + ia} = \frac{\Delta}{\beta - z} $$ somewhere in the equation, where $a$ is proportional to the biexciton decay rate, and I'll use $z$ for brevity in some of the proceeding formulas where there is a complex number. This could be in the nonlinear term (makes sense, since it represents the polariton self interaction and the feshbach resonance would be caused by increased polariton self-interactions, i.e. greater scattering lengths), by modifying the factor $g$ that helps define $\alpha_\pm$ and $G_\pm$ so that it is now $$ g = g_0\left(1 - \frac{\Delta}{\beta - z}\right) $$ This would introduce an imaginary component to $\alpha_\pm$ and $G_\pm$, which may be a bit unexpected, but not necessarily unphysical, it would just be a bit of extra decay rate.

Since the effect is also seen at equally strong magnetic fields with the sign flipped, we should actually have a split $g_\pm$, where the former equation was $g_+$ and $g_-$ is given by $$ g_- = g_0 \left(1 - \frac{\Delta}{\beta + \bar{z}}\right). $$

What's somewhat concerning is if the detuning parameter is $B = B_{ext} + B_{int}$. We would need to include $B_{int}$ since the non-monotonicity observed in the sample arises with different levels of pumping, so there would be no such effect if only $B_{ext}$ was used. The problem is that we now have a circular dependency: $$ B_{int} = \alpha_+\rho_+ - \alpha_-\rho_- + G_+n_+ - G_-n_- $$ $$ \alpha_\pm = \alpha_{\pm,0} \left(1 - \frac{\Delta}{B_{int} - z}\right) $$ $$ G_\pm = G_{\pm,0}\left(1 - \frac{\Delta}{B_{int} - z}\right) $$

where $\rho_\pm = |\psi_\pm|$and it is not clear how it should be resolved. Possibly, just by taking some starting values for $\alpha_\pm$, $G_\pm$ and $B_{int}$ and iterating, one might arrive close to a fixed point. We do also in fact have four equations in four unknowns here, however these are all non-linear equations.

One way to greatly simplify this is to make $B_{int}$ only depend on the unmodified parameters $\alpha_{\pm,0}$ and $G_{\pm,0}$. It should be borne in mind that $(\alpha,G)_{\pm,0}$ are not the same as $(\alpha,G)_\pm$ measured at 0 magnetic field, but rather the limits as the magnetic field goes infinitely far away from the resonances. This may, however, be basically the same as the value at $B_{ext}=0$ depending on the width of the resonance.

One problem I see with having the resonance depend primarily on the magnetic field is that the effect doesn't seem to be observed by just changing the applied magnetic field, as one might expect. Given the equations presented above it shouldn't matter whether the magnetic field is generated externally or by the polaritons, and yet there's not really any non-monotonicity going from 0 to 5 T applied magnetic field aside from just a few data points that don't show a consistent pattern.

Another possibility is that the tuning really depends on $\rho_\pm$. This could make sense as the effect depends on the energy of the +/- spin state going up/down compared to the -/+ state. This is consistent with the population of the +/- state being greater at either magnetic field (iirc the + population was greater at more positive magnetic fields, which I think would actually be contrary to equations 12/13 in Sawicki et al., I would need to make the plot again to be sure, and it might have been the wrong plot anyway), then as the pumping strength increases the effect dissipates, and indeed at greater strengths the population keeps increasing, which would move the state away from the resonance.

We would then be looking at a term like $$ g_\pm \propto 1 - \frac{\Delta}{\rho_\pm - z}. $$ There is still the question of whether one should just use the complex value, or use the absolute value of the term. An alternative would be a term like $$ g_\pm \propto 1 - \frac{\Delta}{(\rho_\pm - \rho_{\pm,0})^2 + a} $$ which has a similar behaviour where the value increases at $\rho_{\pm,0}$, but never goes to infinity while being purely real. However, I don't really see an analogue to such a term in the Feshbach resonance literature and I'm not sure how one would justify it physically.