  We present a calculation of the hindered M$1$ $\Upsilon(2S) \to \eta_b(1S)
\gamma$ decay rate using lattice non-relativistic QCD. The calculation includes
spin-dependent relativistic corrections to the NRQCD action through
$\mathcal{O}(v^6)$ in the quark's relative velocity, relativistic corrections
to the leading order current which mediates the transition through the quark's
magnetic moment, radiative corrections to the leading spin-magnetic coupling
and for the first time a full error budget. We also use gluon field ensembles
at multiple lattice spacing values, all of which include $u$, $d$, $s$ and $c$
quark vacuum polarisation. Our result for the branching fraction is
$\mathcal{B}(\Upsilon(2S)\to\eta_b(1S)\gamma) = 5.4(1.8)\times 10^{-4} $, which
agrees with the current experimental value.
