Frisch–Waugh–Lovell (FWL) Involving Interaction Terms in OLS Suppose that I have two structural equations:
$$ X = \alpha_1W + \epsilon_1$$
$$ Y = \beta_1X + \beta_2W + \beta_3XW + \epsilon_2$$
In other words, $W \sim N(0,1)$ is a confounder and effect modifier. If want to recover $\beta_1$ using the FWL theorem, by the standard definition, I would perhaps compute
$$\frac{X^TMY}{X^TMX} \overset{p}{\to} \beta_1$$
where $M$ is a residual-making matrix. The issue comes in with the interaction term, which I am unsure how to orthogonalize because simply computing $X - E[X|W] = X - \alpha_1W$ and $Y - E[Y|W] = Y - \beta_2W$ seems to be inaccurate if $\beta_3 \ne 0$.
So the questions are 1) can I use FWL in this case? and 2) if so, what would $M$ look like?

I would say as "usual", i.e. project both $Y$ and $X$ on all other regressors in the first stage. In this case, these are the constant, $W$ and the interaction term. Collect the residuals and regress these onto each other.
Hence, denoting by $Z:=(i\quad W\quad X\cdot W)$ the matrix on the regressors except $X$, we have
$$
M=I-Z(Z'Z)^{-1}Z'.
$$
Illustration (since FWL is "just" an algebraical device, it does not matter whether $W$ is or is not a mediator in the underlying DGP):
n <- 20
Y <- rnorm(n)
X <- rnorm(n)
W <- rnorm(n)

summary(lm(Y~X*W))

first.stage.Y <- lm(Y~W+I(X*W))
first.stage.X <- lm(X~W+I(X*W))

summary(lm(resid(first.stage.Y)~resid(first.stage.X)-1))

