  Freeze-in via the axion-photon coupling, $g_{\phi\gamma}$, can produce axions
in the early Universe. At low reheating temperatures close to the minimum
allowed value $T_{\rm reh}\approx T_{\rm BBN}\approx 10\,{\rm MeV}$, the
abundance peaks for axion masses $m_\phi\approx T_{\rm reh}$. Such heavy axions
are unstable and subsequently decay, leading to strong constraints on
$g_{\phi\gamma}$ 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.
