Challenges in Interpreting the NANOGrav 15-Year Data Set as Early Universe Gravitational Waves Produced by ALP Induced Instability

In this paper, we study a possible early universe source for the recent observation of a stochastic gravitational wave background at the NANOGrav pulsar timing array. The source is a tachyonic instability in a dark gauge field induced by an axion-like particle (ALP), a known source for gravitational waves. We find that relative to the previous analysis with the NANOGrav 12.5-year data set, the current 15-year data set favors parameter space with a relatively larger axion mass and decay constant. This favored parameter space is heavily constrained by $\Delta N_{\rm eff}$ and overproduction of ALP dark matter. While there are potential mechanisms for avoiding the second problem, evading the $\Delta N_{\rm eff}$ constraint remains highly challenging. In particular, we find that the gravitational wave magnitude is significantly suppressed with respect to the gauge boson dark radiation, which implies that successfully explaining the NANOGrav observation requires a large additional dark radiation, violating the cosmological constraints. Satisfying the $\Delta N_{\rm eff}$ constraint will limit the potential contribution from this mechanism to the observed signal to at most a percent level.

Axions are pseudo-Goldstone bosons first proposed as the solution to the strong CP problem [58,59], and their quanta were soon realized to also serve as viable dark matter (DM) candidates [60][61][62][63].ALPs, as suggested by their name, inherit and generalize the cosmological phenomenology of axions, while do not necessarily have a connection to the strong CP due to its shift symmetry-breaking mass.From now on, we will use axions and ALPs interchangeably, with no connection to the strong CP problem.ALPs can be responsible for the cosmic inflation and reheating [64][65][66], address the hierarchy problem [67][68][69] and also generically arise in various string theory models [70,71].Specifically, when coupled to massless or very light gauge modes like the dark photons, the ALP field's rolling can exponentially amplify a specific helicity of these gauge modes, a process known as "tachyonic instability".This can serve as an efficient mechanism for particle production in the universe.This mechanism has been studied in the context of reheating [66,72], production of dark photon DM [73][74][75] and depletion of axion DM to avoid over-closure [76].Recent studies show that the drastic enhancement of the vector fields through this mechanism can conveniently serve as a source of GWs that may be visible on interferometers or pulsar timing arrays (PTAs) [77,78], and generate CMB B-modes if the ALP is light enough [79].
In this work, we analyze the interpretation of the recent GW observation by the NANOGrav collaboration in their 15-year data set [1] (hereafter NG15) as primordial GWs induced by tachyonic instability at around MeV temperatures 1 .This interpretation was considered by Ref. [83] for the earlier NANOGrav 12.5-year data set (NG12) [84] and the IPTA's second data release [85].The authors find that the combined data favors a region that is in mild tension with ∆N eff constraints, i.e. the constraints on the effective light degrees of freedom at the big-bang nucleosynthesis (BBN) [86] and cosmic microwave background (CMB) [87].In this paper, we will show that with the updated NG15 data, this interpretation is in severe tension with ∆N eff constraints, which is virtually unavoidable as we explain below.This is explained by the NG15 data preferring a higher peak magnitude, especially with the improvement of the data at higher frequency bins.
In Ref. [83], the analysis is carried out by using a formula fitted to their lattice result (see also Ref. [88]).In this analysis, we use both the linear calculation and the fitted lattice spectrum, which differ significantly in the predicted frequency of the gravitational wave, but not in its magnitude.Therefore, while the parameter regions favored by the NG15 data in our two approaches are different, the expected ∆N eff values for both parameter regions are very similar.The exclusion we derive here is therefore robust.
In the following, we explain the two major constraints on this scenario.These are the overclosure constraints from the massive ALP field, and the ∆N eff constraints from the massless gauge modes.These two constraints arise from the same origin: to produce the observed GW signal, the initial energy density of the ALP field must be large enough.While a substantial portion of this energy density is eventually converted into dark gauge field quanta, the remaining energy density of the ALP field undergoes matter-like redshift.However, across the entire parameter space, this energy density significantly exceeds the observed DM relic abundance, resulting in an over-closure constraint.It is conceivable that this constraint could be circumvented by transferring the matter-like energy density into radiation or kinetic energy, which undergoes faster redshift, through additional model building.The ∆N eff constraint comes from the energy budget in the dark gauge field quanta that behave as dark radiation.As we show in this work, the energy density of gauge quanta required to produce the GW signal far exceeds the constraint on ∆N eff .Since both the GW and the dark radiation densities redshift in the same way, it is much harder to avoid this constraint.Consequently, the aforementioned explanation of the NG15 signal is disfavored, and potential GW signals from the ALP-induced tachyonic instability are expected to be a few orders of magnitude below the observed signal.
Despite no significant evidence being found by the analysis of the NANOGrav Collaboration at the current experimental sensitivity [89], an anisotropic GW background remains a very interesting possibility and could naturally arise from, e.g. the fluctuations of the GW sources like supermassive black hole mergers [90].Here we ignore this possibility and base our analysis on the assumption of an isotropic GW signal.
The organization of the paper is as follows: In Sec. 2 we first briefly review the ALP model, the tachyonic production of dark radiation from the dynamic of ALPs, and how GWs are sourced in this process.Then in Sec. 3 we describe the setup of our numeric calculation and show the model parameter space where the produced GW reasonably explains the NG15 results.In this context, we discuss the severity of the constraints from the cosmic dark matter relic abundance and ∆N eff .After this, we discuss some future perspectives on this model and conclude.
2 The Instability-induced Gravitational Waves

Tachyonic production of dark photons
Consider the Lagrangian of an axion coupling to the gauge boson of a U (1) gauge symmetry where ϕ is the axion field, f is the axion constant, X µν = ∂ µ X ν − ∂ ν X µ is the exterior derivative of the massless dark photon field X µ and α is the axion-dark photon coupling.The axion potential is chosen to have the generic quadratic form of V (ϕ) = 1 2 m 2 ϕ 2 where m is the axion mass, which can be seen as the leading order expansion of the axion-like potential Λ 4 cos(ϕ/f ).At the leading order, we ignore the spatial inhomogeneity of the axion field so that the axion field ϕ does not carry a momentum dependence.We discuss the full solution including the ALP inhomogeneities in Sec.2.3.The equation of motion of ϕ can then be written as where a is the scale factor of the FRW metric, and H is the Hubble parameter.The prime ′ indicates a derivative with respect to the conformal time τ .The E and B on the right-hand side are the corresponding electric and magnetic fields of the dark photons.As in the case of free axions, ϕ receives a Hubble damping when the Hubble parameter is large, and the field only starts rolling when H ∼ m.As a convention, we define the oscillation time τ osc to be the conformal time where H(τ osc ) = m, and a osc to be the scale factor at that time.At the beginning of the rolling, the E • B term is negligible, and the axion dynamics is the same as that of free axions.After the initiation of tachyonic production of the dark photon, as explained in detail below, the E • B term will introduce increased friction to the axion motion, leading to energy transfer from axions to dark photons.The dark photon field is quantized as where we take the Coulomb gauge and write Dk ≡ d 3 k/(2π) 3 .The polarization vector ϵ ± (k) [72].Due to the spatial isotropy, the mode functions v ± should depend only on the magnitude of the wave number k = |k|, and they have the equation of motion with the dispersion relation ω 2 ± (k, τ ) = k 2 ∓ kαϕ ′ /f .In Eq. ( 4), it is evident that when the axion is rolling, one helicity of the gauge modes with momentum 0 < k < α|ϕ ′ |/f will become tachyonic, i.e. ω 2 < 0, resulting in an exponential amplification.The amplification of the v +(−) modes occurs when ϕ ′ is positive (negative), and the enhancements on the two helicities alternate as the axion field rolls around the minima of V (ϕ).In this context, it is reasonable to expect that the helicity experiencing initial enhancement during the ALP field's rolling motion would be more amplified compared to the other, as it spends a longer duration within the tachyonic band.
However, the tachyonic enhancement of the dark photons is not perpetual.As the axion field undergoes redshift and spends less time rolling on either side of the potential, the available time for the gauge modes to amplify decreases, even if their momenta fall within the tachyonic band.The growth time scale of a gauge mode in the tachyonic band is about |ω ± | −1 , while the axion oscillation time scale is about 1/(am) as estimated from free rolling.For a gauge mode to be enhanced we should require the former time scale to be shorter, such that there will be enough growth time for the gauge mode.This suggests that the tachyonic production process should cease at the time τ s where a(τ s )/a osc ∼ (α/2) 2/3 [77].
The left panel of Fig. 1 illustrates the evolution of this coupled ALP-dark photon system, where the energy density evolution of the ALPs (ρ ϕ ), dark photons (ρ X ) and GWs (ρ GW ) are shown.The energy densities are normalized against the cosmic energy density of the standard ΛCDM cosmology.The ALP field starts rolling at a osc and behaves as a freely rolling field before efficient radiation production occurs.It takes a few oscillations for the axion field to release a significant amount of its energy into GW and dark photon radiation (DR), depending on the size of α.After the tachyonic production ceases, the remaining ALPs and dark photons start to redshift like ordinary matter and radiation.When solving the equations of motion for the ALP and dark photons, we adopt the assumption of a standard ΛCDM evolution for the Hubble parameter H and scale factor a. For the majority of the parameter space favored by the NG15 GeV, α = 50.Left: The evolution of the ALP (blue), dark photon (orange), and GW (red) energy density, normalized by the total energy density of the standard ΛCDM cosmology.The dashed green line is proportional to a, showing that ALPs now evolve like matter.Soon after the ALP starts rolling, at a/a osc ≈ 20 the tachyonic production of gauge modes is maximized and the ALP quickly dumps its energy into dark radiation.After the tachyonic production, the ALP and dark photon continue to behave like ordinary matter and radiation.On the other hand, this set of parameters has the ALP relic abundance overclose the universe and has too many radiation degrees of freedom at the BBN and CMB.Right: The GW generated from the tachyonic instability (red curve), compared with the first 14 frequency bins of NG15's free spectrum density (gray shades), taken from Ref. [1].
measured spectrum, this assumption holds reasonably when the tachyonic production is active.When f and m are large and α is small, on the other hand, there may be a mild violation where ρ ϕ /ρ ΛCDM ≳ O(1) during or at the end of the production.While note that this is also the parameter space suffering more strongly from the over-closure issue (see Fig. 4 in later texts), we do not expect our final conclusion to change significantly.

The instability-induced GW spectra
As the gauge modes undergo tachyonic amplification, they begin to act as sources for the SGWB.The generated metric tensor perturbation is given by where Π ij is the transverse traceless energy-momentum tensor of the gauge modes, and G(k, τ, τ ′ ) is the Green's function of the tensor perturbation's equation of motion.As our signal is generated during radiation domination, we explicitly have G(k, τ, τ ′ ) = sin(k(τ − τ ′ ))/k.The GW spectra are then calculated as where ⟩ is the unequal time correlator of the gauge modes' energy-momentum tensor.It can be shown that [77 where is the projector function, and v λ are the quantized mode function as defined in (3).We are thus well-equipped to calculate the GW spectra.
The observed signals on PTAs correspond to the SGWB at the present time.Therefore, in Eq. ( 6), we should in principle integrate until the current age of the Universe, τ = τ 0 .On the other hand, as discussed in the previous subsection, the tachyonic production of the gauge modes ceases at around τ ∼ τ s , and the active sourcing of GW stops thereafter.Hence, when evaluating Eq. ( 6), we only need to integrate up to a suitable point τ max after τ s .In other words, we approximate Π , where θ(z) represents the step function.After τ = τ max the GW simply behaves like normal relativistic degrees of freedom and redshift as a −4 , i.e., ρ GW (k, τ ) = (a(τ max )/a(τ )) 4 ρ GW (k, τ max ), as indicated by the pre-factor in Eq. ( 6).The GW frequency f GW redshifts as a −1 and is related to the co-moving wave number k as k = 2πaf GW .
In the right panel of Fig. 1 we present the generated GW, and compare it with the first 14 frequency bins in the NG15 common-spectrum spatially-uncorrelated red noise-free spectrum, following the analysis in Ref. [1].If we only consider the GW spectrum and ignore other constraints, then the instability-induced GW is well capable of explaining the NG15 result.The mass m and the coupling parameter α determine the peak momentum mode of the tachyonic band which translates to the peak frequency of the GW power spectrum.Ref. [78] suggests the peak frequency to be given by Using the above formula and the peak of the observed NG15 spectrum one can infer the relevant axion mass for instability-induced GWs to be a SGWB explanation.For 40 < α < 100 of our interest, we expect the ALP mass to lie approximately around m ≈ [1.5 − 4.0] × 10 −12 eV.The peak amplitude of the GW spectrum on the other hand depends on the decay constant f but not on the mass.At the onset of ALP rolling the total energy in the sector is given by ρ ϕ = m2 f 2 /2, and the cosmic energy density is ρ tot = 3H 2 M 2 pl ∝ m 2 as H ∼ m, and hence ρ ϕ /ρ tot ∝ f 2 is independent of m.During tachyonic production, the ALP energy is dumped into dark photons, and thus we may reasonably expect ρ X /ρ tot ∝ f 2 .This implies that the peak amplitude of Ω GW , which is the two-point correlator of the tensor metric perturbation, scales as f 4 and is mostly independent of m [77,79], as it is sourced by ρ 2 X (as shown in Eq. ( 6)).This gives a nice complementarity between these two parameters, i.e., f mainly probes the amplitude, and m sets the peak frequency.

Going beyond linear calculation
The axion equation of motion ( 2) is at the linear order, which assumes the axion field to be homogeneous, and thus the evolution only depends on time.However, the tachyonic production of the dark photon is momentum-dependent and thus inhomogeneous.The back reaction from the exponentially enhanced gauge modes could then generate inhomogeneity in the axion field through the ϕX X term and make ϕ momentum dependent.The inclusion of the ALP momentum dependence significantly complicates the equations and is therefore usually solved on the lattice [88], while as a result, the GW spectrum sourced during the process may differ compared with the linear order results.
From the particle point of view, the effect of the back reaction is to introduce backscattering (two dark photons into one ALP) and multiple scattering into the linear system, which shuffles energy between different modes.The spectra of ALPs and dark photons are thus broadened.In terms of the generated GW spectrum, these effects will lead to higher peak frequency and a more "relaxed" spectral shape: having a smaller power-law index at lower frequencies and dropping off slower at higher frequencies.The amplitude of the GW spectrum, on the other hand, receives no drastic changes.An empirical GW profile can be obtained from fitting to the lattice calculation [83], which reads A comparison between this empirical lattice fit GW spectrum (denoted as "ELF" hereafter) and our linear order numeric spectrum is presented in Fig. 2. It is important to note that this empirical expression is fitted up to 60 < α < 100. 2 Also, for our linear order calculation, we find that the GW signal strength drops for α < 50 due to a weakened energy transfer from the GeV, α = 50, the same as those in Fig. 1.The linear calculation gives an integrated Ω GW h 2 about 1.5 times larger than that from the empirical fitted spectrum for this benchmark point.
ALPs into dark photons.We therefore will restrict our use of the empirical fit within this range of α.
In the next section, we will examine our linear order numeric GW spectrum as well as the empirical profile above against the NG15 observation and compare the fitting results.

Calculation setup
To solve the coupled ALP-dark photon equations of motion, at the leading order we replace the friction term on the right-hand side of Eq. ( 2) to its expectation value, i.e.
We discretize the gauge modes v ± (k, τ ) in the wave number k such that we can solve the coupled axion and dark photon equations of motion.The range of k is chosen to be [10 −9 , 10 −6 ] Hz such that the GW signals of interest are covered (note that we have f GW = k/(2π) for frequency today), which we discretize into N = 200 log-uniform modes.The dark photons are assumed to be non-thermal such that all their relic abundance is generated via tachyonic production.We then use the Bunch-Davis vacuum as the initial conditions of the gauge modes, i.e. v ± (k, τ ) = e ikτ / √ 2k.The ALP field, on the other hand, has an initial value of |ϕ 0 | = f and no initial velocity.Without loss of generality, we choose ϕ 0 = −f such that we have the gauge modes v + to be more enhanced.Furthermore, for simplicity, we ignore all the v − modes as their energy density is very suppressed than v + modes [77,79].The differential equations are solved from the time τ min where H(τ = τ min ) > m, until the time τ max = 10τ s .We have tried various choices of k range, gauge mode number N and the ending time of the differential equation solution τ max and found the results to be consistent.
With the GW profiles numerically obtained as above, we run Markov Chain Monte-Carlo (MCMC) with emcee [91] on the parameter region m ∈ [10 −12.5 , 10 −10.5 ] eV, f ∈ [10 17.5 , 10 18.3 ] GeV to get the posterior distribution, with the step width chosen as 0.1 and 0.05 in the exponent of 10 (logarithmic internal), respectively for the first two parameters.Whereas, for α direction we choose an interval of 2 for α ∈ [40,50] and 5 for α ∈ [50,100].The finer scan for the first α range is motivated as the GW production from this system becomes efficient around α ∼ 45.
The prior is chosen to be log-uniform for m and f , and uniform for α.The likelihood is obtained from the NG15 free spectrum density shown in Ref. [1], where we treat the posterior distributions in the violin plot as the likelihood distribution in excess timing delay of the corresponding frequency bin, which can be easily converted to log 10 (Ω GW h 2 ) [44].For each frequency bin, we then have a normalized likelihood function of log 10 (Ω GW h 2 ).We assume there is no correlation between different frequencies, such that the overall likelihood of the GW spectra is a simple multiplication of its likelihood in each frequency bin.We have used GetDist to analyze the MCMC samples and generate the plots.[92]

Results and Constraints Posterior Distributions
In this section, we present the results of the MCMC analysis for the parameter regions favored by the NG15 observation and derive the constraints on this scenario.The favored regions are shown for the linear (blue) and lattice fitted analysis (orange) in Fig. 3 where we plot the posterior distribution triangle plots.We first see that the two analyses differ by an order of magnitude in the 1d posterior of the ALP mass, with masses above 10 −13 − 10 −12 preferred.The ALP mass mostly determines the peak frequency, see Eq. (11), which changes between the two analyses by around an order of magnitude.In contrast, the preferred f regions are very similar -as expected from the similar GW magnitudes in the two analyses (recall that Ω GW ∝ f 4 as derived at the end of Sec. 2).In both cases, the favored f values lie around 8 × 10 17 GeV, almost independently of the other parameters.
We do not find a strong dependence on α, however, two comments are in order.Firstly, for small α ≲ 50 signal from the LO calculation disappears, and so the plot cannot be extended into the small α region.For large α, we cannot trust the linear order (LO) calculation due to the large back reaction, and the empirical lattice fit Eq. ( 12) is only valid for α < 100 as explained in Sec.2.3.The higher α region is therefore the only loophole to our conclusions here.
Comparing with the results of [83,88] for the same model, but with the previous data, we find that the preferred values of m and f by the NG15 data are significantly higher than those preferred by NG12.This can be explained by the inclusion of higher frequency bins in the NG15 analysis [1] compared with the NG12 [84] data, and moreover, the signal strengths in these bins imply GWs of a higher magnitude.As a result, both the peak frequency and the inferred magnitude are higher, resulting in higher m and f (see Eq. ( 11)).Following the same argument we can infer the fitting results on data from the other PTAs, i.e., EPTA, and PPTA.The data used for SGWB analysis in EPTA and PPTA extends to even larger GW frequencies and spectral amplitudes than NG15 and therefore would favor an even larger f and m.See The darker and lighter regions are for the 68% and 95% CL regions.The blue and orange shaded areas correspond to our linear order (LO) numeric GW spectra and the empiric GW spectra for the lattice calculation [83] (empirical lattice fit, ELF), respectively.
Appendix A for some more details.

Cosmological Constraints
We now move to discuss the constraints on this scenario from cosmological observations, mainly arising from the energy density in the ALP and gauge fields.After the tachyonic production stops being efficient, the ALPs start to behave like matter.Given the axion-dark photon coupling, the decay width of ϕ → XX can be shown to be Γ ϕ→XX = α 2 m 3 /(64πf 2 ), implying the ALPs to be cosmologically long-lived in the parameter space of interest.The model is therefore excluded if the dark matter abundance is higher than the observed one (i.e.over-closure) and if the dark radiation exceeds the bounds on ∆N eff .The energy densities of the ALPs and the dark photons are calculated as which we evaluate at the time when the tachyonic production finishes and ALPs and dark photons start to behave like normal matter and radiation, and hence can be redshifted to calculate the fraction of their energy density today.The issue of over-closure was already pointed out in [78,83,93], and we find similarly that the relic abundance of the ALP exceeds that of the observed DM abundance by more than 4 orders of magnitude in all the preferred parameter space.See Fig. 4 where we show the values of the ALP to DM density ratio r for α = 50, 100 in the linear calculation and both the linear and lattice-fitted posterior analyses.In drawing the constant r contours in the figure, we have averaged over ALP energy density from τ max /4 to τ max by properly taking into account the redshift dependence.This results in smoother contours and the averaging period is well beyond the end of the tachyonic band which occurs at τ ≈ τ max /10 signifying the end of particle production.The issue of overproduction of DM may be generic on ALP models if they are to be the major source of nanohertz SGWB [94].There have been some proposals to mitigate the over-closure problem with extra model building that tries to suppress the axion mass (hence energy density) at late times.The models exhibit dilution of the axion energy density at late times which can be achieved via a time-dependent axion potential originating from dark sector phase transition, monopole effects, etc [88,95,96].The NG15 data favors a higher mass of axion, hence a higher T osc ≈ O(10 − 100) MeV.Thus it leaves a wider window in temperature to model away the overdensity problem before the matter-radiation equality at T ≈ 1eV.
We now move to the ∆N eff constraint.The current measurement on the CMB requires ∆N eff ≲ 0.3 [87], and BBN gives a slightly weaker constraint of ∆N eff ≲ 0.4 [86].In Fig. 5 we show the posterior distribution of m and f for α = 50 and 100 together with the contours of ∆N eff .The black ∆N eff contours are obtained from linear order calculation.Similar to the r contours we have employed similar averaging of DR energy densities for smoother contours.We also show the results for the ∆N eff calculation (in dashed lines) from Ref. [83,88], where the authors assumed instantaneous deposition of the total energy from ALP to DR at scale factor a * where ρ ALP ≈ ρ DR .This approximation gives similar results to our LO ∆N eff calculation for small α.For higher values of α, in our linear calculation, not all the ALP energy is transferred to the dark photons at a = a * , and a substantial portion is released gradually in the later evolution.∆N eff is then expected to increase by this process as the residual Ω ALP keeps growing (∝ a) untill it decays to Ω DR .We find that our linear result is higher than the one from the instantaneous assumption by a factor of 2 − 2.5, as we can see on the right side of Fig. 5.We stress that for these values of α, the lattice calculation can significantly alter the dynamics and we do not know whether this effect persists in the lattice.Therefore, we will show both results and treat the instantaneous assumption as a conservative estimate of the ∆N eff as it can only underestimate it.A further lattice study is required to confirm our results here.We see that for the preferred region, in both the linear and the lattice fitted analyses, ∆N eff ≳ 1−3, significantly above the current constraint.
Another way to understand the ∆N eff constraint is to isolate the contributions from the   The plot range is chosen to focus on the favored region of the linear order calculations, and a more complete distribution of the empirical spectra can be found in Fig. 5.The ratio r is obtained from the linear order calculations.The NG15 favored mass regions have the axion rolling and tachyonic production before BBN, as can be seen from the axion oscillation temperature T osc (the cosmic temperature where H = m) shown on the top axes.In both panels, the NG15 preferred parameter space will generate too large an r and thus overclose the universe, and we have checked the situation to be the same for the other values of α.
background DR and GW spectrum.For our linear calculations, we find that across all the viable parameter spaces, the ∆N eff contribution from GWs (∆N eff,GW ) is related to the ∆N eff of DR (∆N eff,DR ) by ∆N eff,DR (∆N eff,GW ) where ∆N eff,GW is calculated as (dΩ GW /d ln k)(dk/k).The value inside the parenthesis corresponds to the instantaneous assumption with α close to 100.This approximate scaling relation holds for all parameter points and for the α range [50,100] used in this paper.This is because the GW spectrum is the two-point power spectrum for the tensor fluctuation in the DR and is thus related to the square of DR energy density.For the GW power spectrum to explain the NG15 observation (like the benchmark on Fig. 1), it would require ∆N eff,GW ≈ 0.003.This in turn gives ∆N eff,DR ≈ 1 − 4 using Eq. ( 15) and thus rules out the parameter point from the cosmological ∆N eff measurement.Conventional mechanisms that dilute ∆N eff through entropy injection will not work since that will also dilute the GW signal by the same fraction.From Eq. (15) it naively seems that one can decrease their ratio ∆N eff,DR /∆N eff,GW by increasing both ρ GW and ρ DR while holding ρ DR /ρ 1/2 GW fixed at the end of tachyonic production.However, any such significant increase would result in the dark sector dominating the energy density of the universe, where we don't expect the scaling in Eq. ( 15) to hold.The analysis of this scenario is beyond the scope of this work.Given the tight cosmological constraints from the dark matter relic abundance and ∆N eff , one may inspect the maximum SGWB amplitude that could be generated by the tachyonic instability mechanism.If the over-abundant axions can be depleted with additional model building, we will need to decrease ∆N eff,DR by an order of magnitude to be compatible with the current CMB measurements, which renders the GW amplitude to reduce by about two orders of magnitude, making up about 1% of the measured NG15 signal.However, if it turns out that the over-closure bound cannot be modeled away, the current setup will imply the induced GW can only correspond to a minuscule < O(10 −8 − 10 −12 ) fraction of the observed GW spectrum depending on the α values.

Conclusion and Discussion
In this work, we explore the possibility that the reported SGWB found in the NANOGrav 15-year data set [1] is generated from the coupled dynamics between the ALPs and the dark photons.The ALPs' rolling exponentially amplifies the massless gauge modes, and the drastic change of the gauge boson mode functions then sources GWs which can potentially explain the observed SGWB spectrum.However, as we have shown in this paper, the parameter region that can explain the NG15 observation is incompatible with the observed dark matter relic abundance in the current universe, as well as the ∆N eff constraints from BBN and CMB.
We have additionally checked the robustness of our results.We have carried out numerical calculations at the linear order, and have compared with the results from lattice-fitted formulae from [83].We find that the results remain independent.We do not, however, explore the large back-reaction parameter region with α > 100 for which both analyses are invalid and new lattice calculations are required, so our argument may not apply in that region.A naive extension of our fit into this region, however, would only strengthen the constraints due to the positive correlation between α and f in the GW peak position.
When generating the MCMC chains, we use the posterior "violins" of the NG15 free spectrum density as the likelihood function of the sample points.It is known that this statistics approach usually results in a slightly smaller favored GW amplitude compared with going through the whole pipeline, as it ignores the correlation between different frequency components [97].The tension between the instability-induced GW explanation of NG15 and the cosmological constraints from dark matter relic abundance and ∆N eff could therefore be even stronger than what we have pointed out in the main text.Thus, a more sophisticated analysis most likely will strengthen the bounds and will not change the conclusion of the paper.
The stringent constraints on the vanilla ALP scenario warrant additional model building to address the over-closure and ∆N eff constraint.Although several mechanisms exist to address the ALP over-closure, reducing ∆N eff without comprising the GW signal is highly challenging.Thus the NG15 observation in conjunction with other cosmological measurements strongly rules out the instability-induced GW from the vanilla ALP model to be the entire signal in the m ≈ 10 −14 − 10 −11 eV ALP mass range.It can still potentially be an at most percent-level subdominant signal if another source of GW, such as supermassive black hole binaries, is responsible for the NG15 signal, and it may be of particular interest to study how such a signal can be potentially probed.Figure 6: Left: Our fitted result to the EPTA released HD free spectrum [98] for fixed α = 100, companion to the right panels of Fig. 4 and 5. Right: Comparison between the posterior "violins" used for the two fittings.
of 28 frequency bins are involved, which extends to even larger frequency and amplitude [4].Therefore, following the same argument, PPTA cannot improve our situation either.

Figure 1 :
Figure 1: An illustration of the tachyonic production mechanism and the induced GWs from the linear-order calculation, where the model parameters are set to m ≈ 8 × 10 −12 eV, f ≈ 8 × 10 17GeV, α = 50.Left: The evolution of the ALP (blue), dark photon (orange), and GW (red) energy density, normalized by the total energy density of the standard ΛCDM cosmology.The dashed green line is proportional to a, showing that ALPs now evolve like matter.Soon after the ALP starts rolling, at a/a osc ≈ 20 the tachyonic production of gauge modes is maximized and the ALP quickly dumps its energy into dark radiation.After the tachyonic production, the ALP and dark photon continue to behave like ordinary matter and radiation.On the other hand, this set of parameters has the ALP relic abundance overclose the universe and has too many radiation degrees of freedom at the BBN and CMB.Right: The GW generated from the tachyonic instability (red curve), compared with the first 14 frequency bins of NG15's free spectrum density (gray shades), taken from Ref.[1].

Figure 2 :
Figure2: A comparison between the benchmark GW spectrum shown in Fig.1(red) and the empirical GW profile (black).The model parameters are set to be m ≈ 8×10 −12 eV, f ≈ 8×1017  GeV, α = 50, the same as those in Fig.1.The linear calculation gives an integrated Ω GW h 2 about 1.5 times larger than that from the empirical fitted spectrum for this benchmark point.

Figure 3 :
Figure3: The posterior distribution for the parameters of the instability-induced GW model for NG15 observation.The darker and lighter regions are for the 68% and 95% CL regions.The blue and orange shaded areas correspond to our linear order (LO) numeric GW spectra and the empiric GW spectra for the lattice calculation[83] (empirical lattice fit, ELF), respectively.

Figure 4 :
Figure 4: The r ≡ ϕ /ρ DM contours of various values shown on top of the posterior distribution of fixed α, where the α values are chosen to be 50 and 100 for the left and right panel, respectively.The plot range is chosen to focus on the favored region of the linear order calculations, and a more complete distribution of the empirical spectra can be found in Fig.5.The ratio r is obtained from the linear order calculations.The NG15 favored mass regions have the axion rolling and tachyonic production before BBN, as can be seen from the axion oscillation temperature T osc (the cosmic temperature where H = m) shown on the top axes.In both panels, the NG15 preferred parameter space will generate too large an r and thus overclose the universe, and we have checked the situation to be the same for the other values of α.

Figure 5 :
Figure 5: Similar to Fig.4, but the solid contours are for ∆N eff obtained in the linear calculations.The plot range is zoomed out to include a more complete parameter space of the empirical spectral.∆N in the plot is short for ∆N eff .The dashed gray lines are the ∆N eff contours obtained by assuming instantaneous energy transfer from the ALP to the dark photons[88].Models with parameters favored by NG15 severely violate the current constraints on ∆N eff .