  We determine the masses, the singlet and octet decay constants as well as the
anomalous matrix elements of the $\eta$ and $\eta^\prime$ mesons in $N_f=2+1$
<EMAIL_ADDRESS>The results are obtained using twenty-one CLS ensembles of
non-perturbatively improved Wilson fermions that span four lattice spacings
ranging from $a\approx 0.086\,$fm down to $a\approx 0.050\,$fm. The pion masses
vary from $M_{\pi}=420\,$MeV to $126\,$MeV and the spatial lattice extents
$L_s$ are such that $L_sM_\pi\gtrsim 4$, avoiding significant finite volume
effects. The quark mass dependence of the data is tightly constrained by
employing two trajectories in the quark mass plane, enabling a thorough
investigation of U($3$) large-$N_c$ chiral perturbation theory (ChPT). The
continuum limit extrapolated data turn out to be reasonably well described by
the next-to-leading order ChPT parametrization and the respective low energy
constants are determined. The data are shown to be consistent with the singlet
axial Ward identity and, for the first time, also the matrix elements with the
topological charge density are computed. We also derive the corresponding
next-to-leading order large-$N_{c}$ ChPT formulae. We find $F^8 =
115.0(2.8)~\text{MeV}$, $\theta_{8} = -25.8(2.3)^{\circ}$, $\theta_0 =
-8.1(1.8)^{\circ}$ and, in the $\overline{\mathrm{MS}}$ scheme for $N_f=3$,
$F^{0}(\mu = 2\,\mathrm{GeV}) = 100.1(3.0)~\text{MeV}$, where the decay
constants read $F^8_\eta=F^8\cos \theta_8$, $F^8_{\eta^\prime}=F^8\sin
\theta_8$, $F^0_\eta=-F^0\sin \theta_0$ and $F^0_{\eta^\prime}=F^0\cos
\theta_0$. For the gluonic matrix elements, we obtain $a_{\eta}(\mu =
2\,\mathrm{GeV}) = 0.0170(10)\,\mathrm{GeV}^{3}$ and $a_{\eta^{\prime}}(\mu =
2\,\mathrm{GeV}) = 0.0381(84)\,\mathrm{GeV}^{3}$, where statistical and all
systematic errors are added in quadrature.
