Constraining Postinflationary Axions with Pulsar Timing Arrays

Models that produce Axion-Like-Particles (ALP) after cosmological inflation due to spontaneous $U(1)$ symmetry breaking also produce cosmic string networks. Those axionic strings lose energy through gravitational wave emission during the whole cosmological history, generating a stochastic background of gravitational waves that spans many decades in frequency. We can therefore constrain the axion decay constant and axion mass from limits on the gravitational wave spectrum and compatibility with dark matter abundance as well as dark radiation. We derive such limits from analyzing the most recent NANOGrav data from Pulsar Timing Arrays (PTA). The limits are similar to the $N_{\rm eff}$ bounds on dark radiation for ALP masses $m_a \lesssim 10^{-22}$ eV. On the other hand, for heavy ALPs with $m_a\gtrsim 0.1$ GeV and $N_{\rm DW}\neq 1$, new regions of parameter space can be probed by PTA data due to the dominant Domain-Wall contribution to the gravitational wave background.


I. INTRODUCTION
Pulsar Timing Arrays (PTA) offer a new window to observe the Universe through gravitational waves (GW) in the nano-Hertz frequency range [1][2][3][4][5][6].A potential source of GWs at these frequencies is a population of supermassive black-hole binaries (SMBHBs) in the local universe [6,7].Besides, cosmic strings, which may have been produced in the early Universe during a spontaneous U (1) symmetry-breaking event [8][9][10][11], generate a stochastic gravitational-wave background (SGWB) down to these low frequencies as part of a vast spectrum spanning many decades in frequency; see [12,13] for recent reviews.In fact, given the very wide and nearly scaleinvariant GW spectrum from cosmic strings, the PTA limits are very relevant to anticipate the prospects for probing a cosmic-string GW signal at LISA [14] or Einstein Telescope [15].Cosmic strings can either be local or global depending on whether the spontaneously broken symmetry is a gauge or global U (1).Models of local strings have been confronted to PTA data in [16][17][18][19], and most recently to the 15-year NANOGrav (NG15) data in [5] and the EPTA data release 2 [6,20].
This paper focuses instead on GW from global strings [12,[21][22][23][24][25], which were not analyzed in [5].Many Standard-Model extensions feature such additional global U (1) symmetry that gets spontaneously broken by the vacuum expectation value of a complex scalar field, thus delivering a Nambu-Goldstone Boson.A famous example is the Peccei-Quinn U (1) symmetry advocated to solve the strong CP problem and its associated axion particle [26][27][28][29].Because the U (1) symmetry gets also broken explicitly at later times, the axion acquires a mass.At that moment, domain walls can also populate the Universe [30].
This paper considers this broad class of models of socalled Axion-Like-Particles (ALPs) with mass m a and decay constant f a , corresponding to the energy scale of spontaneous symmetry breaking.If the cosmic-string and domain-wall formations happen before inflation, those are diluted away.On the other hand, if the U (1) is broken at the end or after inflation (in this case, the ALP is dubbed postinflationary), cosmic strings give rise to a population of loops that generate a SGWB throughout the cosmic history.At the same time, they also generate axion particles [31][32][33][34][35][36][37], while domain walls bring an additional contribution to the GW spectrum [38][39][40][41][42][43][44][45].
We aim to use the most recent limits on the SGWB from NG15 data set to derive independent bounds on the parameter space of postinflationary ALPs.Given that a GW signal has been observed [1], any further improved sensitivity from future PTA observatories will not enable pushing down the constraints.Therefore, the PTA constraints presented in this paper on the axion mass and decay constant are not expected to change by more than a factor of a few from future PTA experiments.On the other hand, future GW experiments operating in other frequency ranges will serve as complementary probes to PTA.
Our approach is the following.We analyze the recent NG15 data via the code PTArcade [46,47], first considering the two SGWB from global cosmic strings and domain walls without the astrophysical background.We compare the interpretation of data in terms of SMBHBs and of global cosmic strings and domain walls by calculating the Bayes Factor (BF).Next, we set constraints on the new physics contribution, leading to a SGWB that is too strong and conflicts with the data.The results of the best fit and the constraints on the SGWB from domain walls have been presented in the recent analysis with NG15 data by the NANOGrav collaboration [5].Regard-ing the analysis of previous data release, Refs.[48][49][50] fitted the domain-wall and/or global-string signal to the PPTA second data release (DR2) or IPTA DR2 and/or NANOGrav 12.5-year data, however did not derive the exclusion region.(See Sec.III for more details on the best fit and constraint.)We further translate these bounds into constraints in the ALP parameter space.In addition, this work presents a similar analysis (determining best fits and setting constraints) for global strings for the first time with NG15.
Sec. II of this paper summarizes the postinflationary axion scenarios and their corresponding GW signals, separated into two cases: either cosmic-string or domainwall SGWB dominates.In Sec.III, we confront these cases with the NG15 data and derive, for each case, the constraints on axion parameter space, illustrated in Fig. 3.We conclude in Sec.IV.Supplemental material contains miscellaneous details, such as the priors for analysis, the case assuming no astrophysical background, and the result for the global strings in the m a → 0 limit.

II. POSTINFLATIONARY AXION AND ITS GRAVITATIONAL WAVES
The ALP can be defined as the angular mode θ of a complex scalar field Φ ≡ ϕ exp(iθ) with ϕ the radial partner.It has the Lagrangian density, symmetry restoration and trapping Φ → 0 at early times.The potential has three terms: where f a is the vacuum expectation value of the field, m a ≡ m a (T ) is the axion mass as a function of the Universe's temperature T , N DW is the number of domain walls, and V bias is some further explicit U (1) breaking term.The first term is responsible for U (1) spontaneous breaking, while the second and third terms explicitly break the U (1).These three terms are ranked according to their associated energy scales (large to small) corresponding to their sequences in defect formations: from cosmic strings to domain walls and then their decays.During inflation, the complex scalar field is driven to the minimum of the potential V (Φ) if V c ≪ V (Φ).Quantum fluctuations along the axion direction due to the de Sitter temperature O(H inf ) can generate a positive quadratic term in the potential and restore the U (1) symmetry, which gets eventually broken at the end of inflation, leading to cosmic strings if H inf /(2πf a ) ≳ 1 [51][52][53].However, the current CMB bound [54] on the inflationary scale H inf < 6.1×10 13 GeV implies that f a is too small to generate an observable cosmic-string SGWB.Still, there are several other ways in which U (1) can get broken after inflation even for large f a : i) A large and positive effective ϕ-mass can be generated by coupling ϕ to the inflaton χ (e.g., L ⊃ χ 2 ϕ 2 ) which, for large χ, traps ϕ → 0 during inflation 1 .ii) ϕ could couple to a thermal (SM or secluded) plasma of temperature T that would generate a large thermal V c correction, restoring the U (1) 2 .iii) Lastly, non-perturbative processes, such as preheating, could also lead to U (1) restoration after inflation [57][58][59][60][61].
When V c drops, the first term of V (Φ) breaks spontaneously the U (1) symmetry at energy scale f a , leading to the network formation of line-like defects or cosmic strings with tension µ = πf 2 a log λ 1/2 f a /H [11].As U (1) symmetry is approximately conserved when the axion mass is negligible, the cosmic strings survive for long and evolve into the scaling regime by chopping-off loops [62][63][64][65][66][67][68][69][70][71][72][73][74][75].Loops are continuously produced and emit GW and axion particles throughout cosmic history.The resulting GW signal corresponds to a SGWB entirely characterized by its frequency power spectrum.The latter is commonly expressed as the GW fraction of the total energy density of the Universe h 2 Ω GW (f GW ).
A loop population produced at temperature T quickly decays into GW of frequency [12], where α ∼ O(0.1) is the typical loop size in units of the Hubble horizon 1/H.If the network of cosmic strings is stable until late times, i.e., in the limit m a → 0, its SGWB is characterized by [12,76], where G(T ) ≡ [g * (T )/g * (T 0 )][g * s (T 0 )/g * s (T )] 4/3 with T 0 the temperature of the Universe today.The logarithmic correction is defined by 1 As the inflaton field value relates to the Hubble parameter, this mass is called Hubble-induced mass. 2 For example, the KSVZ-type of interaction couples ϕ to a fermion ψ charged under some gauge symmetry with Aµ: L ⊃ yϕ ψψ + h.c.+ g ψγ µ ψAµ, that can generate thermal corrections: Vc = y 2 T 2 ϕ 2 for yϕ < T and Vc = g 4 T 4 ln y 2 ϕ 2 /T 2 for yϕ ≳ T [55,56].When Vc > λf 4 a , the ϕ-field is trapped at the origin at temperature T ≳ √ λfa/y for yfa < T and T ≳ λ 1/4 fa/g for yfa > T .For couplings of order unity, fa < T < Tmax ≃ 6.57×10 15 GeV is the maximum reheating temperature bounded by the inflationary scale and assuming instantaneous reheating.Nonetheless, if λ is small (corresponding to a small radial-mode mass), the bound can be weakened.
As the Universe cools, the axion mass develops due to non-perturbative effects (like strong confinement in the case of the QCD axion).The second term in V (Φ) breaks explicitly the U (1) discretely, leading to sheet-liked defects or domain walls, attached to the cosmic strings.The domain wall is characterized by its surface tension σ ≃ 8m a f 2 a /N 2 DW [42].The axion field starts to feel the presence of the walls when 3H ≃ m a .The domainwall network can be stable or unstable depending on the number of domain walls attached to a string.The value of N DW is very UV-model-dependent.It can be linked to the discrete symmetry Z NDW [86][87][88] that remains after the confinement of the gauge group that breaks the global U (1) symmetry explicitly and generates the axion mass.This occurs at the scale Λ ≃ √ m a F a , where F a = f a /N DW , that is when the domain walls are generated, attaching to the existing cosmic strings.
For N DW > 1, the string-wall system is stable and long-lived.Its decay may be induced by V bias , the biased term [89][90][91], which could be of QCD origin [30,40].This decay is desirable to prevent DW from dominating the energy density of the Universe at late times.V bias is therefore an additional free parameter beyond m a and f a that enters the GW prediction in the case where N DW > 1.

A. Case i) NDW = 1
If only one domain wall is attached to a string, i.e., N DW = 1, the string-wall system quickly annihilates due to DW tension when3 H(T dec ) ≃ m a [41].The cosmic string SGWB features an IR cut-off corresponding to the The symbols correspond to the benchmark points in the axion parameter space in Fig. 3 (with T⋆ = {128 MeV, 10 2 GeV} for { , ♣}).The bestfitted spectra to the PTA data are the blue-⋆ curve for global strings, corresponding to {fa, ma} ≃ {9.9 × 10 15 GeV, 4.8 × 10 −15 eV} which is excluded by the axion dark matter abundance [see Eq. ( 14)], and the green-curve for domain walls (with maF 2 a = 2.6 × 10 15 GeV 3 ).The power-law integrated sensitivity curves of GW experiments [14,[92][93][94][95][96][97][98][99][100][101][102][103][104][105] are taken from [12,106].For fixed {ma, fa} values, the peak of the DW-GW spectrum moves along the dashed line as T⋆ varies; see Eq. ( 12). ( The cut-off position -frequency and amplitude -can be estimated with Eqs. ( 2)-( 5).At f GW < f CS GW (T dec ), the spectrum scales as Ω GW ∝ f 3 GW due to causality.Note that for m a ≪ 10 −16 eV, the cut-off sits below nHz frequencies, and within the PTA window, we recover the same GW spectrum as the one in the limit m a → 0. Our analysis applies the numerical templates of the globalstring SGWB -covering the ranges of f a and T dec priors.We calculated these templates numerically by solving the string-network evolution via the velocity-dependent onescale model [68,[108][109][110][111], shutting off the loop production after H = m a , and calculating the SGWB following Ref.[12].
Attached to a string, N DW walls balance among themselves and prevent the system from collapsing at H ≃ m a [41,112].The domain-wall network later evolves to the scaling regime where there is a constant number of DW f  per comoving volume V ≃ H −3 .The energy density of DW is ρ DW ≃ σH −2 /V ≃ σH and it acts as a longlasting source of SGWB [89,[113][114][115][116][117][118]; cf.[119] for a compact review.The network red-shifts slower than the Standard Model (SM) radiation energy density and could dominate the Universe.The biased term V bias -describing the potential difference between two consecutive vacua -explicitly breaks the U (1) symmetry and induces the pressure on one side of the wall [8,89].Once this pressure overcomes the tension of the wall4 , the string-wall system collapses at temperature, The fraction of energy density in DW is maximized at this time and reads, The energy density emitted in GW is [42] where we fix ϵ ≃ 0.7 from numerical simulations [117].It reaches its maximum at T ⋆ .The spectrum exhibits the broken-power law shape and reads, where the normalized spectral shape is, The f 3 GW -IR slope is dictated by causality, the UV slope f β GW is model-dependent, and the width of the peak is δ.The peak frequency corresponds to the DW size, i.e., the horizon size f ⋆ GW ∼ H ⋆ [117].Its value today reads, 10.75 g * s (T ⋆ ) From Eqs. ( 7), (9), and (11), each value of m a f 2 a corresponds to a degenerate peak position of the GW spectrum, which are shown as the dashed line in Fig. 1.
The DW can decay into axions, which either behave as dark radiation or decay into SM particles.When DW decay into dark radiation, the ∆N eff puts a bound α ⋆ ≲ 0.06 [48], i.e., the peak of GW spectrum has h 2 Ω GW ≲ 10 −9 (which cannot fit the whole 14 bins of NG15 data).As α ⋆ controls the amplitude of the GW spectrum (9), we consider a larger range of α ⋆ , up to α ⋆ = 1 when the energy density of DW starts to dominate the Universe.To get around the ∆N eff bound, we will therefore consider the case where the axions produced by domain walls eventually decay into SM particles.
In this paper, we confront the most recent PTA data for both cases: i) N DW = 1 where the SGWB in the PTA range dominantly comes from the cosmic strings, and ii) N DW > 1 where the SGWB in the PTA range comes from the domain walls.These two cases correspond to axions of two utterly different mass ranges.For case i), the cosmic strings live long; that is, m a is small.Instead, the case ii) corresponds to the large m a region.We compare the GW spectra in Fig. 1 for different benchmark points, corresponding to locations in the {m a , F a } plane are shown in Fig. 3.

III. SEARCHING AND CONSTRAINING SGWB WITH PTA
This work analyzes the recent NG15 data set [125] covering a period of observation T obs = 16.03 years [1].From the pulsar timing residuals, the posterior probability distributions of the global-string and domain-wall model parameters are derived.We consider 14 frequency bins of NG15 data, where the first and last bins are at 1/T obs ≃ 1.98 nHz and 14/T obs ≃ 27.7 nHz, respectively.The analysis is done by using ENTERPRISE [126,127] via the handy wrapper PTArcade [46,47].The priors for the model parameters are summarized in Tab.I in Appendix A. We refer readers to Ref. [5] for a short review of Bayesian analysis.
This work considers the SGWB in the two scenarios discussed above together with the astrophysical background.Fig. 2-middle and -right show the 68%-CL (or 1σ) and 95%-CL (or 2σ) in dark and light blue regions, respectively.We obtain the best-fit values f a = 9.87 MeV for domain walls.The global-string and domainwall SGWB are preferred over the SMBHB signal implemented by PTArcade, as suggested by their Bayes Factors (BF) larger than unity (BF CS = 26.0,BF DW = 44.7)when compared to the SMBHB interpretation; cf.Eq. ( 9) of [5].We show the best-fitted spectra for these two new-physics cases in Fig. 2-left.Translating into axion parameters via Eq.( 4) and ( 7), the best fits correspond to {f a , m a } = {9.87× 10 15 GeV, 4.78 × 10 −15 eV} for global strings (excluded by the axion overabundance) and m a F 2 a = 2.6×10 15 GeV 3 for domain walls.For completeness, we show the case without the SMBHB contribution in App. C. Because the two new-physics cases explain the data well by themselves, we see that the 1σ and 2σ regions of Fig. 2 match those without the SMBHB in Fig. 5.The values of the best fits, given in App.C, only change slightly.
Although the two scenarios could explain the signal, this work aims to set bounds on the model parameter space associated with a too strong SGWB in conflict with the NG15 data.Following [5], we identify excluded regions of the new-physics parameter spaces using the posterior-probability ratio (or K-ratio).Specifically, the excluded gray regions in Fig. 2-middle and -right correspond to the areas of parameter spaces where the Kratio between the combined new-physics+SMBHB and the SMBHB-only models drops below 0.15 , according to Jeffrey's scale [128], due to a too-strong SGWB from the new-physics model.We emphasize that the values of the BFs strongly depend on the modeling of the SMBHB signal as it is the ratio of evidence of the considered model and the SMBHB template.However, the constrained regions depend only slightly on it [5].
We emphasize that the constraints on the axion parameter space presented in this paper are not the same as the regions of best-fit obtained in the literature using the previous dataset, e.g., [48,50].For fitting the PTA data, a particular part of the GW spectrum is preferred; thus, the best-fits region is allowed within a tight parameter space (the blue blobs in Fig. 2).On the other hand, the constraint can be drawn from any part of the spectrum if the GW signal becomes too large and disfavored by the data.So, the constraint can be extended over a vast parameter space (the grey regions in Fig. 2).Now we discuss, in turn, the NG15 constraints -on global strings (N DW = 1) and domain walls (N DW > 1) -and translate them into the constraints in the axion parameter space.

A. Result for NDW = 1, implications for light axions
We fit the PTA data with the global-string SGWB, varying {f a , T dec }.The 2D posterior result is shown in Fig. 2, and the dark-blue region is where the cosmicstring SGWB dominates and fits the data to the significance of 1σ with the best fit {f a , m a } ≃ {9.9 × 10 15 GeV, 4.8 × 10 −15 eV}, shown as the benchmark case ⋆ in Figs. 1 and 3. Note that this benchmark point is excluded by the axion overabundance constraint [see Eq. ( 14)].A too-large global-string SGWB is constrained by PTA in the grey region of Fig. 2-middle.For small f a , the GW from cosmic strings cannot fit the data as its amplitude becomes too small.
As T dec ≪ 0.1 MeV (m a ≪ 10 −17 eV), the cut-off (5) associated with T dec moves below the PTA window (f GW (T dec ) < nHz).The constraint in this case, Fig. 2middle, reads f a < 2.8 × 10 15 GeV (m a -independent), which is stronger than the LIGO bound6 (f a ≲ 8 × 10 16 GeV).For completeness, we also analyzed the case of stable global strings (i.e., m a → 0) in App.D, and we obtained a similar bound.For T dec ≫ 0.1 MeV (m a ≫

Z a
FCC-ee FIG. 3. PTA limits (in green) on postinflationary axions, compared to existing experimental constraints as compiled from AxionLimits [120] and to theoretical bounds: dark radiation overabundance ∆N eff bound (13) as dashed horizontal line and ALP overabundance (14) in the shaded grey region.Fa = fa/NDW.The orange dotted lines in the ma ≳ 1 GeV region are the projections of future collider experiments, LHC (h → Za) and FCC (e + e − → ha), obtained from [121,122] with the maximally allowed ALP-SM coupling.The red region denoted ma > Fa is where the axion effective field theory is not valid.The comparison with experimental bounds uses g θγγ = 1.02αEM/(2πFa) ≈ 2.23 × 10 −3 /Fa for the relation between the photon coupling and Fa, as motivated by KSVZ models [123,124].The recent PTA data [1] excludes the green small-ma region due to cosmic-string SGWB (NDW = 1).It also potentially excludes the high-ma region due to domain-wall SGWB for NDW > 1, depending on the value of T⋆.The other green band at large ma is the region that PTA can constrain if T⋆ varies in the range MeV < T⋆ < 302 MeV, as illustrated in Fig. 4. The two benchmark points {⋆, ♠} correspond to cosmic-string SGWB, and the two black benchmark lines { , ♣} correspond to the domain-wall SGWB, whose spectra are shown in Fig. 1.The green dot-dashed line is explained in App.B.     4. The PTA-DW constraint (in green) changes with T⋆.For fixed T⋆ and ma, the constrained range of Fa in green is derived from the α⋆ constrained region of Fig. 2-right, using Eq.(7).The yellow region corresponds to α⋆ > 1, which corresponds to the DW domination and can change the GW prediction; we do not extend the constraint into this region.For T⋆ ≳ 302 MeV (cf.Fig. 2-right), NG15 data constrains α⋆ > 1; that is, the green band overlays part of the yellow region.The blue region is where the axions -produced from DW annihilations -dominate the Universe before they decay prior to BBN.In this case, the theoretical prediction for the GW spectrum also has to be re-evaluated.
10 −17 eV), the cut-off sits at a frequency higher than the PTA window, and the SGWB signal is dominated by the IR tail signal, which scales as Ω GW ∝ f 3 GW .From Eqs. (1) and (2), we obtain the asymptotic behavior of a ), up to the log correction in Eq. ( 2), toward large f a limit.We show this bound (green-region) in the usual axion parameter space in the bottom-left corner of Fig. 3.The NG15 constraint on f a values for N DW = 1 corresponds to f a > H inf /(2π).Therefore, it does not apply to cosmic strings linked to quantum fluctuations during inflation.Note that Eqs. ( 1) and (2) assume a standard cosmological history, i.e., a transition between the radiation era and the matter era occurring at T eq ∼ 1 eV.In the region of parameter space where the axion abundance from the string network exceeds the dark matter abundance [see Eq. ( 14)], the matter era starts earlier, and the cosmological evolution is not viable.The non-standard cosmological history will modify the PTA data (e.g., the calibration of pulsar timing data and the dispersion measure) and also the SMBHB modeling [129].Ignoring its impact on PTA data, we can still estimate how the axion overabundance affects our constraint, just from the dilution effect on the GW spectrum [12,21,22]; see Eq. (B3) in App.B. In Fig. 3, the dot-dashed green line shows the modified PTA constraint due to the diluted GW spectrum from the axion overabundance; see App.B for the estimate of the scaling.
∆N eff & dark matter constraints.-Althoughthe PTA constraint excludes a large region of the axion parameter space, there exist other theoretical bounds.Axionic strings are known to emit axion particles dominantly [31].Depending on its mass, the axion can contribute to either dark radiation or cold dark matter.Axions that are relativistic at the time of Big Bang Nucleosynthesis (BBN) are subject to the dark radiation bound expressed as a bound on the number of extra neutrino species, ∆N eff < 0.46 [130].There are uncertainties in deriving this bound linked to the log-correction to the number of strings in the global-string network evolution [35][36][37]84].In this paper, we quote two bounds: the one relying on the semi-analytic calculation [22] by Chang and Cui (CC), and the lattice result [23] by Gorghetto, Hardy, and Nicolaescu (GHN): where we implicitly assume λ ∼ 1 for the GHN bound and H BBN ≃ 4.4×10 −25 GeV is the Hubble parameter at BBN scale (T BBN ≃ MeV).Since ALPs have a small mass at late times, they behave as cold dark matter.Subject to the uncertainty in simulations [23,37,81], the abundance Ω a of axion dark matter from strings predicted by GHN sets a constraint on the axion, typically ξ * ≈ 25 and x 0,a ≈ 10 [23].Note that the collapse of the string-wall system7 at H ≃ m a produces an axion abundance of the same order as the one from strings [36], therefore an O(1) correction is expected in Ω a in Eq. 14.We show both dark radiation and dark matter bounds in Fig. 3.We see that the PTA constraint becomes competitive with the equivocal ∆N eff bound for m a ≲ 10 (−22,−23) eV.
Effects of non-standard cosmology.-Sofar, the standard ΛCDM cosmology [130] has been assumed.On the other hand, alternative expansion histories to the usually assumed radiation era are not unlikely above the BBN scale, such as a period of matter domination or kination resulting in a strongly different spectrum of GW for cosmic strings [12,21,22,76,132].Nonetheless, the non-standard cosmology modifies the cosmic-string GW spectrum in the high-frequency direction.From Eq. ( 1), the non-standard era must end below the MeV scale to substantially change the SGWB in the PTA window.We have checked the effects of matter and kination eras with PTArcade and found that such SGWB distortion cannot improve the global string interpretation of PTA data.Besides, we expect only a negligible effect on the PTA bound obtained in this work.
QCD axion.-FromFig. 3, the PTA data can exclude some parts of the QCD axion (red line).However, this region of parameter space is already excluded due to the overabundance of axion dark matter or due to ∆N eff bounds.To relax these bounds, one can invoke a scenario where cosmic strings decay during a matter-domination era (or any era with the equation-of-state smaller than that of radiation), which efficiently dilutes these relics but still allows for a GW signal in the PTA frequency range [22,24,25].Interestingly, such matter-domination era at early times can imprint a specific signature in the SGWB from global strings, which can be observed in future-planned GW experiments at frequencies above nHz frequencies [12,21,22,133].

B. Result for NDW > 1, implications for heavy axions
We fit the DW SGWB, varying {T ⋆ , α ⋆ , β, δ}, to the PTA data.Because the posteriors of β and δ are unconstrained, we show only the 2D posterior of {T ⋆ , α ⋆ } in Fig. 2-right.The DW SGWB can fit the PTA data in the dark-blue region to 1σ.The best fit value of {T ⋆ , α ⋆ } is translated via Eq.( 7) into m a F 2 a ≃ 2.6 × 10 15 GeV and corresponds to the benchmark spectrum and line in Figs. 1 and 3, respectively.However, for large enough α ⋆ , DW generates a GW signal well stronger than the PTA signal, leading to a constraint in the gray region in Fig. 2-right.The constraint is the strongest α ⋆ ≲ 0.02 at T ⋆ ≃ 13.8 MeV when the peak of the SGWB is centered in the PTA window; see also Eq. ( 11).For T ⋆ > 13.8 MeV (< 13.8 MeV), the GW spectrum has its IR (UV) tail in the PTA range; thus, the constraint on α ⋆ becomes weaker.
For heavy axions with Z NDW -symmetry whose mass depends on the explicit-symmetry-breaking scale Λ ≃ √ m a F a where F a = f a /N DW , the PTA constraint in Fig. 2-right is translated via Eq.( 7) into a bound on {F a , m a } with the degeneracy among them.For a fixed T ⋆ , we obtain the excluded region on the axion parameter space, i.e., the green region of Fig. 4. Very large f a corresponds to α ⋆ > 1; the DW-domination era occurs before it decays and should affect the GW prediction.We do not extend our PTA bound in the DW domination regime, shown in the yellow of Fig. 4. In fact, Eq. ( 9) assumes a radiation-dominated Universe.Constraining the DW-domination region requires computing the evolution of the DW network and its SGWB in a Universe with a modified equation of state.We leave this nontrivial task for future investigation; see also [134].To be conservative, we mark this region unconstrained for now, although we expect some constraints will prevail there.
Because the PTA constraint on α ⋆ is not linear in T ⋆ , the width of the PTA band is maximized only for T ⋆ ≃ 13.8 MeV where the bound on α ⋆ is the strongest.In Fig. 3, we also show the ability to constrain axion parameter space with the PTA-DW signal.We obtain the constraint by summing the excluded regions for the range MeV ≤ T ⋆ ≲ 302 MeV, where T ⋆ ≃ 302 MeV is where the constraint has α ⋆ ≥ 1 in Fig. 2-right.The upper limit of the green region (large-m a ) of Fig. 3 is set by the constraint at T ⋆ = MeV: α ⋆ ≳ 0.2; see Fig. 2right.Using Eq. ( 7), this upper bound is defined as m a F 2 a ≳ 2 × 10 11 GeV 3 .Some regions above and within the green band (smaller m a F 2 a ) will be probed by future particle physics experiments [121,122,[135][136][137].
Other than the PTA bound, the {T ⋆ , α ⋆ } parameter space is subject to theoretical constraints related to the DW decay and its by-products.In this work, we consider that the heavy axion produced from the DW decay subsequently decays into SM particles, e.g., photons via L ⊃ − g 4 F F θ with the decay rate Γ θγ = m 3 a g 2 /(64π) [138].Using this to F a = 1.92αEM /(2πg θγ ), the decay is efficient when Γ θγ > H(T ), which is equivalent to, The bound T BBN < T θγ is similar to the BBN bound from [139] in Figs. 3 and 4.Moreover, the heavy axion which behaves nonrelativistically might decay after it dominates the Universe if T ⋆ > T dom > T θγ where the temperature T dom corresponds to the heavy-axion domination, i.e., We mark this region in the blue region of Fig. 4. For the sum of PTA constraints varying T ⋆ in Fig. 3, we omit showing the color of the axion matter-domination (MD) region, which cuts the PTA region from the low-m a region 8 .This heavy axion induces a matter-domination era that would change the GW prediction, e.g., the causality tail of the spectrum gets distorted [140][141][142].Although this spectral distortion would change the fitting of the data, it would affect the constraint derived here minimally for two reasons.First, the blue region in Fig. 4, leading to the axion-MD, is smaller than the constrained region.Second, within this region, we find T dom /T θγ ≲ 10 which leads to Ω GW ∝ f GW for frequencies in the range [10 −2/3 , 1] × f DW p , using Eq.(4.5) of [140] where f DW p is the peak frequency (11).Other effects.-Thefriction from axionic DW interactions with particles of the thermal plasma could change the network's dynamics [143] and potentially the SGWB spectrum.Another effect that could change the bounds is the potential collapse of DW into primordial black holes [44,45,[144][145][146][147][148].Nonetheless, since the prediction is based on the spherical collapse, we would need a largescale numerical simulation of DW to check whether the PBH formation can be realized.Lastly, further QCD effects can impact the DW decays relevant for PTA [134,149,150].

IV. CONCLUSION
We analyzed the consequences of the 15-year NANOGrav data on the parameter space of postinflationary axions.The bounds in Fig. 3 come in two distinct regimes: the low and large axion mass ranges, which are respectively associated with signals from axionic global strings (N DW = 1) and domain walls (N DW > 1).In the low-axion-mass region, the constraint on f a is strongest for m a ≪ 10 −17 eV, and reads f a < 2.8 × 10 15 GeV.It is competitive with the ∆N eff bound.At high masses, 0.1 GeV ≲ m a ≲ 10 3 TeV, a substantial region, corresponding to m a (f a /N DW ) 2 ≳ 2 × 10 11 GeV 3 , can be excluded for DW decaying in the T * ∝ √ V bias ∼ 1 − 300 MeV range.
This study motivates the investigation of the SGWB in the regime of DW domination, as this knowledge could lead to substantial new constraints at large m a and f a values.Once the network of DW dominates the Universe, the scaling regime might be lost.DW would instead enter the stretching regime [151] where the energy density scales as ρ ∝ a −1 , the equation of state of −2/3 leading to the accelerated cosmic expansion could be in tension with several cosmological observations [134,152].Moreover, a period of early DW domination together with the axion matter domination can also affect the SGWB spectra from DW and cosmic strings [12, 140-142, 153, 154].
To conclude, GW is a promising tool to probe axion physics.PTA measurements have opened the possibility of observing the Universe at the MeV scale, enabling us to constrain several classes of axion models.By combining NG15 with other data sets from EPTA, InPTA, PPTA, and CPTA collaborations, the constraints on axions can become more stringent, similar to what has been shown for other cosmological sources [155,156].Other planned GW observatories will permit the search for different parts of the predicted SGWB from axion physics and probe the axion parameter spaces uncharted by the PTA.Moreover, the synergy of GW experiments over a wide frequency range will allow us to distinguish the axion-GW signals from other SGWB from astrophysical and cosmological sources [157].
where we used a −3 ∝ g * s (T )T 3 .The domination before the radiation-matter equality T ′ dom > T eq ≃ 0.75 eV leads to dark matter overabundance.This bound is similar to Eq. ( 14) and the gray region denoted "DM strings" in Fig. 3.A universe with T ′ dom > T eq cannot resemble the standard ΛCDM model.Below, we compute the modified GW spectrum from global strings due to the earlier matter era, although we do not use it for analyzing the PTA data which relies on the standard cosmology assumption, for e.g. the calibration of pulsar timing data and the dispersion measure [129].
We set today's time as when the photon temperature matches the CMB observation.The GW signal emitted with frequency f emit GW at temperature T (a emit ) has the frequency today [f emit GW (a emit /a 0 )] which is the same for ΛCDM and non-ΛCDM cases, i.e., (a emit /a 0 ) = (a emit /a 0 ) ΛCDM .On the other hand, the emitted GW energy density [Ω GW ∼ (ρ emit GW /ρ tot,0 )(a emit /a 0 ) 4 ] gets diluted as ρ tot,0 > ρ ΛCDM tot,0 [12].We define the dilution factor Υ as

94 <
l a t e x i t s h a 1 _ b a s e 6 4 = " q n FIG. 2. left:The SGWB spectra from global strings and domain walls + SMBHBs, providing the best-fits to the PTA data and corresponding to {fa, ma} ≃ {9.9 × 10 15 GeV, 4.8 × 10 −15 eV} for global strings and maF 2 a = 2.6 × 10 15 GeV 3 for domain walls (in violins, taken from[5]).middle and right: 1σ (dark blue) and 2σ (light blue) regions of the likelihood of the globalstring/domain-wall parameters, assuming the template of global-string/domain-wall + SMBHB backgrounds.The gray region is excluded due to too strong GW signals from global strings/domain walls that conflict with PTA data.The region above the black dashed line in the middle panel (including the best fit) conflicts with the dark matter abundance [see Eq. (14)].

6
GeV f a /N DW .