Does Frisch–Waugh–Lovell theorem work on panel data? We consider the following panel regression model
\begin{align} 
Y_{i} = X_{1,i}\beta_1 + X_{2,i}\beta_2 + \epsilon_{i} , \ i=1,...,N,
\end{align}
where $Y_i := [y_{i1},...,y_{iT}]'$, $X_{k,i}:= [x_{k,i1},...,x_{k,iT}]'$ is a $T \times p_k$ matrix for $k=1,2$, $\epsilon_{i}$ is the $T \times 1$ idiosyncratic error.
It is known that the OLS estimate of the above model is given by
\begin{align*}
\begin{bmatrix}
\widehat{\beta}_1\\
\widehat{\beta}_2
\end{bmatrix}
= (\sum_{i=1}^{N}X_i'X_i)^{-1}\sum_{i=1}^{N}X_i'Y_i,
\end{align*}
where $X_i:= [X_{1,i},X_{2,i}]$.
My question is the following:
whether the OLS estimate $\widehat{\beta}_1$ in the above expression is the same as
\begin{align*}
\widetilde{\beta}_1= (\sum_{i=1}^{N}X_{1,i}'M_{2,i}X_{1,i})^{-1}\sum_{i=1}^{N}X_{1,i}'M_{2,i}Y_i,
\end{align*}
where $M_{2,i} := I_T - X_{2,i}(X_{2,i}'X_{2,i})^{-1}X_{2,i}'$.
I thought $\widetilde{\beta}_1$ is an extension of the Frisch–Waugh–Lovell theorem to the panel data, but I couldn't find the relevant reference.  I would appreciate it if you can share with me your thoughts on this issue, or let me know any related paper/theory/lemma that you might know.

In principle, yes, since FWL is "just" OLS regression algebra and hence does not "care" if the underlying model is cross-sectional, time series, panel etc. Indeed, fixed effects (which you do not seem to have in your model) is one key motivation (at least in my courses...) of the FWL, see e.g. Utility of the Frisch-Waugh theorem
...but not exactly how you describe it, as far as I can see:
What pooled OLS (which is what you state) does is run an OLS regression of the stacked data of $Y$ on stacked data of $X_1$ and $X_2$. What FWL says is that it is equivalent to first collect the residuals from regressing both stacked $Y$ and stacked $X_1$ on stacked $X_2$ and regress these residuals onto each other.
Your notation however suggests to produce these residuals for each unit $i$, which is not equivalent.
Here is a numerical illustration:
library(plm)

# create some artificial panel data
set.seed(324)
m = 8
n = 12

beta1 = -1
beta2 = 2
y = X1 = X2 = r.y = r.X1 = matrix(NA, nrow=m, ncol=n)
for (i in 1:n) {
  X1[,i] = rnorm(m,i)
  X2[,i] = rnorm(m,i)
  y[,i] = X1[,i]*beta1 + X2[,i]*beta2 + rnorm(m,sd=.75)  
  
  # already collect residuals for "your" approach (I suppress a constant in all regressions for simplicity)
  r.y[,i] = resid(lm(y[,i]~X2[,i]-1))
  r.X1[,i] = resid(lm(X1[,i]~X2[,i]-1))
}
stackX1 = as.vector(X1)
stackX2 = as.vector(X2)
stackY = as.vector(y)
unit = rep(1:n,each=m) # first two columns are for plm to understand the panel structure
paneldata = data.frame(unit, rep(1:m,n), stackY , stackX1, stackX2) 

# estimators
coef(plm(stackY~stackX1+stackX2-1, data = paneldata, model = "pooling")) # same as pooled OLS in the line below
coef(lm(stackY~stackX1+stackX2-1))
res.y <- resid(lm(stackY~stackX2-1))
res.x1 <- resid(lm(stackX1~stackX2-1))
coef(lm(res.y~res.x1-1)) # FWL gives same results

# "your" approach
stack.r.X1 = as.vector(r.X1)
stack.r.y = as.vector(r.y)
coef(lm(stack.r.y~stack.r.X1-1)) # slightly different results

Your approach appears to be a mixture of pooled OLS (in that it runs a single regression to produce $\tilde\beta_1$) and what is known as a mean group estimator
\begin{align*}
\frac{1}{N}\sum_{i=1}^{N}\tilde\beta_{1,i}=\frac{1}{N}\sum_{i=1}^{N}(X_{1,i}'M_{2,i}X_{1,i})^{-1}X_{1,i}'M_{2,i}Y_i,
\end{align*}
(in that it runs separate regressions to produce the FWL residuals).
