NANOGrav signal from axion inflation

Several pulsar timing arrays including NANOGrav, EPTA, PPTA, and CPTA have recently reported the observation of a stochastic background of gravitational wave spectrum in the nano-Hz frequencies. An inflationary interpretation of this observation is challenging from various aspects. We report that such a signal can arise from the Chern-Simons coupling in axion inflation, where a pseudoscalar inflaton couples to a (massive) $U(1)$ gauge field, leading to efficient production of a transverse gauge mode. Such tachyonic particle production during inflation exponentially enhances the primordial perturbations and leads to a unique parity-violating gravitational wave spectrum, that remains flat near the CMB scales but becomes blue-tilted at smaller scales. We identify the parameter space consistent with various cosmological constraints and show that the resultant gravitational wave signals can provide extra contribution on top of the standard astrophysical contribution from inspiraling supermassive black hole binaries towards explaining the observed excess at NANOGrav. The parity-violating nature of the signal can be probed in future interferometers, distinguishing it from most other new physics signals attempting to explain the NANOGrav result.

In this paper we focus on a possible origin of this signal from an axion inflation model [75][76][77][78][79][80][81][82][83][84][85][86][87][88], in which the pseudoscalar inflaton ϕ with an approximate shift symmetry [82] has a Chern-Simons coupling ϕF F to a U (1) gauge field, where F is the field strength of the gauge field and F is its dual.This coupling enables a tachyonic production of a transverse mode of the gauge field [89][90][91][92][93][94][95][96][97][98][99][100][101][102][103][104], which enhances the primordial scalar and tensor perturbations.On one hand, it introduces significant non-Gaussianity in the scalar perturbations, on the other hand, it yields a characteristic parity-violating GW spectrum that remains flat at low CMB frequencies but rises at smaller scales probed by the gravitational wave interferometers.Intriguingly, the backreaction of the gauge fields on the inflationary dynamics tames the growth of the GW signal, helping it evade stringent constraints from cosmological considerations and non-observation of SGWB at LIGO scales [105].The degree of backreaction, which determines at what frequency the signal starts to rise, depends on the inflationary potential.Here we consider the T-model [106] as a specific example currently favored by Planck 2018 data [107], and show that the GW signal becomes blue-tilted at small enough frequencies to reach the amplitude range reported by the PTAs, above the expected contribution from inspiraling supermassive black hole binaries (SMBHBs), while remaining consistent with various cosmological constraints.
The production of massive gauge fields during inflation has at least two more important aspects.First, it generates a parity-violating GW spectrum coming from two polarizations of the graviton, a unique feature that distinguishes it from most other sources of SGWB.Such parityviolating nature can be probed in future interferometers [108,109].Second, it has interesting implications for "cosmological collider" physics , where the mass and spin of a particle produced during inflation leave imprints on the oscillatory three-point correlation function of the scalar perturbations in the "squeezed limit".The fact that the same process generates gravitational wave signals potentially observed at PTAs offers a remarkable opportunity to probe cosmological collider physics in a scale well beyond the CMB/LSS scales.
The paper is organized as follows.In section 2, we review the phenomenology of gauge field production in axion inflation and discuss relevant constraints on the model.We discuss the evolution of the inflationary variables from CMB scales to smaller scales in the context of the T-model potential in section 3.In section 4, we discuss the gravitational wave signals in this model and show that the results provide extra contribution on top of SMBHBs towards explaining the NANOGrav 15-yr signal at face value.We conclude in section 5.

Phenomenology of gauge boson production
We contemplate a single field slow-roll inflation scenario, where the inflaton ϕ is an axion-like pseudoscalar with an approximate shift symmetry ϕ → ϕ + const.Due to the shift symmetry, the inflaton couples to a U (1) gauge field A µ via a Chern-Simons term, −ϕF F /(4Λ), where is the field strength of the gauge field, and The action can be written as where ϵ 0123 = +1 is a tensor antisymmetric in any two indices.We assume a quasi-de Sitter space with a slowly varying Hubble rate H, with the scale factor given by a(t) = e Ht .The metric can be expressed as where t and τ are the physical and conformal time, respectively.We assume that the inflaton field has a homogeneous background, on top of which we consider an inhomogeneous perturbation, ϕ(t, x) = ϕ 0 (t) + δϕ(t, x).This allows us to study the background dynamics and the correlation functions of the perturbations separately.
The dynamics of the background inflaton field can be expressed in terms of the following coupled equations where the physical 'electric' and 'magnetic' fields corresponding to the gauge field are given by E = −A ′ /a 2 and B = (∇ × A)/a 2 .The terms on the right hand side act as sources of backreaction of the gauge fields on the inflationary dynamics, and can have a significant impact on observables at and beyond the CMB scale.Backreaction effects can no longer be neglected when the source terms on the right hand side are comparable to the terms on the left hand side.
For the general case of massive gauge bosons, we employ the constraint ∂ µ ( √ −gA µ ) = 0 derived from the equation of motion and decompose the gauge field in the helicity basis λ = ±, 0, (2.5) where the polarization vector ϵ λ (k) and the annihilation operator a λ (k) obey the usual commutation and orthonormality relations.A λ (τ, k) is the mode function, where the longitudinal mode and two transverse modes are denoted by λ = 0 and λ = ±, respectively.
The dominant vector field production is governed by the field equations of the transverse modes, where ξ ≡ φ0 /(2ΛH) and we have used a ≈ −1/(Hτ ) during inflation.Taking φ0 > 0 without loss of generality, only the A + mode experiences tachyonic instability and is dominantly produced.We assume the usual Bunch-Davies initial condition and treat the inflaton's rolling speed φ0 to be a constant during inflation.The solution of the mode function can then be written in terms of the Whittaker-W function, where µ ≡ (m A /H) 2 − 1/4.The A − mode is exponentially suppressed, and the longitudinal mode is unaffected by this production mechanism, but can otherwise be produced purely from gravitational effects [147][148][149][150][151]. The mode function solution for the massless case is derived by setting m A = 0 in eq.(2.7).For a particular k mode, the energy density of the gauge mode function A + receives an exponential enhancement for ξ > m A /H when −kτ ∼ O(1), which overcomes the Boltzmann suppression and enhances particle production, so that the physical observables impacted by the gauge boson production have an approximately exponential dependence on ξ − m A /H [75].We treat ξ and m A /H as free parameters of the model that are later constrained from various phenomenological considerations.
We now turn to the effect of gauge field production on primordial scalar and tensor perturbations.The inflaton perturbations follow the equation of motion [75] where β ≡ 1−2πξ⟨E • B⟩/(3ΛH φ0 ).From the inflaton perturbations, the curvature perturbation on uniform density hypersurfaces is defined as To calculate the tensor perturbation, we write the perturbed metric in terms of the tensor perturbation h ij using the scalar-vector-tensor decomposition as where h ij is transverse (∂ i h ij = 0) and traceless (h ii = 0).The equation of motion of h ij is given by [152] h where M Pl ≃ 2.4 × 10 18 GeV is the reduced Planck mass, and T T T ij is the transverse and traceless part of the stress-energy tensor.We decompose the tensor perturbation into two helicity modes 12) The correlation functions of the curvature and tensor perturbations are calculated at τ 0 = 0 after the end of inflation.For details of the calculation using the in-in formalism [153], we refer the interested readers to refs.[75,76].Here we focus on the phenomenological observables derived from the two and three-point correlation functions of the primordial perturbations.The scalar power spectrum is proportional to the two-point correlation function of the curvature perturbation and can be written as (2.13) comes from the impact of the gauge field production.Here ′ denotes that the δ-function ) is stripped off.The amplitude of the scalar power spectrum at the CMB scale is measured to be P ζ ≃ 2.4 × 10 −9 [154,155], which accounts for the contribution of the vacuum modes as well as the extra degrees of freedom (gauge modes in this case).Conservatively taking the gauge field's contribution to be subdominant at the CMB, we can ignore P [A] ζ and fix This assumption would be valid as long as Gauge field production during inflation may introduce significant non-Gaussianity in the curvature perturbations.It can be parametrized by a dimensionless quantity, where the superscript 'eq' represents the equilateral shape k 1 = k 2 = k 3 ≡ k of the three-point correlation function.Current bound on equilateral non-Gaussianity gives f eq NL = −25 ± 47 at 68% CL [156].
Another important constraint comes from the ratio of the tensor to scalar power, r ≡ P h /P ζ at the CMB scales, where the tensor power spectrum is chiral and can be decomposed as where ± corresponds to the two polarizations of the graviton.Here the vacuum contribution has been equally distributed to the two polarizations.In our case P + h ≫ P − h , implying a parity violating tensor power spectrum.At the CMB scale, tensor-to-scalar ratio is constrained by current data r * ≤ 0.056 at 95% CL [107].
In fig. 1, we show the parameter space violating the above constraints, along with the region where the assumptions of negligible backreaction and negligible contribution of the gauge modes to the scalar power spectrum at CMB scales breaks down.
Evidently, the upper bound on scalar non-Gaussianity puts the strictest constraint on the parameter space of the model at this stage.It implies that ξ is restricted not too far from m A /H at the CMB scale.However, as we will show in the next section, even for ξ − m A /H ∼ O(1) at the CMB scales, the power spectrum of tensor perturbations can evolve to generate observable GW spectrum at PTA and other interferometer scales.

Evolution beyond the CMB scale
Modes responsible for CMB scale observables can be assumed to experience constant ξ and H as long as slow-roll conditions prevail.However, modes contributing to the GW spectrum at smaller scales experience the time evolution of these parameters and are subject to strong backreaction from the inverse decay of the gauge field.These effects can be studied from a simultaneous solution of the coupled eqs.(2.3) and (2.4), ignoring the source term in eq.(2.4) as it is negligible compared to the source term in eq.(2.3).For convenience, we change variables from time t to the e-folding number N , defined as dN = −Hdt, where N decreases as we  approach the end of inflation.Eqs.(2.3) and (2.4) can then be expressed as ) Solving eqs.(3.1) and (3.2) numerically for a given potential V , we get H(N ) and ϕ(N ), which can be used to calculate ξ(N ) ≡ 1 2Λ dϕ dN .We now specialize to the example of the T-model potential, which yields a consistent combination of the spectral index, n s , and the tensor-to-scalar ratio, r with respect to the combined Planck 2018 data analysis [107].The inflaton potential in this model is given by where V 0 and α T are free parameters, which can be constrained from CMB measurements of n s = 0.9649 ± 0.0042 (at 68% CL) and r < 0.056 (at 95% CL) [107].We choose α T = 4 and V 0 ≈ 1.7×10 −9 which are consistent with all of the above constraints.The shape of the potential is shown in fig. 2. We illustrate the evolution of ξ and m A /H as a function of N in fig. 3 for three benchmark points for this potential.
At N ∼ 60 when the CMB modes leave the horizon, standard slow-roll condition prevails and backreaction effects can be neglected.Furthermore, the selected benchmark points are in the allowed region in fig. 1, so that large non-Gaussianities can be avoided.Beyond that, ξ increases rapidly for about 20 e-folds after which backreaction effects slow down its rise.Closer to the end of inflation, slow-roll condition is re-established, and ξ rises rapidly.In this regime, however, the Hubble rate also rises swiftly, compensating the parameter m A /H with respect to ξ, so that particle production remains under control.

Gravitational wave signatures
We now turn to the gravitational wave spectrum of the model potentially observed at NANOGrav and within the sensitivity of various upcoming interferometers.The tensor perturbations sourced by the gauge modes leave the horizon during inflation and can source gravitational waves after re-entering the horizon at later stage.Taking the redshift into account, the amplitude of the GW signal observed today is given by Here Ω R,0 ≃ 8.6 × 10 −5 denotes the radiation energy density today and P h (f ) is the frequency dependent power spectrum of the tensor fluctuations at the time of horizon exit.Since P + h ≫ P − h , we have a parity violating GW spectrum.
The power spectrum depends on the model parameters ξ and m A /H, whose time evolution was discussed in sec.3. We can relate the time parameter N to frequency of the GW signal f as [100] where the CMB pivot scale is k CMB = 0.002 Mpc −1 and we take N CMB = 60.
Here we specifically focus on the parameter space that can potentially explain the observed excess at the NANOGrav 15-yr data.The effect of the gauge field creation on the tensor fluctuations is minimal for CMB scales and the power spectrum is dominated by the vacuum fluctuations.Current bound on scale-invariant stochastic gravitational wave at the CMB scales [157] implies a tensor-to-scalar ratio r < 0.056 [107], which gives H/M Pl ≲ 2.6 × 10 −5 , and Ω GW < 1.2 × 10 −16 .Larger frequencies correspond to modes which left the horizon later than the CMB modes.By that time the rolling speed of the inflaton increases and the Hubble rate decreases, the combined effect of which implies a larger value of ξ.This dramatically enhances the power spectrum of the tensor perturbations sourced by the gauge field and it quickly supersedes the contribution from the vacuum fluctuations.Gravitational wave amplitude that eludes observation at the CMB scale now offers the possibility of detection at the interferometer scales.
Although the blue-tilting of the tensor power spectrum at higher frequencies is a known feature of the axion inflation scenario, achieving a signal at the nano-Hz range probed by PTAs is significantly challenging.It requires the signal to start to rise at sufficiently smaller frequencies, yet such signals should not violate the upper bound at the LIGO-VIRGO-KAGRA (LVK) [105,158,159].Furthermore, such a wideband signal with a blue-tilted spectrum is likely to introduce significant contribution to ∆N eff .These characteristics depend essentially on the inflaton potential one considers, as the rolling speed of the inflaton is determined by the potential.Intriguingly, we find that such peculiar features can be nicely accommodated in the T-model potential, but are difficult to achieve in the broader class of α-attractor models [160][161][162][163][164][165].
We show the posterior distributions of the model parameters m A /H CMB and ξ CMB − m A /H CMB with dashed lines, and 68% (darker blue) and 95% (lighter blue) Bayesian credible regions for the NANOGrav 15yr results using PTArcade [166] in fig. 4. The 2D posterior regions shrink for larger m A /H CMB as backreaction from gauge fields becomes so strong for larger ξ CMB that the amplitude of the GW signal at NANOGrav frequencies is diminished.In the same plot we also show the constraints from f NL in the equilateral limit (gray region) and excess GW amplitude at LVK scales (orange region).Interestingly, these constraints completely rule out the parameter space for m A /H ≲ 0.6.We show the GW spectrum for three benchmark points from the viable parameter space in fig. 5.In each case we keep the parameter m A /H same at the CMB scales but let the parameter ξ vary.For comparison we show the upper limits from CMB [157] and LVK [105,158,159], and the projected sensitivites of future detectors SKA [167], µ-Ares [168], LISA [169], DECIGO [170,171], BBO [172], AEDGE [173], AION [174], CE [175], ET [176] and future upgrades of LVK.The recent results from NANOGrav and EPTA are shown with blue and orange violins, and the expected astrophysical background from inspiraling supermassive black hole binaries (SMBHBs) [14,15] is shown in the PTA frequencies with a black dashed line.The GW spectrum from SMBHB best fit parameters fall short both in amplitude and slope in fully explaining the power law fit of the observed GW background in the NANOGrav data, which could be improved if a GW signal with larger amplitude and steeper slope were present [14].All the three benchmark signals originating from our model satisfy this criteria.This leads to an improved fit to the spectral characteristics of the NANOGrav signal.For all three benchmark points, the signal remains flat at smaller frequencies as the contribution from gauge boson production remains subdominant compared to the vacuum fluctuations.The signals begin to rise at f ≳ 10 −13 Hz, but the growth is slowed down near the PTA frequencies.For higher frequencies, backreaction effects saturate the signals so that they remain below the upper bound at LVK frequencies.These signals can be detected in various terrestrial and space-based interferometers from nano-Hz to kilo-Hz range.
While wideband stochastic GW signals can appear from some other sources, for example, topological defects like cosmic strings, an important distinction of the signals produced in the present case is their parity-violating nature.Such parity-violating nature is expected to be detected in the ET-CE [109] and the LISA-Taiji network [108].
We now calculate the contribution of the GW spectrum to the radiation energy budget of the Universe in terms of the effective number of relativistic species N eff .The extra contribution is given by where f min and f max depend on the era of interest and the maximum temperature reached.One needs to make sure that the GW spectrum does not spoil the successful prediction of BBN, corresponding to the frequency of a mode crossing the horizon during the time of BBN when the temperature of the radiation bath was T ∼ O(MeV).We take f max = 10 4 Hz for a sufficiently large reheating temperature.The current upper bound on ∆N eff from BBN and CMB probes is ∆N eff ≃ 0.4 [177].We find that all the benchmark points satisfy this bound, with ∆N eff ∼ 0.013 in all three cases.In fact, we can do a naive estimate to show that any signal that does not violate the LVK bound would also not violate the ∆N eff bound in our model.Since the spectrum does not have any peak, suppose we take it to be flat at Ω GW h 2 ∼ 10 −8 for all relevant frequencies.Even then, we find ∆N eff ≲ 0.1, well below the upper bound from BBN and CMB probes.Both the tensor and scalar perturbations are enhanced by gauge fields.Large scalar perturbations can lead to the creation of primordial black holes (PBH).The mass of a PBH is related to the e-folding number N when a fluctuation responsible for the creation of the PBH leaves the horizon.An upper bound on the scalar power spectrum as a function of N has been derived in ref. [97] using the estimates of refs.[178,179], as shown in fig.6 with a dashed curve.In the same plot we show the scalar power spectrum for the three benchmark points.In all cases the power spectrum suffers from strong backreaction at smaller scales.Nevertheless, we see some tension of the spectrum with the PBH bound, which is within the theoretical uncertainties.Furthermore, recent lattice simulations [180] indicate that the combined effect of non-Gaussian perturbations become more Gaussian overall at smaller scales, which weakens the bound (shown with a dot-dashed line in fig.6) and allows the model to avoid overproduction of PBHs.The problem can be further alleviated by considering N copies of the gauge field, which reduces the scalar power spectrum by a factor of N [90], but does not significantly affect the tensor power spectrum at smaller scales.Large scalar perturbations can also generate second order tensor perturbations [181], which will be discussed in a future work.

Conclusion and outlook
We have pointed out that an explanation of the observed excess in the NANOGrav 15-yr dataset is possible from the gravitational waves generated in axion inflation from axion coupling to (massive) gauge bosons.Such a coupling, natural for a shift-symmetric inflaton, leads to copious particle production during inflation, leaving an indelible imprint on the primordial scalar and tensor perturbations.This leads to a unique parity-violating contribution to the GW spectrum, that remains flat at CMB scales, but blue-tilts at smaller scales and can become audible to pulsar timing arrays.The growth of the spectrum at higher frequencies, potentially dangerous for LVK scales, is controlled by the backreaction of the gauge quanta on the inflationary dynamics, the details of which depend somewhat on the inflaton potential.We have specifically discussed the case of the T-model potential, a model currently favored by data, and have identified the parameter space that can potentially produce a GW spectrum statistically consistent with the NANOGrav result.More specifically, we have shown three benchmark points for which the amplitude and spectral slope of the GW signals provide better fit to the NANOGrav result than the standard astrophysical background from inspiraling SMBHBs, suggesting an interpretation of the NANOGrav excess from GW waves generated in axion inflation is quite likely.The GW spectrum can also be probed in various detectors from nano-Hz to kilo-Hz frequencies, while its parity-violating nature would clearly distinguish it from other wideband signals.

Figure 1 .
Figure 1.Shaded regions denote exclusion of the gauge boson's parameter space from various constraints.Tensor-to-scalar ratio bound is drawn for H/M Pl = 10 −5 .For f eq NL < −25 ± 47 [156], the left (right) part corresponds to the positive (negative) bound.

Figure 3 .
Figure 3. Evolution of model parameters ξ and m A /H for three benchmark points in the context of the T-model inflaton potential.

Figure 4 .
Figure 4. Posterior distributions for the axion inflation model parameters m A /H CMB and ξ CMB − m A /H CMB .The diagonal subplots represent 1D marginalized distributions with vertical lines showing 68% and 95% Bayesian credible intervals.In the lower left subplot, the darker and lighter blue regions show the 68% and 95% Bayesian credible regions in the 2D posterior distributions.The gray region is excluded by the upper bound on equilateral f NL .The orange region shows the parameter space excluded by the upper bound on GW amplitude at LVK.

Figure 5 .
Figure 5. Gravitational wave spectrum for gauge boson production in the context of the T-model potential.

Figure 6 .
Figure 6.Scalar power spectrum for the T-model potential.Gray regions represent upper bounds from the overproduction of primordial black holes.