Whispers from the dark side: Confronting light new physics with NANOGrav data

The NANOGrav collaboration has recently observed first evidence of a gravitational wave background (GWB) in pulsar timing data. Here we explore the possibility that this GWB is due to new physics, and show that the signal can be well fit also with peaked spectra like the ones expected from phase transitions (PTs) or from the dynamics of axion like particles (ALPs) in the early universe. We find that a good fit to the data is obtained for a very strong PT at temperatures around 1 MeV to 10 MeV. For the ALP explanation the best fit is obtained for a decay constant of $F \approx 5\times 10^{17}$ GeV and an axion mass of $2\times 10^{-13}$ eV. We also illustrate the ability of PTAs to constrain the parameter space of these models, and obtain limits which are already comparable to other cosmological bounds.


I Introduction
With the first direct observation of gravitational waves (GWs) by LIGO [1], a new era in astrophysics and cosmology has started. Since GWs travel almost undisturbed through spacetime, they can carry information from before the time of CMB emission, which is where our direct observations using electromagnetic radiation end. GWs therefore open a new window to the early Universe.
Pulsar Timing Arrays such as EPTA [2], PPTA [3] and NANOGrav [4] are sensitive to GWs with frequencies of 10 −8 Hz and below. A stochastic background of such low frequency GWs could be produced in the early universe by a variety of processes, such as inflation, cosmic strings, phase transitions, or scalar field dynamics [5]. The most recent data release of the NANOGrav collaboration [6] for the first time shows evidence for such a stochastic GW background, which is well described by a f −2/3 power law spectrum with a GW strain amplitude of 2×10 −15 , or equivalently a GW energy density Ω GW h 2 of order 10 −10 . This is indeed consistent with the GW density one expects from a variety of cosmological sources, as was discussed for the case of cosmic strings [7][8][9], phase transitions [10,11], or primordial black hole formation [12,13].
So far these studies have focussed on demonstrating that a sufficiently large GW density can be achieved in these models in the required frequency range. Here we perform the first fit to the frequency binned NANOGrav data. Since most cosmological sources of GWs have specific spectral features, it is important to verify that indeed they agree well with the data. In doing this, we are able to obtain best fit parameter regions for two classes of models that produce primordial GWs, namely phase transitions in the early universe [14][15][16][17][18] and audible axions [19][20][21]. We also show that the NANOGrav data already puts constraints on the parameter space of these models, which are comparable to the ones coming from other astrophysical observations such as big bang nucleosynthesis (BBN) or the constraint on the number of relativistic degrees of freedom, N eff .
With more precise data it will become possible to distinguish between different cosmological sources and from the expected background due to supermassive black hole binaries. Our work presents a first step in this direction. It is organised as follows: In the next section, we describe our effort at recasting the NANOGrav data, and re-derive the best fit regions for single power law fits. The following two sections introduce the parameterisation of the stochastic GW background produced by audible axions and phase transitions, respectively, and the best fit regions for the model parameters, before we present our conclusions.

II Refitting the NANOGrav data
The magnitude of a stochastic GW background is typically described by the dimensionless, frequency dependent characteristic strain amplitude h c (f ). For a single power law it can be written as where A GW is the amplitude, α is the slope and f y = 1/year is a reference frequency at which the amplitude is fixed. An important related quantity is the energy density in GWs as a fraction of the critical energy density, Ω GW , which is given by [4] Ω where H 100 = 100 km/s/Mpc and H 0 = h H 100 is the Hubble rate today with h ≈ 0.7. In Fig. 1 of Ref. [6] the NANOGrav collaboration provides the results of different fits to the data, namely a free spectrum fit of the individual frequency bins, a fit of a single power law to the lowest 5 frequency bins or to all 30 bins, and a broken power law with different slopes for the low and high frequency part of the data. The high frequency bins are expected to be dominated by white noise with slope α = 3/2, which is corroborated by the broken power law fit. Instead the 5 lowest frequency bins contribute 99.98% of the significance of the potential GW signal.
In the following, we will therefore fit our signal models to the 5 lowest frequency bins, assuming that the remaining data points are explained by white noise. The results of the free spectrum fit are given in terms of the timing residual, which is related to the characteristic strain as in units of seconds. Note that we have chosen the prefactor in this formula such that by fitting a single power law to the data, we can reproduce the best fit contours of [6], see Fig. 1. In the following sections, we will fit this data with signal templates motivated by concrete new physics scenarios.

III Audible axions and NANOGrav
The audible axion is a simplified model where an axion-like particle a couples to a dark photon X through a term of the form where F is the axion decay constant, i.e. the scale where the global symmetry in the UV is broken and gives rise to the light pseudoscalar a, q is a dimensionless charge, and X µν andX µν are the dark photon field strength tensor and its dual. The axion has a potential V (a) = m 2 a F 2 (1 − cos(a/F )), such that its mass is given by m a . As usual in the axion misalignment mechanism, we assume that after the end of inflation, the axion is displaced from the minimum of V (a) by θF , with θ an order one angle. The axion remains displaced until the Hubble rate becomes of order m a , at which point it starts to oscillate around the origin. It was shown in [22] that the presence of a dark photon leads to a suppression of the axion dark matter abundance, making larger values of F consistent with observations. An efficient energy transfer to the dark photons is possible due to a tachyonic instability that develops while the axion rolls. The same process also amplifies quantum fluctuations in the dark photon field, which grow to macroscopic scales and source a detectable GW background [19].
The GW spectrum produced by audible axions is peaked at the frequency corresponding to the dark photon momentum mode that grows the fastest, and is closely related to the axions mass m a . In terms of the model parameters, the peak frequency, redshifted to today, can be estimated as (III. 2) The amplitude of the GW signal is determined by the strength of the source, i.e. the energy that is initially carried by the axion. This is mostly influenced by the size of the decay constant F . The peak amplitude of the signal can be estimated as To perform our fits we use the signal shape provided in [20] Ω 0 In Fig. 2 we show on the left the best fit of an audible axion compared to the five first frequency bins from NANOGrav. On the right we show the one and two sigma contours in the F -m a plane with θ = 1 and q = 50 fixed. To get such a strong signal the energy in the axion that is transmitted to the dark photon has to be quite significant. The dark photon is a form of dark radiation and therefore contributes to the number of relativistic Values of F and m a which lie above the green contours predict a GW signal which is too large, i.e. this region is excluded by the NANOGrav data. While the N eff is slightly stronger, it is worth noting that PTAs are already able to put competitive bounds on this scenario.

IV Phase transitions and NANOGrav
It has been known for many years that a cosmological phase transition (PT), such as from the spontaneous breaking of a global or gauge symmetry through a scalar field that acquires a vacuum expectation value, produces a stochastic GW background if the transition is strongly first order [14][15][16]. While a large variety of models exists that predict such a transition at different scales, the GW signal of a strong first order PT is universally described by only four parameters, the ratio between the vacuum and total energy density α = ρ vac /ρ tot , the time scale of the transition β/H, where H is the Hubble scale at the time of the transition, the temperature T * at which the transition takes place and the bubble wall velocity v w [17,23].
We use the signal templates in terms of these parameters as given in [24]. The peak frequencies and amplitudes of the two most important contributions to the signal scale as  where n = 1 for the sound wave contribution and n = 2 for the scalar field contribution, and we neglect order one numbers which are not relevant for the qualitative discussion. Very strong transitions are characterised by α > 0.1 and a wall speed approaching the speed of light, v w → 1. The NANOGrav signal corresponds to an energy density Ω GW h 2 > 10 −10 at a frequency around 10 −8 Hz, so that only a strong transition will be able to explain the data. Furthermore we immediately see that T * should be of order 10 −3 − 10 −2 GeV, i.e. the PT should happen at a very low scale. The implications of this for concrete models will be discussed in more detail below.
We consider two scenarios. If the PT takes place at a temperature significantly below the critical temperature, the Universe will be dominated by vacuum energy, i.e. the α dependence drops out of Eq. (IV.2). In such a supercooled PT, no friction acts on the bubble wall, so that v w = 1. Furthermore in the absence of a plasma, the only source of GWs is the scalar field itself, i.e. n = 2 in Eq. (IV.2). In that case, a good fit to the data requires relatively small values of β/H 50, and transition temperatures around or below the MeV scale, as shown in Fig. 3. Above the peak frequency, the GW strain amplitude of the PT signal falls as f −3/2 . Therefore if the peak frequency lies below the lowest frequency probed by NANOGrav, the signal will look like a single power law to the detector. This explains the flat direction in the fit towards lower temperatures and lower values of β/H. However lower values of β/H are increasingly difficult to obtain in realistic models, therefore this region should be considered less favoured.
If the PT is very strong but not supercooled, the bubble walls will still reach a relativistic terminal velocity, so for simplicity we again set v w = 1. In this case sound waves in the plasma induced by the PT are the dominant source of GWs, and the amplitude is only suppressed by one power of β/H. As expected, in Fig. 3 we see that a good fit to the data in the T * − α plane is found both for β/H = 10 and β/H = 100, where in the second case the suppression of the signal is compensated by a larger energy budget α. Again we also find a flat direction, where the peak of the PT signal is shifted below the NANOGrav frequency range, and data is fit by the high frequency tail.
In both scenarios, we find that the PT should happen at a temperature around 1 MeV, with only a small viable region slightly above 10 MeV. Since extensions of the SM at such low scales are almost impossible to hide from laboratory experiments, it is clear that the PT should take place in a dark sector, with only very weak interactions with the SM [24][25][26][27][28][29][30][31].
Nevertheless it was shown in [24] that also PTs in a dark sector are subject to strong constraints, in particular if they happen close to the scale of BBN. The reason is that BBN is a sensitive probe of the Hubble scale at temperatures below the MeV scale, which in turn depends on the total energy density in the Universe, since gravity is universal. Either the energy density in the hidden sector should be transferred to the SM before the onset of BBN at T ∼ 1 MeV, which essentially prohibits PTs below that scale, or the energy should be converted into dark radiation, in which case the dark sector temperature is constrained by N eff .
Viable models should therefore have few degrees of freedom, and still feature a very strong first order PT. The simplest scenario is probably a single scalar field with a nonrenormalizable potential, such as a very light radion or dilaton. Indeed for these models it is known that a strongly supercooled first order PT can occur and produce a large GW background [32][33][34][35][36]. For renormalizable scenarios, the most minimal models that were found in [24] consist of either two real singlet scalars or a U (1) gauge boson with a complex scalar charged under the gauge symmetry. While the majority of the parameter space of these models features a weaker PT, there are benchmark points with α > 0.5 and β/H 100, while still being consistent with constraints from BBN or N eff .
Finally also here it should be noted that PTs with T * ∼ 1 MeV which produce a GW signal stronger than the observed one are now excluded by the NANOGrav data. We are therefore finding the first non-trivial constraints on the dynamics of potential dark sectors around these scales. Of course, to obtain robust limits on concrete models, a reduction of the large theoretical uncertainties in the prediction of the GW signals would be desirable. For some recent progress in this heroic task, see e.g [37][38][39][40].

V Discussion and Outlook
The first hint of a GWB observed by NANOGrav is very intriguing. While the data can be well explained with a single power law, consistent with the expected background from supermassive black hole binaries (SMBHBs), we show here that also broken power law spectra, which are predicted in various extensions of the SM, can well describe the signal.
In both new physics scenarios we considered, the peak of the GW signal is strongly correlated with the relevant mass scale of the new physics, either the axion mass or the mass scale of the new sector that undergoes a phase transition. The PTA data therefore already allows us to narrowly constrain the potential mass range.
Since the data suggests very light new physics, it is already clear that these new particles have to be part of a dark sector that is only very weakly coupled to the SM, otherwise laboratory experiments would have uncovered them already. Yet astrophysical data on BBN and N eff constrain the parameter space of such dark sectors.
For the audible axion scenario, we find parameter regions consistent with N eff for masses around 10 −13 eV and a decay constant of 5 × 10 17 GeV. This region may be probed in the future by the CASPEr-wind experiment [41], and also by future black hole binary merger data through the superradiance effect [42].
A first order PT can explain the data if the transition is very strong and happens at temperatures between 1-10 MeV, or slightly below, if BBN and N eff constraints can be evaded. We have briefly illustrated some dark sector models that are known to satisfy all requirements. Here it will of course be interesting to ask whether concrete realisations can also explain the observed dark matter abundance, and whether they leave observable imprints elsewhere.
Already this first hint of a stochastic GW background in the PTA range provides us with a deep insight into possible new physics explanations of the signal. With more precise frequency binned data it will be possible to distinguish between different models and astrophysical backgrounds such as the one from SMBHBs. It would also be interesting to directly fit a broader range of GW templates to the pulsar timing data, possibly including polarised signals such as the one expected from audible axions. Exciting times lie ahead!