(See UPDATE below)
Inspired by @AlexGhorbal 's attempt,
I present a trigonometric derivation of Alex's result with an approximation used only at the last step.
(I was playing around with the equations... and the result fell out.)
The methods shown here are useful for other problems in special relativity,
not just the question in the OP.
I'm going to use rapidity $\theta$, where $v=\tanh\theta$ and thus
the Lorentz factor $\Gamma=\cosh\theta$ and $v\Gamma =\sinh\theta$ and the Doppler factor $k=\exp\theta$.
Using components in the lab frame, energy-momentum conservation yields
\begin{align}
N\cosh\phi &=M\cosh\theta+\epsilon\\
N\sinh\phi &=M\sinh\theta-\epsilon
\end{align}
By addition, we obtain
$$\boxed{
\begin{align}
N\exp\phi &=M\exp\theta,
\end{align}
}
$$
which is essentially the Doppler Effect, where $K=\exp(\theta-\phi)=(N/M)$ is relative-Bondi-Doppler factor.
By subtraction, we obtain
$$
\boxed{\begin{align}
N\exp(-\phi) &=M\exp(-\theta)+2\epsilon,
\end{align}}
$$
[The above boxed equations represent
"energy-momentum conservation in light-cone coordinates" (in the eigenspace)]
We wish to eliminate $N$ (the mass of the final particle resulting from the absorption of the photon)
By multiplication,
$$\boxed{N^2=M^2+2\epsilon M\exp\theta}$$
which is the related to Alex's $\Rightarrow$ equation, but without using any approximation.
(Of course, this is the square-norm of the 4-momentum equation
$\tilde N=\tilde M+\tilde \gamma$.)
Re-write this as
\begin{align}
\left(\frac{N}{M}\right)^2 &=1+\frac{2\epsilon}{M}\exp\theta \\
K^2 &=1+\frac{2\epsilon}{M}\exp\theta \\
(\exp(\theta-\phi))^2 &=1+\frac{2\epsilon}{M}\exp\theta \\
\end{align}
(This is the a-ha moment!) Taking square-roots,
\begin{align}
\exp(\theta-\phi) &=\sqrt{1+\frac{2\epsilon}{M}\exp\theta} \\
\end{align}
and solving for $\exp\phi$, we have
\begin{align}
\exp\phi &=\frac{\exp{\theta}}{\sqrt{1+\frac{2\epsilon}{M}\exp\theta}} \\
\end{align}
This is a physics relationship between the lab frame's Doppler factors of the initial and final particle, the mass $M$ of initial particle and the photon energy $\epsilon$ in the lab frame. That is,
$$
\boxed{\boxed{
\begin{align}
k_N &=\frac{k_M}{\sqrt{1+\frac{2\epsilon}{M}k_M}} \\
\end{align}
}}
$$
which is an exact result.
UPDATE
The double-boxed result is obtained more directly by division of the first pair of boxed equations, rather than multiplication then some additional algebra to eliminate $N$.
\begin{align}
\boxed{\exp(2\phi)
= \frac{M\exp\theta}{M\exp(-\theta)+2\epsilon}}
&= \frac{\exp(2\theta)}{1+\frac{2\epsilon}{M}\exp\theta}\\
\exp\phi
&= \frac{\exp\theta}{\sqrt{1+\frac{2\epsilon}{M}\exp\theta}},
\end{align}
which is the double-boxed equation (when the Doppler-$k$ is introduced)!
In summary,
- express energy-momentum conservation in light-cone coordinates
(by adding and subtracting energy-momentum conservation in rectangular coordinates).
- divide those equations (which relate the square-of-$N$'s-Doppler-factor with $M$'s-Doppler-factor, mass $M$, and the lab-frame energy $\epsilon$ of the beam
-- this is Mermin's aspect-ratio of the "causal diamond of N" (Mermin's "light rectangle")
to obtain the square of the double-boxed-equation
- In the large $\Gamma$ limit, the Doppler-factor $k \approx 2\Gamma$ (as shown below).
- If desired, one could carry out the expansion for $k$ further to get the next order in the approximation in terms of the Lorentz factors $\Gamma$.
Note
\begin{align}
k=\exp\theta &=\cosh\theta+\sinh\theta\\
&= \cosh\theta+\sqrt{\cosh^2\theta-1}\\
&=\Gamma+\sqrt{\Gamma^2-1}
\end{align}
Now, because the particles are assumed to have large Lorentz factors,
we can make the approximation that:
for large $\Gamma$, we have $k\approx 2\Gamma$,
and thus,
\begin{align}
2\Gamma_N &\approx \frac{2\Gamma_M}{\sqrt{1+\frac{2\epsilon}{M}(2\Gamma_M)}} \\
\Gamma_N &\approx \frac{\Gamma_M}{\sqrt{1+\frac{4\epsilon}{M}\Gamma_M}} \\
\end{align}
...in agreement with Alex's final result...
but
using only one approximation, once [for $k_N$ and $k_M$]
at this last step here at the very end.
We see the approximate relation between Lorentz-factors
is based on a similar but exact relation between Doppler factors.
