  The search for axion-like particles $X$ in muon decays is an excellent
opportunity for the MEG II and Mu3e experiments to extend their horizons beyond
$\mu^+ \to e^+ \gamma$ and $\mu^+ \to e^+ e^- e^+$. A suitable process for both
experiments is the two-body decay $\mu^+ \to e^+ X$, whose only signature is a
monochromatic peak close to the kinematic endpoint of the positron energy
spectrum of the $\mu^+ \to e^+ \nu_e \bar\nu_\mu$ background. The hunt for such
an elusive signal in a vast amount of irreducible background requires extremely
accurate theoretical predictions to be implemented in a Monte Carlo event
generator. This work presents a new state-of-the-art computation of $\mu^+ \to
e^+ \nu_e \bar\nu_\mu$ for polarised muons, accomplished with the McMule
framework. The calculation includes next-to-next-leading order QED corrections
and logarithmically enhanced terms at even higher orders. The results are also
used to estimate the sensitivity of both experiments on the branching ratio of
$\mu^+ \to e^+ X$, in order to evaluate the impact of the theoretical error.
