We study the dynamical scaling of long-range $O(N)$ models after a sudden quench to the critical temperature, using the functional renormalization group approach. We characterize both short-time aging and long-time relaxation as a function of the symmetry index $N$, the interaction range decay exponent $\sigma$ and the dimension $d$. Our results substantially improve on perturbative predictions, as demonstrated by benchmarks against Monte Carlo simulations and the large-$N$ limit. Finally, we demonstrate that long-range systems increase the performance of critical heat engines with respect to a local active medium.