I want to numerically solve a bilevel optimization problem of the form

$$ \min_y f(y, \hat x(y)), \qquad \hat x(y) = \arg\min_x g(x, y) $$

(for simplicity assume that $\min_x g(x, y)$ exists and is unique for each $y$).

Or, to be more precise I want to find a stationary point $(x^*,y^*)$ with $\frac{\partial}{\partial x} g = 0$ and $\frac{\partial}{\partial y} f = 0$.

I am looking for literature on state-of-the-art numerical optimization schemes (gradient based) for this problem.

So far I am only aware of the "approximate solution method" which consists of approximating the optimal solution $\hat x(y)$ of the inner objective with the value returned by applying $T$ steps of some optimization routine $\Phi$ (e.g. gradient descent) $$\hat x(y) \approx x^{(T)}(y), \qquad x^{(t+1)}=\Phi(x^{(t)}, g), \quad x^{(0)}=x_0 $$

Then one can optimize $y$ directly by 'backpropping' through $x^{(T)}$. This seems to work somewhat in practice but convergence to the stationary point is tricky. (in the sense of getting both $\frac{\partial}{\partial y} f$ and $\frac{\partial}{\partial x} g$ small simultaneously). A good reference to a theoretical treatment of this method would also be greatly appreciated.