  Gravitational waves (GWs) at ultra-low frequencies (${\lesssim
100\,\mathrm{nHz}}$) are key to understanding the assembly and evolution of
astrophysical black hole (BH) binaries with masses $\sim
10^{6}-10^{9}\,M_\odot$ at low redshifts. These GWs also offer a unique window
into a wide variety of cosmological processes. Pulsar timing arrays (PTAs) are
beginning to measure this stochastic signal at $\sim 1-100\,\mathrm{nHz}$ and
the combination of data from several arrays is expected to confirm a detection
in the next few years. The dominant physical processes generating gravitational
radiation at $\mathrm{nHz}$ frequencies are still uncertain. PTA observations
alone are currently unable to distinguish a binary BH astrophysical foreground
from a cosmological background due to, say, a first order phase transition at a
temperature $\sim 1-100\,\mathrm{MeV}$ in a weakly-interacting dark sector.
This letter explores the extent to which incorporating integrated bounds on the
ultra-low frequency GW spectrum from any combination of cosmic microwave
background, big bang nucleosynethesis or astrometric observations can help to
break this degeneracy.
