  For a discrete symmetry that is anomalous under QCD, the domain walls
produced in the early universe from its spontaneous breaking can naturally
annihilate due to QCD instanton effects. The gravitational waves generated from
wall annihilation have their amplitude and frequency determined by both the
discrete symmetry breaking scale and the QCD scale. The evidence of stochastic
gravitational waves at nanohertz observed by pulsar timing array experiments
suggests that the discrete-symmetry-breaking scale is around 100 TeV, assuming
the domain-wall explanation. The annihilation temperature is about 100 MeV,
which could naturally be below the QCD phase transition temperature. We point
out that the QCD phase transition within some domains with an effective large
QCD $\theta$ angle could be a first-order one. To derive the phase diagram in
$\theta$ and temperature, we adopt a phenomenological linear sigma model with
three quark flavors. The domain-wall explanation for the NANOGrav, EPTA, PPTA
and CPTA results hints at a first-order QCD phase transition, which predicts
additional gravitational waves at higher frequencies. If the initial formation
of domain walls is also a first-order process, this class of domain-wall models
predicts an interesting gravitational wave spectroscopy with frequencies
spanning more than ten orders of magnitude, from nanohertz to 100 Hz.
