  We calculate the rate for thermal production of axions via scattering of
quarks and gluons in the primordial quark-gluon plasma. To obtain a finite
result in a gauge-invariant way that is consistent to leading order in the
strong gauge coupling, we use systematic field theoretical methods such as hard
thermal loop resummation and the Braaten-Yuan prescription. The thermally
produced yield, the decoupling temperature, and the density parameter are
computed for axions with a mass below 10 meV. In this regime, with a
Peccei-Quinn scale above 6x10^8 GeV, the associated axion population can still
be relativistic today and can coexist with the axion cold dark matter
condensate.
