Thermal axion production in the early universe goes through several mass thresholds, and the resulting rate may change dramatically across them. Focusing on the KSVZ and DFSZ frameworks for the invisible QCD axion, we perform a systematic analysis of thermal production across thresholds and provide smooth results for the rate. The QCD phase transition is an obstacle for both classes of models. For the hadronic KSVZ axion, we also deal with production at temperatures around the mass of the heavy-colored fermion charged under the Peccei-Quinn symmetry. Within the DFSZ framework, standard model fermions are charged under this symmetry, and additional thresholds are the heavy Higgs bosons masses and the electroweak phase transition. We investigate the cosmological implications with a specific focus on axion dark radiation quantified by an effective number of neutrino species and explore the discovery reach of future CMB-S4 surveys.