I am studying a problem arising in physics, and I managed to simplify it to a differential system (initial value problem) of the form: $$ \begin{cases} \dot{x} = \epsilon f_1(x,y,t) + \epsilon^2 f_2(x,y,t) + \mathcal O\left(\epsilon^3\right),\\ \epsilon \dot{y} = g_0(x,y,t) + \epsilon g_1(x,y,t) + \epsilon^2 g_2(x,y,t) + \mathcal O\left(\epsilon^3\right) \end{cases} $$ where $0<\epsilon\ll 1$ is a small parameter, $x$ and $y$ are real-valued functions, the $f_i$ and $g_i$ are $\mathcal O(1)$ and $2\pi$-periodic in $t$. I have the actual expressions for these functions, but I don't think this is relevant here, they are just very nonlinear in all the variables.
I want to find an asymptotic solution of this system valid up to $\mathcal O\left(\epsilon^2\right)$, under the initial condition $x(0)=1$, $y(0)=0$. I actually am interested in the long-time behaviour of $\langle x \rangle(t) = \frac{1}{2\pi}\int_t^{t+2\pi} x(s)~{\rm d}s$, and I have reasons to believe this requires an asymptotic analysis up to $\epsilon^2$.
After some research, I found that this is a mix of 1) a singular system of ODE, what is sometimes referred to as a Tikhonov problem in the literature (see e.g. this question); and 2) a problem that could be suited for the averaging method (Wikipedia).
The only solution for the degenerate system for $y$ at $\mathcal O(1)$, $g_0(x,y,t)=0$, is $y(t)=0$. The associated equation for $x$, $\dot{x} = \epsilon f_1(x,0,t) + \epsilon^2 f_2(x,0,t)$, remains tricky to deal with. But I could easily use the averaging method to $\mathcal O\left(\epsilon^2 \right)$ to find an approximation of $\langle x \rangle$. Of course, this is inconsistent, and I must deal with the system properly. This just illustrates that the averaging method -- which essentially uses a separation of timescales between the "fast" time $t$ and a "slow" time associated with $\dot x$ -- seems relevant.
I would like to know if there is already a framework or specific methods that discuss asymptotic solutions of such time-periodic Tikhonov problems with a slow dynamic for the first variable. I feel like this may have specific properties that could well have been studied already; though I have not been able to find any reference. However, I am venturing very far from my comfort zone and field (experimental physics), and I may have overlooked important literature or keywords.