<p>We compute the QCD axion contribution to the energy density of dark radiation, parameterized by ∆<i>N</i>_eff, 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 ∆<i>N</i>_eff independently of 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 ∆<i>N</i>_eff 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 ∆<i>N</i>_eff 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.</p>