Precise predictions for the QCD axion contribution to dark radiation with full phase-space evolution

We compute the QCD axion contribution to the energy density of dark radiation, parameterized by $\mathrm{\Delta }{N}_{\mathrm{e}\mathrm{f}\mathrm{f}}$, by solving Boltzmann equations for the momentum distribution functions including the effects of quantum statistics for all particles involved in the axion production processes. This approach gives precise prediction for $\mathrm{\Delta }{N}_{\mathrm{e}\mathrm{f}\mathrm{f}}$ independently on whether axions are produced via freeze-out or freeze-in. We focus on axions produced via flavor-conserving and flavor-violating interactions with leptons. Our precise predictions for $\mathrm{\Delta }{N}_{\mathrm{e}\mathrm{f}\mathrm{f}}$ can differ from those assuming thermal shape for the momentum distribution functions, as commonly done in the literature, by more than the experimental precision of future Cosmic Microwave Background (CMB) observations. Current lower limits on the axion couplings from Planck constraints on $\mathrm{\Delta }{N}_{\mathrm{e}\mathrm{f}\mathrm{f}}$ are also affected by our precise computation which, in particular, results in a strongly relaxed bound on flavor-violating axion couplings to tau lepton and muon or electron.