"Classical" adiabatic approximation: please help me understand better In this1 paper they use an adiabatic approximation to reduce two differential equations to one. Could you please recommend some alternative reading for this (semi) classical adiabatic approximation, or offer some additional explanation. I can only find discussions in the context of quantum Hamiltonians which don't help me understand this.
Here I present the equations stripped back to the essentials needed for my question. The evolution of an exciton reservoir with density $N(t)$, and polariton condensate $\psi(t)$, are described by the equations
$$\begin{aligned}
\frac{\text{d}}{\text{d}t} \psi &= \left( R N - \gamma \right) \psi, \\
\frac{\text{d}}{\text{d}t} N &= P - \left( \Gamma + R |\psi|^2\right) N,
\end{aligned}$$
where all the other parameters are time-independent. We then assume that $N$ adiabatically follows $|\psi|^2$. In other words $|\psi|^2$ changes on a much longer time-scale than $N$; $N$ is able to quickly adjust to the condensate density $|\psi|^2$. In this case we can treat $|\psi|^2$ as a constant in the equation for $\text{d}N/\text{d}t$, whose solution is
$$
N = \frac{P}{\Gamma + R |\psi|^2} + C \text{e}^{-(\Gamma + R |\psi|^2)t}.
$$
This can simply be substituted back in (ignoring the second term in the equation for some reason I don't understand) to give just one equation
$$
\frac{\text{d}}{\text{d}t} \psi = \left( \frac{R P}{\Gamma + R |\psi|^2} - \gamma \right) \psi.
$$
For this approximation to hold, the terms in the equation for $\psi$ (i.e. $\gamma$) need to be much smaller than the inverse time-scale $\tau^{-1} = (\Gamma + R|\psi|^2)$.
Whilst this all seems somewhat reasonable to me. I am lacking an understanding of exactly why this works and why this approximation is ok.

1 N. Bobrovska and M. Matuszewski, "Adiabatic approximation and fluctuations in exciton-polariton condensates", Phys. Rev. B 92, 035311 (2015), arXiv:1505.06663.
 A: Let us consider a simpler equation:
$$\dot{N}=-\Gamma N + f(t),$$
where $f(t)$ is an external force (i.e., independent on $N$). We can solve this equation by the method of variation of constants, i.e., by solving first
$\dot{N}=-\Gamma N \rightarrow N(t)=N_1e^{-\Gamma t},$
and then treating $N_1$ as a time-dependent function:
$N(t)=N_1(t)e^{-\Gamma t}$. When plugging this into the original equation we obtain for $N_1$ equation $\dot{N}_1=f(t)e^{\Gamma t}$, whereas the full solution is obtained as
$$N(t)=N(0)e^{-\Gamma t}+\int_0^t dt_1f(t_1)e^{-\Gamma (t-t_1)} = N(0)e^{-\Gamma t}+\int_0^t d\tau f(t-\tau)e^{-\Gamma \tau}.$$
We can now expand the function in the last integral in Taylor series:
$$
\int_0^t d\tau_1f(t-\tau)e^{-\Gamma \tau}=\sum_{n=0}^{+\infty}\frac{(-1)^nf^{(n)}(t)}{n!}\int_0^td\tau\tau^ne^{-\Gamma \tau}
=\sum_{n=0}^{+\infty}\frac{(-1)^nf^{(n)}(t)}{n!\Gamma^{n+1}}\int_0^{\Gamma t}du u^n e^{-u}.$$
The last integral can be evaluated exactly, but in can use the limit $\Gamma t\rightarrow+\infty$, so that we have $$\int_0^{\Gamma t}du u^n e^{-u}=n!\Rightarrow \frac{1}{n!}\int_0^{\Gamma t}du u^n e^{-u}\leq 1.$$
Thus our series above can be approximated by
$$
\int_0^t d\tau_1f(t-\tau)e^{-\Gamma \tau}\approx\sum_{n=0}^{+\infty}\frac{(-1)^nf^{(n)}(t)}{\Gamma^{n+1}}$$
The condition of $f(t)$ changing slowly in comparison with $e^{-\Gamma t}$ is the condition that
$$\left|\frac{f^{(n)}(t)}{\Gamma^{n+1}}\right|\ll 1,$$
and we neglect all the terms except the first one ($f^{(0)}(t)=f(t)$).
The analysis for the system in the OP proceeds along the same lines, but mathematically is a bit more involved - it could be  agood homework problem ;)
