  We investigate the masses and mixing of $\eta$, $\eta'$ mesons in the
framework of twisted mass lattice QCD with $N_f=2+1+1$ dynamical quark flavors.
For the first time we perform a controlled chiral and continuum extrapolation
to obtain physical results. To this end, we have analyzed 17 gauge ensembles
with pion masses ranging from $220\,\mathrm{MeV}$ to $500\,\mathrm{MeV}$ and at
three values of the lattice spacing. This calculation became feasible through
the application of a powerful variance reduction technique for the computation
of quark-disconnected diagrams together with several improvements in the final
analysis, including a correction of the relevant correlation functions for
residual topological finite volume effects. We obtain physical results for the
masses $M_{\eta}=557(11)_\mathrm{stat}(03)_{\chi\mathrm{PT}}\,\mathrm{MeV}$ and
$M_{\eta'}=911(64)_\mathrm{stat}(03)_{\chi\mathrm{PT}}\,\mathrm{MeV}$, as well
as the mixing angle in the quark flavor basis
$\phi=38.8(2.2)_\mathrm{stat}(2.4)_{\chi\mathrm{PT}}^\circ$, in excellent
agreement with other results from phenomenology. Similarly, for the physical
values of the decay constant parameters we find
$f_l=125(5)_\mathrm{stat}(6)_{\chi\mathrm{PT}}\,\mathrm{MeV}$ and
$f_s=178(4)_\mathrm{stat}(1)_{\chi\mathrm{PT}}\,\mathrm{MeV}$. Finally, we
present a test of the Veneziano Witten formula using our lattice data.
