Searches for the exclusive decays of Higgs and Z bosons into a vector quarkonium state and a photon are performed in the $\mu ^+\mu ^-\,\gamma $ final state with a proton–proton collision data sample corresponding to an integrated luminosity of 139 fb$^{-1}$ collected at $\sqrt{s}=13$ TeV with the ATLAS detector at the CERN Large Hadron Collider. The observed data are compatible with the expected backgrounds. The 95% confidence-level upper limits on the branching fractions of the Higgs boson decays into $J/\psi \,\gamma $, $\psi (2S)\,\gamma $, and $\Upsilon (1S,2S,3S)\,\gamma $ are found to be $2.0\times 10^{-4}$, $10.5\times 10^{-4}$, and $(2.5,4.2,3.4)\times 10^{-4}$, respectively, assuming Standard Model production of the Higgs boson. The corresponding 95% CL upper limits on the branching fractions of the Z boson decays are $1.2\times 10^{-6}$, $2.4\times 10^{-6}$, and $(1.1,1.3,2.4)\times 10^{-6}$. An observed 95% CL interval of $(-133,175)$ is obtained for the $\kappa _c/\kappa _\gamma $ ratio of Higgs boson coupling modifiers, and a 95% CL interval of $(-37,40)$ is obtained for $\kappa _b/\kappa _\gamma $.