Motivated by applications in thermal QCD and cosmology, we elaborate on a general method for computing next-to-leading order spectral functions for composite operators at vanishing spatial momentum, accounting for real, virtual as well as thermal corrections. As an example, we compute these functions (together with the corresponding imaginary-time correlators which can be compared with lattice simulations) for scalar and pseudoscalar densities in pure Yang-Mills theory. Our results may turn out to be helpful in non-perturbative estimates of the corresponding transport coefficients, which are the bulk viscosity in the scalar channel and the rate of anomalous chirality violation in the pseudoscalar channel. We also mention links to cosmology, although the most useful results in that context may come from a future generalization of our methods to other correlators.