Freeze-in via the axion-photon coupling, gϕγ, can produce axions in the early Universe. At low reheating temperatures close to the minimum allowed value Treh ≈ TBBN ≈ 10 MeV, the abundance peaks for axion masses mϕ ≈ Treh. Such heavy axions are unstable and subsequently decay, leading to strong constraints on gϕγ from astrophysics and cosmology. In this work, we revisit the computation of the freeze-in abundance and clarify important issues. We begin with a complete computation of the collision terms for the Primakoff process, electron-positron annihilation, and photon-to-axion (inverse-)decay, while approximately taking into account plasma screening and threshold effects. We then solve the Boltzmann equation for the full axion distribution function. We confirm previous results about the importance of both processes to the effective “relic abundance” (defined as density prior to decay), and provide useful fitting formulae to estimate the freeze-in abundance from the equilibrium interaction rate. For the distribution function, we find an out-of-equilibrium population of axions and introduce an effective temperature for them. We follow the evolution right up until decay, and find that the average axion kinetic energy is larger than a thermal relic by between 20% and 80%, which may have implications for limits on decaying axions from X-ray spectra. We extend our study to a two-axion system with quartic cross-coupling, and find that for typical/expected couplings, freeze-in of a second axion flavour by annihilations leads to a negligibly small contribution to the relic density.