We present a detailed study of the exclusive radiative decays $Z\to\eta^{(\prime)}\gamma$ employing the QCD factorization approach. We derive a factorization formula for the decay amplitudes valid at leading power in an expansion in $(\Lambda_{QCD}/m_Z)^2$, which includes convolutions of calculable hard-scattering kernels with the leading-twist quark and gluon light-cone distribution amplitudes of the mesons. Large logarithms arising in the evolution from the high scale $m_Z$ down to hadronic scales are resummed using the renormalization group, carefully accounting for the effects of the heavy bottom and charm quarks. Our results for the branching ratios are very sensitive to hadronic input parameters, such as the decay constants and mixing angle characterizing the $\eta-\eta'$ system. Using the most recent estimates of these parameters, we obtain the branching ratios $Br(Z\to\eta\gamma)\sim 1.6\cdot 10^{-10}$ and $Br(Z\to\eta'\gamma)\sim 4.7\cdot 10^{-9}$. A measurement of these processes at a future high-luminosity Z factory could provide interesting information on the gluon distribution amplitude.