Does NANOGrav observe a dark sector phase transition?

Gravitational waves from a first-order cosmological phase transition, at temperatures at the MeV-scale, would arguably be the most exciting explanation of the common red spectrum reported by the NANOGrav collaboration, not the least because this would be direct evidence of physics beyond the standard model. Here we perform a detailed analysis of whether such an interpretation is consistent with constraints on the released energy deriving from big bang nucleosynthesis and the cosmic microwave background. We find that a phase transition in a completely secluded dark sector with sub-horizon sized bubbles is strongly disfavoured with respect to the more conventional astrophysical explanation of the putative gravitational wave signal in terms of supermassive black hole binaries. On the other hand, a phase transition in a dark sector that subsequently decays, before the time of neutrino decoupling, remains an intriguing possibility to explain the data. From the model-building perspective, such an option is easily satisfied for couplings with the visible sector that are small enough to evade current collider and astrophysical constraints. The first indication that could eventually corroborate such an interpretation, once the observed common red spectrum is confirmed as a nHz gravitational wave background, could be the spectral tilt of the signal. In fact, the current data already show a very slight preference for a spectrum that is softer than what is expected from the leading astrophysical explanation.


Introduction
Gravitational wave (GW) observatories provide us with an additional pair of eyes -or rather ears -to explore the Universe [1][2][3].Since the first detection of GWs in 2015 [4] they help us to better understand our cosmic neighbourhood.With more observatories joining in international collaborations and new experiments being consistently proposed, we can eventually hope to not only detect local, astrophysical sources of gravitational radiation, but also GWs of cosmological origin.To provide context, observables related to primordial nucleosynthesis currently constitute the earliest robust test of the cosmological concordance model [5].GWs could in principle be used to directly probe the cosmological history at even earlier times, because they propagate essentially unperturbed after their production.The corresponding cosmological GW signals will generally form a stochastic gravitational wave background (GWB), in many ways similar to the cosmic microwave background (CMB) [6].
In 2020 the North American Nanohertz Observatory for Gravitational Waves (NANO-Grav) collaboration announced strong evidence for a so-called common red process in the nHz frequency range [3], i.e. a stochastic signal that features the same spectrum of temporal correlations in the pulse arrival times among the set of analysed pulsars.To verify that the signal is really the first detection of a stochastic GWB, the characteristic quadrupolar "Hellings-Downs" (HD) inter-pulsar correlation [7] will have to be confirmed.While this correlation is not expected to be seen in the current data set in view of the still relatively low signal-to-noise ratios [8], its confirmation may be imminent with the upcoming data release.Following NANOGrav's announcement, also the Parkes Pulsar Timing Array (PPTA) [9], the European Pulsar Timing Array (EPTA) [10] as well as the International Pulsar Timing Array (IPTA) [11] published comparable hints for a GWB (see however Zic et al. [12]).Follow-up data releases including more correlation data are thus both anticipated and eagerly awaited, and will contribute to a better understanding of the origin of the observed common red signal [13].
Assuming a confirmation of the inter-pulsar correlations as a GWB, it will be fascinating to further study its properties and to determine the underlying physical processes that generated the GWs.Expected astrophysical sources in this frequency range are mergers of supermassive black hole binaries (SMBHB) [3].In order to explain the observed signal amplitude, however, the local SMBHB density would need to be higher by a factor of a few compared to previous estimates [14 -16], and the question of whether realistic astrophysical models could give rise to a sufficiently strong GWB signal remain the subject of an ongoing debate [17][18][19][20].An alternative and arguably even more exciting possibility is that the signal could be of truly cosmological origin, i.e. the GWB could have been generated by physical processes at redshifts much larger than the z ≲ 2.5 associated with the SMBHB hypothesis.Cosmological sources that have been investigated in the context of the tentative NANOGrav signal include GW emission due to inflation [21], first-order phase transitions [22][23][24], cosmic strings [25][26][27] or the formation of primordial black holes [28,29].
In this work we examine how and under which circumstances first-order cosmological phase transitions can provide a viable explanation of the signal.Such a phase transition is not present in the standard model (SM) of particle physics [30], but there are many well-motivated options for how it could be induced by symmetry breaking or confinement connected to physics beyond the SM [31].In order to match the putative GW signal at nHz frequencies, the preferred temperature for such a phase transition must be at the MeVscale [24].As new physics at this energy scale is very strongly constrained by a large variety of direct experimental searches [32], this immediately implies that the associated new states should only couple very weakly to the SM.In other words, such a phase transition would have to take place in a more or less secluded dark sector (DS) -which could, in fact, also directly be related to the dark matter puzzle [33][34][35].Importantly, even if the DS is only very weakly coupled to the SM, a DS phase transition (DSPT) could impact the successful predictions of Big Bang Nucleosynthesis (BBN) and the CMB due to the extra energy density that is present in the DS (conventionally parametrised as an effective number of new relativistic neutrino degrees of freedom, ∆N eff ) [36,37] or through the late decay of the additional states [38][39][40][41][42].In this work we will consider the two main generic possibilities, namely where the additional energy density 1. fully remains within the DS ("stable DS"), or 2. is subsequently injected into the SM sector ("decaying DS").
It is worth stressing that these options simply refer to different regimes of the inter-sector coupling(s), and hence are a priori equally viable from a phenomenological point of view.In both cases, in particular, we do not assume those inter-sector couplings to be sufficiently large for the DS to thermalise with the SM heat bath.As a result, for strong enough couplings within the DS, the visible and dark sector will generally have different temperatures [43].
Figure 1 illustrates in a nutshell the need to consistently combine cosmological and pulsar timing information when interpreting the NANOGrav results in terms of a DSPT.The blue contours show the results of a naive fit of the DSPT parameters -to be introduced in of the NANOGrav 12.5 year data to a GW spectrum emitted in a DS phase transition, ignoring cosmological constraints.To illustrate the tension with BBN and CMB, the 95% C.L. excluded regions corresponding to ∆N eff > 0.22 are shaded in red [37], cf. the discussion in section 4. The regions in which no percolation (β/H < 3) and an overestimation of the GWB amplitude (β/H < 10) are expected are shaded in grey, see section 2 for further details.
section 2 -to NANOGrav data, without taking into account physically motivated priors on the rate β/H of the phase transition or cosmological constraints.We discuss the former in more detail in section 2, and the latter in section 4. Here, we simply wish to demonstrate that these considerations (as indicated by grey and red shadings, respectively) will necessarily have a major impact on the naively inferred parameter space.One of our main results from a full statistical treatment, including information from cosmology, is indeed that an astrophysical explanation of the common red noise signal is much more credible than a GWB due to a phase transition from a stable DS.When allowing the DS to decay at pre-BBN temperatures, on the other hand, we find that the viable parameter space of DSPTs opens up; in this case, the NANOGrav data can be explained without violating BBN constraints, fitting the pulsar timing data as good as SMBHBs.For earlier work on cosmological constraints on phase transition interpretations of the NANOGrav results, see Refs.[22,44,45].This article is organised as follows.We start by discussing, in section 2, how to predict the GW spectra expected from a DS phase transition.In section 3, we continue with a detailed description of our statistical procedure to analyse PTA data, remarking also on pitfalls and limitations of simpler or more heuristic methods sometimes adopted in the literature.We describe the cosmological constraints on DS dynamics in section 4, and explain how to construct global likelihoods simultaneously taking into account both pulsar timing and cosmological information.We present our results in section 5, before concluding in section 6.In three appendices, we provide further technical details about our analysis.

Gravitational wave backgrounds from dark sector phase transitions
The detailed prediction of the GW signal resulting from a first-order phase transition is a highly non-trivial task.Even though first scaling relations for GWB spectra have been found decades ago [46], the form of the GW power spectrum is still subject to active investigations.Today, there exist a handful of broadly consistent semi-analytical approximations for predicting GWB spectra emitted by first-order phase transitions that are neither too strong nor percolate too slowly [31,47].In deriving these spectra, a source of inspiration has often been a detailed investigation of SM extensions that turn the electroweak symmetry breaking into a first-order phase transition.
Here we are interested in a DSPT whose characteristic energy scale is unrelated to the electroweak scale.We therefore need to adapt the usual approach for calculating GWB spectra from visible sector phase transitions [48][49][50][51][52][53][54][55].We define a DS as an extension to the SM whose particles are coupled so feebly to the SM species that the two sectors do not equilibrate in the early universe.If the couplings within the DS are sufficiently strong, the DS particles instead form a separate thermal bath with temperature T DS that is different from the temperature T SM of the SM photon bath [43].In general, the DS temperature ratio ξ DS ≡ T DS /T SM is a time-dependent quantity.We use ξ perc DS to denote the ratio of the DS temperature in the old phase to that of the SM bath at the time of percolation.Assuming that the energy injection into the DS bath happens instantaneously after percolation, this means that ξ perc DS corresponds to the ratio just before the phase transition.To simplify our analysis, we also assume that the DS reheats instantaneously and that the DS energy density after the transition is dominated by at least one relativistic particle species, such that the speed of sound is given by c DS s = 1/ √ 3 throughout the transition. 1he amount of vacuum energy released in the transition is encoded in the difference ∆θ of the trace of the energy momentum tensor between the old and the new phase [57,58].We can thus introduce two transition strength parameters (2.1) Here, ρ perc SM and ρ perc DS are the energy densities of the SM and DS bath, respectively, at the time of percolation (before reheating).They are related through ξ perc DS as given above, with g perc SM,ρ (g perc DS,ρ ) being the number of relativistic SM (DS) energy degrees of freedom.While α tot determines the amplitude of the GWB spectrum [59], α DS is the relevant quantity to use when calculating the efficiency κ(α DS , v w ) of conversion from vacuum energy to the kinetic energy of the bulk plasma motion.To calculate the efficiency factors we use the semi-analytical approximations given in the appendix of ref. [60] for luminal bubble wall velocities v w .We note that α DS ≫ α tot for small temperature ratios ξ perc DS ≪ 1.Hence, the conversion efficiency of vacuum energy to bulk fluid motion is generally large, κ → 1, such that ultra-relativistic bubble wall velocities v w → 1 are obtained.
The gravitational waves emitted during and after a first-order phase transition are sourced by bubble collisions, sound waves, and turbulences in the perturbed plasma.We consider the cases of dominant bubble wall collision-and sound wave-induced GWB spectra separately.This is a good assumption since bubble wall contributions are dominant for the case of runaway bubbles, whereas sound wave contributions dominate in case of bubble walls with terminal velocities.Such a terminal, yet ultra-relativistic, velocity exists for example when the symmetry being broken is gauged.In that case, the emission of ultra-soft gauge bosons at the bubble walls leads to friction from transition radiation [61], i.e. pressure on the bubble wall that grows with the wall's Lorentz factor γ w , unlike the friction from particle kinematics that saturates for large velocities [62].We do not take into account contributions to the GW spectrum from turbulent motion of the DS plasma after the phase transition [63][64][65][66][67] as this source still requires a better theoretical understanding.
For relativistic bubble walls the differential energy density in the GWB today, with respect to the GW frequency f , can be parametrized as [47,48,59] where β is the inverse timescale of the transition and H is the Hubble expansion rate.The first line of the above equation refers to the case of a GWB dominated by sound wave production and the second line to the case of dominant production through bubble wall collisions, for which we adopt spectra normalized to the respective peak frequencies f p,i as obtained by Refs.[48,51], At the time of production, the peak frequency is given by f em p,bw /β = 0.23 for ultra-relativistic bubble-wall collisions [48], whereas for sound wave-induced spectra it is f em p,sw /β = 0.53 [51].Since the emitted GWs redshift, one has to instead use in eq.(2.2), for the frequencies observed today.Here, we introduced the total entropy degrees of freedom g tot,s , in analogy to the energy degrees of freedom g tot,ρ , as well as a dilution factor D ≡ S 0 SM /S perc tot ≃ 1 with S 0 SM being the SM comoving entropy density today and S perc tot being the total comoving entropy density just before the transition, respectively (see below for a discussion) [59,68].Further quantities entering in eq.(2.2), to be introduced next, essentially fix the overall normalization of the signal.
Starting with the overall prefactor, we use [ which is essentially the same as F gw,0 in ref. [47], but taking into account the changed redshift history that is induced by the presence of a DS.In the second step we adopted the standard ΛCDM values for today's quantities g 0 SM,s = 3.93, g 0 γ,ρ = 2, h 2 Ω γ = 2.473 • 10 −5 , and T 0 SM ≃ 2.35 • 10 −13 GeV [36,69], noting that the total degrees of freedom at early times can potentially receive large DS contributions [70]: g perc tot,s = g perc SM,s + g perc DS,s ξ perc DS 3 . (2.7) Next, the factors Ωsw = 0.012 [51] and Ωbw = 0.11 v 3 w /(0.42 + v 2 w ) = 0.077 [48] are spectrum normalizations obtained in the already previously mentioned simulations; the additional factor of 0.687 ensures that the sound wave spectrum is normalized, s sw (x) d log x = 1/0.687(the normalization of the bubble wall spectrum is already absorbed in the prefactor Ωbw ).The factor of (8π) 1/3 in eq.(2.2) comes from the estimated mean bubble separation [47], and we explicitly keep the relative factor 3 in the sound wave spectrum pointed out in ref. [51] (but erroneously missing in ref. [47]).For the case of sound waves, there is an additional suppression of the spectrum by a factor [47,71] if the sound-wave sources last for less than a Hubble time.Finally, we set κ ϕ = 1 for bubble wall spectra, while we compute κ sw (α DS ) using the high-velocity approximation from ref. [60]; for large α DS , this also results in κ sw ≃ 1 for the case of non-runaway bubbles.
The duration of the phase transition is parametrised by β/H, and the two spectra in eq.(2.2) differ in their scaling with this quantity.In particular, the peak amplitude of the bubble wall spectrum is parametrically suppressed by one power of β/H, such that sound wave spectra can typically fit a given high-amplitude GWB more easily for β/H > 1.As the size of this parameter turns out to be very important for the interpretation of our results, let us briefly discuss what is physically expected.The bubble nucleation rate Γ is suppressed by the tunnelling action S 3 .Close to the nucleation time t * of the phase transition, it is approximately Obtaining smaller values for β/H in a thermally induced phase transition is possible in near-conformal potentials [71,[73][74][75] or at the expense of tuning the scalar potential such that the system is close to meta-stable (when it will not tunnel during the lifetime of the Universe).For much smaller values of β/H ∼ O(1), however, the phase transition might not complete and a phase of eternal inflation might occur.Following the example of ref. [76], we demand β/H > 3 for successful percolation.Moreover, the time between the nucleation of the first bubbles to percolation is about 10/β (see, e.g., ref. [55]).As simulations neglect the expansion of the Universe during the phase transition, which suppresses the GW signal, spectra obtained from such simulations are therefore likely overestimated, or at least subject to sizeable uncertainties, for β/H < 10.
Concerning the exact form of the GW spectra that we adopt here, let us mention that recent results seem to indicate sound shell decays leading to an x −3 rather than x −4 scaling in the UV [55].This scaling has little impact on our results since very low phase transition temperatures are disfavoured in our analysis, implying that the signal is not fitted by the (far) UV tail.Moreover, for wall velocities close to the speed of sound, the sound shell thickness (the thickness of the layer that is put into motion for individual bubbles) becomes imprinted in the fluid motion [54,55,77,78].This leads to an additional knee in the power spectrum and an intermediate scaling x 1 .We also neglect this effect since we focus on wall velocities close to the speed of light.
In all scenarios studied in this work, we set the dilution factor appearing first in eq.(2.4) to D = 1, corresponding to the assumption that there is no significant deviation from the standard cosmological redshift history.This is natural in our analysis of stable dark sectors, 3  as the potential dilution gets sourced by an entropy injection, which can only happen for decaying dark sectors.A dilution factor D > 1 would correspond to a faster expansion than in radiation domination, e.g. if the phase transition were followed by an intermediate phase of early matter domination [59].We checked that this dilution is always negligible also in the case of our decaying dark sector scenario, cf.appendix A.2.
We further set g perc DS,ρ = g perc DS,s = 1.This is not a strong assumption, as the DS degrees of freedom always occur as a prefactor of powers of ξ perc DS in the calculation.Since the experimentally preferred temperature ratios ξ perc DS are small, corresponding to cold dark sector phase transitions, the spectrum's dependence on the precise number of DS degrees of freedom is negligible.As discussed in section 4, this choice also does not influence the cosmological constraints that we impose.Plot showing the NANOGrav "violins" (red) [3], their standard explanation through a power-law spectrum from the inspiral of SMBHBs with A SMBHB = 1.53 × 10 −15 and γ SMBHB = 13/3 (grey) and two phase transition spectra.A characteristic bubble wall collision spectrum is shown in orange.A sound wave-induced GWB spectrum is shown in blue.These spectra correspond to the best-fit parameter points found in section 5 (including cosmological constraints and demanding β/H > 1).The blue arrows indicate how an increase in the phase transition parameters α tot , T perc SM or β/H would shift the spectra.
In summary, we calculate the GWB spectrum based on the following set of DSPT parameters: (2.12) In figure 2 we illustrate the generic influence of increasing α tot , β/H, and T perc SM on the sound wave and bubble wall collision spectra. 4The spectra shown here correspond to the best-fit points of the analysis including cosmological constraints presented in section 5 (with a prior of β/H > 1).

PTA data analysis
In this section, we briefly go over the methods used to analyse PTA data for common red processes.We start by commenting on two often adapted approaches that are used to fit arbitrary GWBs to such signals.After that, we discuss in detail why model comparisons based on global fits require a more rigorous analysis.
In order to fit a given spectrum to the PTA data, one first has to define a deterministic timing model for each pulsar.The timing residuals are then analysed for "white" noise (being uncorrelated in time) and time-correlated "red" noise.To make the fit to different red noise models computationally feasible, the timing model parameters are treated as nuisance parameters and the initial likelihood is marginalized over those.The remaining likelihood for the low-frequency red noise is modelled as a Fourier-series at multiples of the inverse observation period 1/T obs of the PTA.Correlations between observed pulsars then help to further distinguish a GWB from the intrinsic noise, as well as from other sources of red noise (like clock errors or a mismodelled solar system ephemeris).Notably, a GW signal would present itself in the data as a common red process before being confirmed through the characteristic Hellings-Downs cross correlation between the pulsars [7,79].
The NANOGrav collaboration has reported strong evidence for a common red noise signal [3].The Bayes factor between this interpretation and the competing no common red noise (nCRN) hypothesis, in which only pulsar-intrinsic red noise contributes, exceeds 10.000 [3].The arguably simplest explanation of this observation is a GWB sourced by the inspirals of SMBHBs.In this case the spectrum is expected to follow a power-law [80], with f yr = 1/year, H 100 = 100 km/s/Mpc, and the slope γ CP = 13/3.Keeping this slope fixed, the common process spectral amplitude needs to be A CP = A SMBHB ≃ 1.53 × 10 −15 to account for the signal [3].
An easy-to-implement possibility to fit arbitrary GWB spectra to the PTA data is to reinterpret the contour plots in the (A CP , γ CP )-plane produced by the experimental collaborations -e.g.figures 1 and 9 in Refs.[3,11], respectively -to also be valid for a GWB whose spectrum is close to a single power-law in a certain frequency interval.This method has been used in many works [25,26,76,[81][82][83] that aim to explain the common red signal by a GWB from cosmic origins rather than astrophysical sources.Since phase transitions result in GWBs with a broken power-law shape,5 this mapping to a specific combination of A CP and γ CP breaks down around the peak of the spectrum.While this method is often sufficient for estimating the approximate amplitude of the signal, it is not powerful enough for making a proper model comparison between different signal hypotheses.
Another tempting possibility to fit an arbitrary GWB spectrum to the common red noise signal is to use the results of a free spectral analysis to the PTA data [3].In that analysis the assumption of a power-law GWB was dropped in favour of free spectral amplitudes at the aforementioned frequencies of integer multiples of 1/T obs .The posteriors of these spectral amplitudes are depicted as the infamous "violins" reproduced in figure 2. Assuming that the total likelihood to fit an arbitrary spectrum to the common red signal factorizes into five dominant likelihoods, at the Fourier frequencies 1/T obs , ..., 5/T obs , then allows to directly fit spectra that can in principle deviate arbitrarily from a power-law (see, e.g., Refs.[23,[85][86][87]).This approach, however, systematically underestimates the uncertainties of the full calculation as it by definition neglects correlations between amplitudes at different Fourier frequencies [87].Another significant limitation is that the posteriors for the respective frequency bins are calculated with a finite prior range of signal amplitudes.Adding to the fact that the tails of these posteriors are not necessarily sampled well enough, this implies that the violins cannot be used in any statistically meaningful way for signal amplitudes too far from their most likely values.For example, the finite extent of the violins shown in figure 2 would strictly speaking imply that the 'no-signal' (nCRN) hypothesis is excluded with arbitrary confidence -while instead it is only disfavoured by a Bayes factor of ∼ 10 4 − 10 5 .
Crucially for us it turns out that cosmological constraints can force a potential DSPT to be so weak that the resulting GWB spectrum only contributes negligibly to the measured common red noise.In that case, the signal is fit by fine-tuning the individual pulsar-intrinsic red noise amplitudes.Since for both the mapping to a single power-law and in the free spectral analysis the pulsar-intrinsic red noise components are already marginalized over, it is therefore no longer possible to use one of these quick methods.Instead, in order to treat correlations consistently, a full evaluation of the likelihood is required.
The evaluation of the full PTA likelihood to fit cosmological GWB spectra was so far used only in a rather short list of works [24,88], due to the large numerical cost of the likelihood evaluations, and -as far as we are aware of -never in a context that included further constraining likelihoods, like in our case from cosmology.
To interpret the NANOGrav data in terms of a DSPT we first construct a likelihood L PTA (θ PSR , θ DS ) for fitting the timing residuals to a given set of pulsar-intrinsic red noise parameters θ PSR and a common red spectrum that depends on the DSPT parameters θ DS .Each pulsar's red noise is fitted by a power-law with an amplitude A red,i and slope γ red,i , in analogy to eq. (3.1).To include constraints from cosmology on the available dark sector parameter space θ DS , we further construct a likelihood L cosmo (θ DS ) in section 4. We multiply this likelihood with the PTA likelihood to obtain a global likelihood, In this work we concentrate on the NANOGrav 12.5 yr data [3], for which the full set of arrival time data, the pulsar white noise parameters as well as a tutorial on how to use these resources publicly available [89]. 6In this data set [3], a total of 47 pulsars were taken into account, out of which we consider those that were observed for at least three years, as done in the original analysis.From the remaining 45 pulsars, we treat the pulsar J1713+0747 as advertised in ref. [89] due to the probably mis-modelled noise process found in the dropout analysis [3].The parameter space we evaluate therefore consists of 90 nuisance parameters θ PSR = {A red,i , γ red,i } for the pulsar-intrinsic red noise, adding to our four (five) DSPT model parameters for the case of a (decaying) DSPT.To evaluate this high-dimensional global likelihood in a numerically feasible way, we implement the DSPT spectra and L cosmo (θ DS ) into the codes enterprise and enterprise extensions [90,91].To sample over the global likelihood we use PTMCMC [92].The Markov chain Monte Carlo (MCMC) chains are evaluated using numpy and scipy [93,94], and triangle plots are generated using matplotlib and a customized version of ChainConsumer [95].For performing model comparisons, finally, we calculate Bayes factors using the product space method [96][97][98][99], which we briefly review in appendix C.1.

Cosmological constraints
The past decades of observational cosmology have provided a large amount of data which allow us to confidently reconstruct the cosmological evolution up to MeV-scale temperatures.Most notably these include observations of the CMB, both in terms of the spectral shape [100] and anisotropies [36], and the primordial light element abundances produced during BBN [101].These observations are in very good agreement with the standard ΛCDM model and with each other [101], implying that any changes to ΛCDM at temperatures below a few MeV can have observational consequences and need to be checked for consistency with available CMB and BBN data [36][37][38][39][40][41][42].
For a phase transition to produce a strong GW signal a sizeable amount of vacuum energy needs to be released in the transition, most of which is subsequently converted into DS energy density as only a small fraction ends up in the GWB.This additional energy density could change the well-tested cosmic expansion history, possibly even long after the transition itself.To understand whether NANOGrav may observe the remnants of such a phase transition we therefore need to include a cosmological likelihood L cosmo when analysing the PTA data.Specifically, we include information from the primordial light element abundances and CMB anisotropies into our analysis as described below. 7

Stable dark sector
If the entire DS energy density after the phase transition is contained in light degrees of freedom, this contributes an additional radiation energy density that can be described by a (potentially large) additional contribution ∆N eff to the effective number of neutrinos N eff = N SM eff + ∆N eff , where N SM eff = 3.044 [107] is the SM prediction for ΛCDM cosmology. 8The effective number of neutrinos affects the predictions of BBN as well as the CMB power spectrum and is constrained by observations to N eff = 2.941 ± 0.143 [37].This bound can be modelled by a Gaussian likelihood L BBN+CMB (N eff ), i.e.
enabling us to implement the cosmological constraints for the case of a stable DS.The number of degrees of freedom in the DS does not have any relevant impact on cosmological constraints as the available latent heat would simply be distributed among the different species, leaving the total injected energy density unchanged. 9We therefore set g perc DS,ρ = 1 in our analysis.For more details we refer to appendix A.1.

Decaying dark sector
If there are additional small inter-sector couplings -which happens very naturally due to possible 'portal' couplings such as Higgs mixing for scalars or kinetic mixing for dark photons -the energy density injected into the DS will subsequently be transferred to the SM heat bath via decays of DS particles.In this case cosmological constraints in general depend on the lifetime, mass, and coupling structure of the decaying particles.As a simple concrete example and for minimality, we consider the DS after the phase transition to consist only of one bosonic degree of freedom ϕ decaying into photons or electrons with a lifetime τ ϕ that we sample over.This is a natural setup as a light scalar degree of freedom with mass below the phase transition temperature is generally expected to exist and e.g.mixing with the Higgs is generally allowed.Given that the NANOGrav data prefer an MeV-scale phase transition temperature, also the mass of ϕ is expected to be of this order, m ϕ ≃ MeV.In the example of Higgs mixing very short lifetimes τ ϕ correspond to a sizeable Higgs mixing angle θ which is constrained by a number of laboratory experiments to be θ ≲ 10 −4 for small masses [109,110].Translating this constraint for MeV-scale masses into the lifetime results in τ ϕ ≳ 10 −2 s whereas cosmological constraints roughly require τ ϕ ≲ 10 −1 s, so that some allowed region remains in this scenario. 10As already indicated, the resulting cosmological constraints largely depend on the lifetime τ ϕ of the particle with only a very mild dependence on the mass as long as the energy density in ϕ is not very suppressed. 11o simplify our analysis we therefore fix the mass of the decaying particle to m ϕ = 5 MeV, noting that the results will be very similar for other choices of the mass in the MeV-range.To implement the cosmological constraints, we use a likelihood incorporating BBN and the CMB power spectrum constructed from the results of [40,44] as detailed in appendix A.2.

Results
We now discuss the results of our work, based on global fits of the model setups described in the previous sections.We start by considering phase transitions in a stable dark sector (section 5.1), and then turn to a decaying dark sector that thermalises with the visible sector some time after the phase transition (section 5.2).We finally compare different GW explanations of the signal in terms of their respective Bayes factors, in section 5.3, and discuss how future data might strengthen the DSPT interpretation.

Stable dark sector phase transitions
Let us first focus on GWs that are mainly emitted as a consequence of the bulk motion of the DS plasma after the phase transition, i.e. a GWB dominantly produced through sound waves.The full set of model parameters to describe such a scenario for a stable DSPT, as introduced in section 2, is given by {α tot , β/H, T perc SM , ξ perc DS } . (5.1) We sample over these input parameters with flat log priors, as well as over 90 nuisance parameters θ PSR for the pulsar-intrinsic red noise, based on the combined PTA and cosmological likelihood given in eq.(3.2).For a full overview over parameters and prior ranges, see also table 1 in appendix C.6.We show the resulting corner plot of posterior distributions for the four model parameters in figure 3, to which we add the derived parameter ∆N eff .Allowing for inverse time scales down to β/H > 1 (light blue) formally results in a good global fit of the pulsar timing residuals, as indicated by the compact ellipsoid-like posterior regions where the NANOGrav signal is explained by the GWB.Such small values of the transition rate would however suppress the GW spectrum w.r.t. the commonly adopted parametrization, cf. the discussion in section 2, and, for β/H ≲ 3, likely not even lead to successful percolation.We therefore also show, in the same figure, the results of a fit with a more restrictive prior of β/H > 10 (dark blue).In this case, the best-fit region moves to a somewhat larger value of α tot , but it also becomes apparent that there is no longer a single preferred region in the model parameter space.Instead, the combined data now shows a similar preference for a very weak DSPT-induced GW signal with correspondingly weak cosmological constraints, where the NANOGrav signal is not explained by the GWB but absorbed in the pulsar-intrinsic noise parameters.This indicates that also the best-fit region no longer corresponds to an equally plausible interpretation of the combined data set.
The reason is that cosmology adds a constraint on ∆N eff which effectively translates into a constraint on the phase transition strength α tot .In terms of fitting the NANOGrav signal, the required lower value of α tot can partially be compensated by a lower value of β/H (see also figure 2).And indeed, comparing the posterior distributions for β/H in figure 3, we can see that β/H always sticks closely to the lower prior boundary -which is qualitatively different from the analysis without cosmological constraints, cf.figure 1, where inverse timescales of β/H = O(10 − 100) are favoured.Increasing the lower prior bound on β/H therefore directly induces a shift towards larger α tot .At the same time, a larger inverse timescale β/H also means smaller bubbles at the time of their collision, and hence a larger peak frequency in the spectrum (see again figure 2).This is compensated for by a lower percolation temperature T perc SM , which by itself would lead to smaller peak frequencies.Finally, we can identify in figure 3 an upper bound on the initial temperature ratio, which again is a direct consequence of the constraint on ∆N eff .For ξ perc DS ≳ 0.7, in particular, eq.(A.7) would imply a violation of this constraint already for a single dark sector species that is not non-relativistic before the transition [72].
Overall, these effects push the posterior for ∆N eff towards higher values.Since this is strongly punished by the cosmological part of the likelihood, however, that also explains the already mentioned appearance of a second preferred parameter region, characterized by a weak DSPT that corresponds to ∆N eff ≃ 0 (at the price of an unobservably small GW signal).We confirmed that this region of parameter space is indeed explained by fine-tuning the pulsar-intrinsic red noise parameters, rather than by a GWB, by directly comparing the marginalized posteriors of the nuisance parameters θ PSR between the two chains depicted in light blue and dark blue.Let us stress that this parameter range would have been impossible to reliably infer with simpler statistical methods, i.e. by re-fitting a power-law common red spectrum described by (A CP , γ CP ) or by using the free-spectrum "violins" (see the discussion in section 3).
From the above discussion, one would expect the Bayes factor between the DSPT and the no nCRN hypothesis to decrease when increasing the lower prior bound on β/H, as higher and higher values of ∆N eff are needed to explain the combined data.We confirm this expectation explicitly in figure 4. Here, each of the coloured dots corresponds to a separate MCMC chain that was employed to determine the Bayes factor between the two models by using the product space method explained in appendix C.1.The corresponding lines serve as a cross-check for the prior dependence of the Bayes factors, see appendix C.4 for further details.Yellow dots and lines refer to the case of a stable DSPT where the GWB production is dominated by sound waves; this corresponds to the same model as shown in figure 3.
For comparison, we also show the case of a GWB that is mostly due to bubble wall collisions (red).Just from the point of view of the resulting spectrum, cf.figure 2, one might expect that this could be a viable alternative.Compared to sound waves, however, bubble wall spectra receive an additional parametric suppression of h 2 Ω peak GW by a factor of (β/H) −1 .This induces the need for a larger α tot and hence an even stronger constraint on ∆N eff , making the GWB hypothesis worse than the nCRN assumption for β/H ≳ 5. To further illustrate these considerations, we refer to figure 6 in appendix B, showing the posterior distribution of the bubble wall spectra for different prior choices on β/H.Note also that neither figure 3 nor figure 4 include the expected suppression in the GWB spectra for β/H ≲ 10, which would further decrease our Bayes factor estimates.
Overall we therefore come to the conclusion that a stable DSPT can hardly compete with the alternative SMBHB explanation of the NANOGrav timing residuals, once one takes into account cosmological constraints from BBN and CMB.For β/H > 10, in particular, the Bayes factor between a DSPT explanation and the nCRN hypothesis is only O(10) even in  ).We also indicate how to translate the Bayes factor to Jeffrey's scale (coloured shadings) as well as Z-scores from frequentist statistics (right y-axis, cf.appendix C.3).In the grey shaded area (β/H < 3), a constant nucleation rate is not sufficient to drive percolation.
the favourable case of sound wave-induced GWBs -much smaller than the factor of ∼ 10 4.5 that is claimed for a GWB from SMBHBs [3].

Decaying dark sector phase transitions
We next consider a DS that couples sufficiently strongly to ordinary matter such that it can decay after the phase transition.A decay long before BBN, in particular, is not subject to the cosmological constraints that played such a decisive role for the case of a stable DSPT.On the other hand, a phase transition happening too early would produce a GWB at too high frequencies to be compatible with the NANOGrav signal.It therefore is a non-trivial question whether any relevant parameter space remains where PTA and cosmological data are indeed compatible.For concreteness, we assume the decay of a dark Higgs boson as detailed in section 4 that decays with a lifetime τ ϕ , such that our model parameters read For the dark Higgs lifetime we adopt a log prior ranging from 10 −6 s to 10 2 s; the remaining parameters we treat as in the previous section (with β/H > 1).We show the resulting triangle plot for this model in figure 5.In a nutshell, we find that the common red spectrum observed by NANOGrav can be explained as long as τ ϕ ≲ 0.1 s and T perc SM ≳ 2 MeV.Larger lifetimes, corresponding to temperatures smaller than 2 MeV, are strongly constrained as the decays occur after the onset of BBN or neutrino decoupling. of the neutrino and photon temperatures, leading to a negative contribution to ∆N eff .If the percolation temperature T perc SM drops to values close to the temperature of neutrino decoupling strong constraints arise independent of the lifetime τ ϕ .Note that we implemented the results from ref. [44] for simplicity as a sharp cut-off enforcing T perc SM > 2 MeV in our likelihood, cf.appendix A.2. Adopting a more accurate likelihood would result in a smoother transition of the posterior in the range T perc SM ≈ 1 − 2 MeV, the main effect being a slight increase of the maximal possible value of β/H.

Such decays change the time-temperature relation of the SM heat bath and alter the ratio
Compared to our analysis of stable DSPTs, the posterior for the inverse timescale β/H is relatively flat up to β/H ∼ 30, implying a very limited prior dependence.The underlying reason for this is the possibility of dumping the liberated energy density into the SM photon bath before the neutrino decoupling at around 2 MeV, thereby evading cosmological constraints and hence opening up for large phase transition strengths α tot ≳ 0.1 to fit the common red signal even for β/H ≳ 10.This however only works up to the point when β/H becomes so large that its effect on the peak frequency can no longer be compensated for by a correspondingly lower percolation temperature, cf.figure 2.
In figure 4 we also indicate the Bayes factor for the decaying DSPT scenario (blue).
As anticipated, the prior dependence on β/H is much less severe than in the scenarios discussed previously.In particular, this shows that a GWB from a decaying DSPT is a viable explanation of the common red noise signal even for β/H > 10.Quantitatively, the model evidence is a factor of ∼ 200 larger than that of the nCRN hypothesis, corresponding to 2.5 σ or a "decisive" evidence on Jeffrey's scale.

Comparison with other GW sources
Let us next address in more detail the question of how a DSPT interpretation of the signal compares to alternative GWB hypotheses, in particular the leading astrophysical explanation of an SMBHB-induced GWB, and how future data will help to further distinguish these two.We start by pointing out that the DSPT spectra actually fit the common red spectrum in the NANOGrav data slightly better than the SMBHB spectra.Naively, this is already expected from figure 2, and we demonstrate this in more detail in appendix C.5.Nevertheless the maximal Bayes factor (with respect to nCRN) that we find is only about 10 2.5 , significantly smaller than the ∼ 10 4.5 typically quoted for SMBHBs.We checked explicitly that the reason is not connected to the goodness of fit, but entirely due to prior volume effects: Bayesian statistics dutifully renders the simple SMBHB explanation of the data more credible than the apparently more complicated DSPT model.It is however important to keep in mind that the amplitude A SMBHB of the astrophysical signal is not necessarily a single fundamental parameter, as assumed in deriving a Bayes factor of 10 4.5 , but effectively derived from more fundamental parameters in a more complicated astrophysical model, such as astrophysical merger timescales or the mass-, redshift-, and spatial distributions of the SMBHBs.A full Bayesian analysis should thus also consider constraints on these fundamental parameters, as e.g.done in ref. [14] without fitting the NANOGrav data.This would decrease the formal evidence for the SMBHBs interpretation because (i) the prior volume increases due to the intrinsic parameters that the amplitude depends on, and (ii) astrophysical models predict amplitudes that are smaller than those inferred from the NANOGrav data [11,[14][15][16].Analogously, of course, the DSPT parameters α tot , β/H, T perc SM , ξ perc DS , and τ ϕ are in reality derived quantities from a given SM extension whose independent parameters are masses and couplings.Evaluating specific models where these underlying parameters are known would be interesting and clearly deserves further study, but is not the aim of our more model-independent analysis.
Concerning near-future perspectives, the first thing to be confirmed is that NANOGrav indeed sees a GWB by verifying that the correlations follow a Hellings-Downs curve [7].This is is expected to happen relatively soon, within the next weeks to years [13], and thus possibly already for the pending release of the 15 yr data set.Even with a confirmed GW signal one may initially expect the relative odds between SMBHBs and an alternative DSPT explanation to remain about the same -as long as the preferred spectral shape of the GWB remains similarly vaguely determined as of today.This is simply because, as discussed above, the prior volume is currently the main driver to distinguish between these two hypotheses.In other words, the curves in figure 4 will likely all shift towards higher Bayes factors in favour of DSPTs, because the relative pull by ∆N eff towards lower values of α tot will reduce, but the Bayes factors for the SMBHB interpretation will improve by roughly the same factor.
A first lever-arm to really distinguish the DSPT from the SMBHB scenario will be a determination of the spectral index α = (3 − γ)/2.This quantity could be determined by up to 40 % relative accuracy once Hellings-Downs is confirmed [13].Measuring a value of γ that is at odds with the value of γ SMBHB = 13/3 expected for SMBHBs [80] would then be a strong indicator for a DSPT instead of an astrophysical interpretation of the GWB.Current data, in fact, already slightly prefer γ > 13/3, cf.figure 2. This could easily be accommodated by the GWB expected from a DSPT, both close to the peak amplitude and in the high-frequency tail of the spectrum.An even clearer signature would be a broken powerlaw, but to test for such a feature will require even higher signal-to-noise ratios and hence take correspondingly longer time to probe observationally.A valuable source of information to further disambiguate between SMBHBs and GWBs of primordial origin may lie in the anisotropies of the GWB [111][112][113], in particular via the polarization of the GWB [114][115][116][117][118][119][120] and possibly the detection of singular GW sources as already searched for [121,122].
Let us finally stress that not only the evidence for a nHz GWB will change in the future; also the competing cosmological constraints on phase transitions are subject to improvements.Upcoming experiments like the Simons observatory and CMB-S4 measurements in combination with surveys of large-scale structure will be able to push the limits on ∆N eff by about an order of magnitude [123][124][125], contributing to the tension on the stable dark sector explanation we investigated above (and reducing the parameter space for a decaying DSPT).Measurements of CMB spectral distortions with PIXIE [126] will give additional and complementary constraints [102], which would be relevant even in the case of a decaying dark sector that avoids constraints on ∆N eff .In fact, a confirmation of spectral distortions in the CMB could in principle even provide supporting evidence for such a scenario.It will therefore remain crucial to include all relevant cosmological observables as well as new PTA data to eventually decide on the most probable origin of a putative GWB signal.

Conclusions
We investigated the appealing possibility that the common red spectrum observed by the NANOGrav collaboration in their 12.5 yr data set [3] is due to a dark sector phase transition (DSPT) just before the onset of Big Bang Nucleosynthesis (BBN).For the first time, we performed a global analysis on pulsar timing array (PTA) data from a (potential) gravitational wave background (GWB) including constraints from BBN and the cosmic microwave background (CMB) anisotropies.
We found that a dark sector (DS) undergoing a phase transition can in principle explain the measured signal with a goodness of fit that is comparable to, or even better than, that of the standard astrophysical explanation in terms of a stochastic GWB from supermassive black hole binaries (SMBHBs).However, if one accounts for 1. the changes in the early element abundances that the energy density released during the phase transition would induce, 2. the impact on the CMB anisotropies trough a contribution to ∆N eff , and 3. possible issues for transitions with β/H < 10, connected to percolation and an overestimation of the produced GWB, the possibility of a stable DSPT no longer gives a good fit to all available data. Figure 1 provides an intuitive illustration of this tension, by directly confronting the above constraints with the results of a naive DSPT analysis of the NANOGrav data that ignores them.Fully including all relevant constraints in a global fit, the available parameter space is indeed significantly reduced, cf.figure 3.
On the other hand, there is no intrinsic reason why a DS should be stable on cosmological timescales.In particular, tiny interactions with the visible sector (e.g. through small portal couplings [127]) could well lead to a decay before neutrino decoupling at T perc SM ≳ 2 MeV.We find that such a decaying DSPT scenario remains a compelling alternative to the more conventional SMBHB hypothesis for lifetimes τ ϕ ≲ 0.1 s, cf.figure 5. We arrived at this conclusion by further taking into account constraints on electromagnetic energy injection from decaying dark scalars [40] and on the reheating of the photon bath after a phase transition [44].Compared to the no-GWB hypothesis, we find a Bayes factor that indicates a decisive evidence for the DSPT interpretation even for a prior of β/H > 10 on the transition rate.The currently maximal value of this quantity that is compatible with the data, β/H ≲ 50, still indicates the need for a relatively slow transition; further model-dependent research will be needed to investigate how this can be implemented in a given SM extension.
We also studied the effect of prior choices on the absolute scale of Bayes factors, finding that prior volume effects are highly relevant when comparing SMBHB and DSPT explanations of the NANOGrav data.The SMBHB interpretation, in particular, seemingly only requires one parameter to fit the signal, namely the amplitude A SMBHB ≃ 1.53 × 10 −15 .We however argue that A SMBHB should rather be treated as a derived quantity that depends on several intrinsic, independently measured astrophysical quantities [14].This would reduce the difference between the Bayes factors above 10 4 for the SMBHB explanation [3] and the Bayes factors of O( 102 ) that we find for the decaying DSPT interpretation.
We remain excited about the pending 15 yr data release of the NANOGrav collaboration, which might already confirm the Hellings-Downs curve needed to state the detection of a GWB [7,13].While we do not expect a definite answer concerning the origin of the signal within the coming months, we are confident that additional PTA data as well as complementary constraints from cosmology [102] will help to disambiguate between different models.The earliest evidence we can hope for in favour of a cosmological origin of the signal will be any deviation from the spectral slope of γ SMBHB = 13/3 expected for SMBHBs.With further data it will be ever more crucial to fully include complementary constraints such as from BBN and CMB when assessing different signal models.In this work we have made a first step in this direction, thereby contributing to moving the realm of testable cosmology to pre-BBN times.

A Details on cosmological likelihood
In this appendix we provide details on the computation of the relevant quantities for the cosmological likelihoods as well as the mapping to previously published results.

A.1 Stable dark sector
For the model of a stable DS after the phase transition we assume that the energy density in the DS after the transition is contained purely in radiation.Note that constraints would only be more stringent if the DS energy density would start to redshift like matter, increasing the ratio to the SM energy density.Furthermore, we focus on transitions before the onset of BBN and neutrino decoupling.At the observationally relevant temperatures we can then model the DS energy density by an additional contribution ∆N eff to the effective number of neutrinos N eff .This assumption is validated a posteriori by the results of our MCMC chains.
At the time of percolation t perc the DS energy density in radiation ρ perc DS can be quantified through the effective number of relativistic degrees of freedom g perc DS,ρ and the DS temperature Assuming an instantaneous reheating of the DS after the phase transition, one finds where the index reh (perc) indicates a point in time immediately after (before) the DS reheating.The strength of the phase transition is characterized by α tot ≡ ∆θ/ ρ perc DS + ρ perc SM , where ∆θ ≃ ∆V for sufficiently strong transitions.We thus find In terms of the temperature ratio ξ perc DS before the phase transition we obtain Energy density in radiation additionally to the SM can be quantified by its contribution to the late time effective number of relativistic neutrino species By subtracting the SM prediction in a ΛCDM cosmology N SM eff = 3.044 [107] we can define the contribution of the DS to the effective number of relativistic neutrino species where we inserted the SM degrees of freedom today and assumed the DS degrees of freedom g DS,ρ = g perc DS,ρ to remain constant.This reproduces eq. ( 1) from ref. [44].Observationally, N eff affects the predictions of BBN as well as for the CMB power spectra.Combining data sets from the Planck satellite with observations of the primordial abundances of deuterium and helium-4 and marginalizing over the baryon-to-photon ratio η = n b /n γ , ref. [37] finds N eff = 2.941 ± 0.143.We approximate this by a Gaussian likelihood. 12 cosmo (α tot , ξ perc DS , T perc SM ) = 1 where µ N eff = 2.941, σ N eff = 0.143, N SM eff = 3.044, and the normalization needed for the correct determination of Bayes factors is given by As ∆N eff is simply a function of the phase transition parameters, we can now evaluate the cosmological likelihood for a DS phase transition with given parameters.By multiplying this likelihood to the one for the timing residuals we can simultaneously fit the NANOGrav data and check for consistency with cosmological observables.

A.2 Decaying dark sector
The constraints from BBN and the CMB on a decaying particle with MeV-scale mass have been calculated in ref. [40] for a general setup.The cosmological evolution was started at an SM temperature 13 T cd SM = 10 GeV and corresponding time t cd with a Bose-Einstein distribution function of ϕ of temperature T cd DS = ξ cd DS T cd SM and zero chemical potential, i.e.
The corresponding number density is To use the results and constraints from ref. [40] we therefore need to map to this scenario.Sufficiently strong intra-sector couplings, i.e. such that the interaction rate is much larger than the Hubble rate, lead to a self-thermalisation of the DS quickly after the phase transition.If the processes are further ϕ-number violating, we can expect the chemical potential of ϕ to vanish.The phase-space distribution function after the DS has reheated at t reh can hence be described by a Bose-Einstein distribution function with corresponding number and energy densities To compute the energy density ρ DS ϕ we again assume an instantaneous reheating with ρ perc DS = (π 2 /30) g perc DS,ρ (ξ perc DS T perc SM ) 4 and compute (cf.eq.(A.3)) For fixed mass and lifetime, cosmological constraints are mostly driven by the comoving number density of a particle.We hence compute the comoving number density n reh ϕ , which can be obtained from ρ reh ϕ by numerically solving eq.(A.14) for T reh DS = ξ reh DS T perc SM and then computing n reh ϕ from eq. (A.13).This last step of mapping ρ reh ϕ to n reh ϕ can be accelerated considerably by utilizing tabled data for zero chemical potential and fixed mass m ϕ .By redshifting number densities using comoving conservation of SM entropy density s SM , and then inverting n cd ϕ numerically to obtain the required temperature ratio ξ cd DS , cf. eq.(A.11), we can map the case of a DS phase transition to the case computed in ref. [40].We comment below on the validity of this mapping.
To construct a likelihood, we use the calculated primordial light element abundances, their theoretical errors from nuclear rate uncertainties, and N eff underlying figure 6 (left) of ref. [40].These are compared to the recommended values of the observed primordial light element abundances of deuterium D/ 1 H obs = (2.547± 0.025) × 10 −5 and the mass fraction of helium-4 Y obs p = (2.45 ± 0.03) × 10 −1 [101] as well as N obs eff = 2.99 ± 0.17 [36].As the BBN calculations use the baryon-to-photon ratio η (or equivalently the baryon density in the Universe) from ref. [36] as an input, but N eff is affected by the DS decays injecting energy in the photon heat bath, we use the best-fit value of η for given N eff from figure 26 (Planck TT, TE, EE+lowE+lensing+BAO) of ref. [36].Note that due to the strong dependence of the prediction of D/ 1 H on η we need to propagate the uncertainty in the determination of η, effectively marginalizing over η, such that the total observational error becomes η obs D/ 1 H = 0.035 × 10 −5 [40].After using the N eff -dependent best-fit value there still remains the aforementioned constraint on N eff .The cosmological likelihood is given by a product of Gaussian likelihoods with mean at the observed values and total errors obtained by summing the observational and the theoretical error (negligible for N eff ) in quadrature.
In the calculation above we started from eq. (A.15), crucially assuming that the DS and the SM only thermalise after reheating, i.e. that (inverse) decays of ϕ only become relevant after reheating.This holds as long as τ ϕ > t reh , where t reh can be obtained in radiation domination via the Hubble rate as H(T reh SM ) ≃ 1/(2t reh ).For shorter lifetimes τ ϕ < t reh the DS and the SM thermalise during or shortly after reheating.Therefore, the temperature ratio after reheating is generally (close to) one, ξ reh DS = 1, but T perc SM is no longer equal to T reh SM as reheating goes into the DS as well as the SM due to the decays of ϕ.We find T reh SM = T reh DS by solving for T reh SM with ξ reh DS = 1 in eq.(A.12).With T reh SM , ρ reh ϕ , and ξ reh DS = 1 we can find n reh ϕ as before and map to the results of ref. [40] using eq.(A.11).This still assumes that there is no change in comoving number density of ϕ or comoving SM entropy density in the calculation of ref. [40] between the temperature T cd SM = 10 GeV and the numerical value of T reh SM , i.e. decays and inverse decays can be neglected up to the corresponding numerical value of the time t reh .Strictly speaking, this assumption is not valid any more due to thermalisation around a time t ∼ τ ϕ < t reh .However, successful thermalisation erases all knowledge of initial conditions, implying that we only have an inaccurate mapping if thermalisation itself has observable consequences, i.e. it occurs during BBN.Hence, also the reheating process would need to occur during BBN, and we, in fact, have to include the effect of the reheating process (also reheating the SM) on BBN.These were studied in ref. [44] (under the assumption ξ perc DS = 0), finding that α * ≳ 0.07 is constrained if T perc SM ≲ 2 MeV for reheating into photons.To compare this to our case we note that in ref. [44]  The complete cosmological likelihood is given by where the normalization is N = L Yp × L D/ 1 H × L N eff for standard ΛCDM cosmology with only the SM contributing to the energy density during BBN.The likelihood given above accurately describes the relevant cosmological constraints on a decaying scalar in a quite model-independent way.We however needed to make some assumptions along the way, either to assure the numerical feasibility of our calculations or to keep the number of parameters describing the decaying dark sector scenario low to allow for a phenomenological interpretation of our results.Note that these assumptions generally do not affect our conclusions that the NANOGrav signal can be explained well by a decaying DSPT as long as the energy from the DS is injected into the SM before the onset of BBN and neutrino decoupling (i.e.τ ϕ ≲ 0.1 s and T perc SM ≳ 2 MeV).Still, we discuss the effect of the assumptions below for completeness.
A first simplification is the choice of a particular mass m ϕ = 5 MeV due to readily available data from ref. [40].Generally, the dependence on the mass is expected to be very mild.Only when the mass is small enough, m ϕ ≲ 2 MeV, for the decaying particle to thermalise with the SM heat bath as a relativistic particle for τ ϕ ≲ 0.1 s and act as an additional relativistic degree of freedom, arbitrarily small lifetimes can be constrained [40].
Next, we use results for decays into photons.The results of ref. [40] show that decays into electron-positron pairs give very similar constraints for m ϕ > 2 m e .For m ϕ < 2 m e decays into electron-positron pairs are kinematically forbidden such that a corresponding coupling will, in fact, mostly lead to decays into two photons and the constraints on the total lifetime are just as in the case for decays into photons (whereas constraints on couplings would become suppressed by loops and additional SM couplings).Other interesting portals of the two sectors that lead to a thermalisation between them are imaginable, e.g.due to the kinetic coupling of a dark photon field with the SM photon [128].These couplings are however beyond the scope of our work due to the additional model-dependent constraints.
When considering a particular DS model, the assumption of instantaneous reheating might not be justified, e.g. if thermalisation of the DS is not sufficiently fast due to smaller intra-sector couplings or if there is additional dilution of the relevant energy densities between t perc and t reh in eq.(A.2) due to a longer reheating process.
By only redshifting the number density of the decaying particle via eq.(A.16) we neglected a possible change in comoving number density due to decays and inverse decays in the calculation of ref. [40] between the starting temperature T cd SM = 10 GeV and the numerical value of T reh SM .While our treatment using the results of ref. [44] redeems this shortcoming adequately, a full model-dependent study should calculate the whole cosmological evolution of the DS including all relevant energy transfers between the DS and the SM.In this way, also the dilution factor D entering the GWB spectrum would be calculated [59], which we simply set to 1.This assumption is valid as long as the decaying particle does not become non-relativistic for an extended period of time before its decay.We checked a posteriori how large the dilution factor would be for the regions favoured by NANOGrav data, cf.figure 5, using the dilution sub-package of TransitionListener from ref. [59].Within the 1 σ contour, the dilution factor deviates by at most 2 % from 1.The maximal dilution factor within the 2 σ contour is D ≈ 2. Hence, our choice of setting D = 1 is not only conservative but also well justified in the relevant parameter space.
Finally, our assumption of a negligible chemical potential might be violated in some setups.As long as number-changing processes are efficient, the chemical potential is driven to zero.However, if the decaying scalar is not completely relativistic after reheating and the scalar self-couplings are low, number-changing processes can cease to be relevant and a chemical potential can develop [129,130].Such a scenario would likely require a full numerical solution of the respective Boltzmann equations and is therefore beyond the scope of this work.

B Posterior distribution of GWB spectra
To further demonstrate the possible tension between cosmological constraints and the interpretation of the NANOGrav data in terms of a DS phase transition, we illustrate the posterior Posterior distribution of GWB spectra for different β/H priors distribution of GWB spectra for different parameter scans in figure 6. 14 The orange curves show the distribution of bubble wall collision spectra from a dark sector phase transition with β/H > 1.For bubble wall collision spectra with β/H > 10 (yellow curves), the tension in the data becomes apparent through the low signal amplitudes.The spectra cannot explain the common red spectrum, as depicted by the red "violins", due to the strong constraints from ∆N eff .When going to a decaying DS (violet), a phase transition with β/H > 10 can clearly still be consistent with the measured signal.
This visualization also clearly shows that a mere fit using the violins is not sufficient for our analysis including cosmological constraints due to the unknown likelihood for small GWB spectra, where the signal is absorbed in finely tuned pulsar-intrinsic red noise parameters, cf.section 3.

C.1 The product space method
The calculation of Bayes factors in this work relies on the product space method, which we briefly review for completeness.To compare two models, a hyper-model H is introduced, whose parameter space contains the Cartesian product θ H of the two sub-models' parameters θ 0 and θ 1 as well as an additional model index n that can run from −0.5 to +1.5.The key idea of this method is that the model index n can be treated as an additional, continuous parameter that is sampled over in an MCMC chain.At any given step of the evaluation of the hyper-model, the underlying algorithm casts the model index to either 0 or 1, corresponding to one of the two sub-models.Given a current model index n, the hyper-model parameter space is partitioned in an active part, which is used to evaluate the respective likelihood of the sub-model, and an inactive part.The posterior odds ratio P 01 is then be relative amount of chain entries in model 1 compared to model 0, from which the Bayes factor B 01 between the two models can be deduced.
To show how this procedure can be used to calculate a Bayes factor consider the posterior probability for the model index n The first equality used that the posterior for n can be obtained by marginalizing the posterior for n and θ H over the hyper-model parameters.Bayes' theorem leads to the second equality, where Z H is the hyper-model evidence, the first integrand is the likelihood and the second one is the prior for the model parameters and the model index itself.The hyper-model evidence is unknown and difficult to obtain, but of no further importance, as we are interested in the posterior odds ratio P 01 between the two models.For a fixed n we can factorize the prior p(θ H , n|H) in eq.(C.1) as This factorization uses the aforementioned distinction between active (first factor) and inactive (second factor) parameters, which are not correlated with each other as the sub-models are distinct.The last factor is a subjective prior for each of the sub-models.The two last factors do not depend on the active parameters.If we insert these expressions into the definition for the posterior odds ratio for the model index, we find In the last step we used that the inactive parameters, denoted by a bar, do not contribute to the sub-model evidence Z n .A marginalization over their priors therefore gives one for both sub-models and the last factor is equal to one.As P 01 is just the ratio of the number of chain entries after the burn-in period of sub-model 1 compared to sub-model 0 and the model weight ratio Π 01 can be set as a model prior ratio when starting the MCMC chain, the Bayes factor can be obtained by multiplication of the posterior odds ratio with the inverse model weight ratio,

C.2 Implementation details and uncertainties of the computed Bayes factors
The goal of our calculation is to calculate the Bayes factor as accurately as possible.This can be achieved when the posterior odds ratio P 01 is close to one, i.e. if the hyper-model scan spends about the same amount of time in the two sub-chains.We therefore set the model weight ratio Π 01 to the inverse of the expected Bayes factor and iterate over different weight ratios until the posterior odds ratio is close to one.In practice, this iterative procedure is still an intricate problem due to long runtimes.When the two compared models favour rather distinct regions of the available pulsar-intrinsic red noise parameter space, jumping from one sub-model to the other is initially unlikely.Hence, the burn-in period is long and uncertainties in the computed Bayes factors increase.Even though our main aim is to calculate Bayes factors for several DSPT scenarios compared to the no common red noise (nCRN) hypothesis, it is faster and more reliable to first compare a given DSPT scenario with the SMBHB hypothesis, see eq. (3.1) (γ = 13/3, A CP kept as a free parameter).This speeds up the burn-in and results in a more precise computation of the Bayes factor because the posterior distributions for the pulsar-intrinsic red noise parameters are very similar between these two scenarios, unlike when directly comparing to the nCRN hypothesis.As the Bayes factor between a GWB induced by SMBHB and the nCRN hypothesis is known, log 10 B SMBHB/nCRN = 4.5 (9) [3], and since Bayes factors are multiplicative B 02 = B 01 × B 12 , we can simply rescale our results comparing to SMBHB to the Bayes factor compared to the nCRN hypothesis.In doing so, the burn-in period is shortened by a factor of O (10).Still, even with a good choice of Π 01 and using this method, the chains take several days to converge to a Bayes factor.Due to these runtime challenges, a publicly available implementation of the described pilot run in ref. [98] to speed up the calculation of Bayes factors would be highly appreciated.
Note however that the method of first comparing to the SMBHB hypothesis does not work if cosmological constraints do not allow a significant GWB in a DSPT model, i.e. for a stable DS with a large lower boundary of β/H.Here, the DSPT models favour similar regions in pulsar-intrinsic red noise parameter space as the nCRN hypothesis and a direct comparison would be more advantageous, even though we still compare to the SMBHB model first for consistency.The Bayes factors in these cases therefore have larger uncertainties.
The computation of Bayes factors in itself also comes with statistical uncertainties, in particular from the finite length of the underlying chains.We explicitly checked this uncertainty to be under control by computing the Bayes factor as a function of the number of drawn samples.Using 5 × 10 6 samples from the hyper-model (including both sub-models) and conservatively discarding the first 25 % due to burn-in, the Bayes factors all converged to a relative uncertainty of a factor of up to ∼ 2. The convergence rate however depends sensitively on the precise value of the model weights Π 01 as mentioned above.We therefore made sure that the number of samples of both models differed by no more than O(10 %), which required the model weights Π 01 to be known up to one decimal place.
The continuous lines in figure 4, depicting the Bayes factors' expected prior dependence, are further prone to uncertainties due to the reduction of points after adapting the prior, cf.appendix C.4.This is particularly relevant when this reduction is large, i.e. when the Bayes factor is reduced by a significant amount e.g for a GWB from bubble wall collisions and a large lower boundary for β/H.These uncertainties together cause the difference between the expected Bayes factors' prior dependence and the individually calculated Bayes factors for a given prior choice in figure 4.

C.3 Relating Bayes factors to p-values and Z-scores
We can express Bayes factors in terms of a Z-score to describe the improbability of the null hypothesis [131,132].The Z-Score is more commonly known as the "number of sigmas" with which some measured quantity deviates from its expectation value.The probability of the null hypothesis can interpreted to be p(0|data) = 1 − p(1|data) in a frequentist's manner, if one asserts that there is no other possible model to explain the data.If one further interprets the posterior odds ratio to be the ratio of these probabilities P 01 = p(1|data)/p(0|data), the probability of the null hypothesis can equally be expressed as p(0|data) = 1/(1 + P 01 ).This can be interpreted as a p-value, i.e. the probability to measure data as extreme as the one observed if the null hypothesis were indeed correct.Assuming equal prior probabilities for the two models under comparison, i.e.Π 01 = 1, the p-value reads p = 1/(1+B 01 ).We convert the obtained p-value to a Z-score using a one-tailed Gaussian, (C.8) The function Φ −1 is the inverse of the cumulative density function of the standard normal distribution with zero mean and unit standard deviation.For Bayes factors B 01 below 1, implying more evidence for the null hypothesis than for a given more complicated model, we replace B 01 by its inverse.We do so to be able to express the tension between the data sets, where the null hypothesis is a better explanation for the observed data than a more complicated model.In that case the Z-score should then rather be interpreted as to express the probability of obtaining a signal as low as observed if the signal model were instead correct.

C.4 Influence of the prior choice on the Bayes factor
The Bayes factor between two models 1 and 0 is given by the evidence ratio where Z i is the evidence of model i.We denote the model parameters of model 1 by θ, L(θ) is the likelihood to explain the data and π(θ) is the chosen prior.We here want to investigate the effect of a change in prior π(θ) → π(θ) on the Bayes factor B 01 → B 0 1.Keeping the likelihood and its normalization as well as the number of model parameters unchanged, the Bayes factor differs by a factor (C.10) approach spans over four orders of magnitude, just as the prior for A SMBHB ∈ 10 −18 , 10 −14 , which allows for an easy comparison of the two models on a similar footing.
In the first analysis we obtained a Bayes factor of O(5) between the decaying DSPT and SMBHB interpretation in favour of the DSPT hypothesis.The second analysis found a Bayes factor of O(1.4) in favour of the decaying DSPT interpretation.We checked explicitly that this matches the expected Bayes factors from an "a posteriori" prior change as described in appendix C.4.Both Bayes factors are also consistent with each other, since the one of the second analysis could be enhanced by a factor of 4 if one decreased the prior range to one decade α tot ∈ 10 −2 , 10 −1 , confirming the O(5) Bayes factor of the first analysis.To make this a fair comparison it should however be noted that also the evidence of the SMBHB hypothesis can be boosted by a relative factor of 4 by decreasing the prior range of A SMBHB from 10 −18 , 10 −14 to 10 −15 , 10 −14 , cf. figure 7 in appendix C. 4.
We therefore conclude that a GWB with a broken power-law spectrum with fixed slopes Ω sw (f ) ∝ f 3 and Ω sw (f ) ∝ f −4 for low and high frequencies respectively, can in principle fit the common red spectrum slightly better than a featureless single power-law with spectral index γ SMBHB = 13/3.The origin of the discrepancy between the highest Bayes factors depicted in figure 4 and the reference value of 10 4.5 is due to prior volume effects.These rightfully decrease the credibility of the phase transition explanation, as the required parameter space to explain the common red spectrum is much bigger than that of the SMBHB interpretation.
In particular, the difference in Bayes factors is not due to a worse fit of the pulsar timing data.In fact, in the considered parameter ranges, a decaying DSPT would yield a comparably good (or even slightly better) fit to the timing residuals.If future data is fitted even better by the DSPT explanation compared to SMBHBs, the prior volume effects could be alleviated.With just the NANOGrav 12.5 yr data set the SMBHB explanation is favoured.

Figure 1 :
Figure 1: Triangle plot showing the 1 σ and 2 σ contours obtained by a naive fit (blue)of the NANOGrav 12.5 year data to a GW spectrum emitted in a DS phase transition, ignoring cosmological constraints.To illustrate the tension with BBN and CMB, the 95% C.L. excluded regions corresponding to ∆N eff > 0.22 are shaded in red[37], cf. the discussion in section 4. The regions in which no percolation (β/H < 3) and an overestimation of the GWB amplitude (β/H < 10) are expected are shaded in grey, see section 2 for further details.
Figure2: Plot showing the NANOGrav "violins" (red)[3], their standard explanation through a power-law spectrum from the inspiral of SMBHBs with A SMBHB = 1.53 × 10 −15 and γ SMBHB = 13/3 (grey) and two phase transition spectra.A characteristic bubble wall collision spectrum is shown in orange.A sound wave-induced GWB spectrum is shown in blue.These spectra correspond to the best-fit parameter points found in section 5 (including cosmological constraints and demanding β/H > 1).The blue arrows indicate how an increase in the phase transition parameters α tot , T perc SM or β/H would shift the spectra.

Figure 3 :
Figure3: Comparison of the triangle plots for two MCMC chains with different priors on β/H, assuming a stable DS and a GWB generated through sound waves.Blue shaded regions mark 1 σ and 2 σ contours, respectively.Regions with β/H < 3 and β/H < 10, in which the phase transition cannot complete and the GW signal is overestimated respectively, are indicated with grey shadings.

3
Bayes factorBayes factors for a phase transition vs.only pulsar-intrinsic red noise Sound waves, decaying dark sector Sound waves, stable dark sector Bubble wall collisions, stable dark sector

Figure 4 :
Figure4: Bayes factor estimates of various DSPT models with respect to the no common red noise-hypothesis.Filled coloured dots correspond to actually performed full model comparisons, while lines in the corresponding colours are derived from an a posteriori reduction of the prior range of chains with minimal β/H of 1 (for details, see appendix C.4).We also indicate how to translate the Bayes factor to Jeffrey's scale (coloured shadings) as well as Z-scores from frequentist statistics (right y-axis, cf.appendix C.3).In the grey shaded area (β/H < 3), a constant nucleation rate is not sufficient to drive percolation.
the transition strength is defined by α * = ∆θ/ρ perc SM such that ρ reh SM = (1 + α * ) ρ perc SM and we need to map α * = α tot + (1 + α tot ) b le w a ll s , s t a b le D S , β / H > 1 B u b b l e w a l l s , s t a b l e D S ,

Figure 6 :
Figure6: Posterior distributions of GWB spectra for parameter scans assuming bubble collision spectra from a stable DS phase transition with β/H above 1 (orange) and above 10 (yellow).Due to cosmological constraints, the spectra with β/H > 10 clearly cannot fit the common red spectrum detected by NANOGrav[3] (cf.red violins).If the dark sector can decay (violet spectra), slightly faster transitions with β/H > 10 can still explain the observed signal.The GWB from SMBHB with A SMBHB = 1.53 × 10 −15[3] is depicted as a grey line for comparison. 59] In the last step we used that quite generally S 3 /T DS is a polynomial in T DS , which implies that β/H is naturally of the same order as S 3 /T DS at nucleation.The first bubble in a Hubble volume nucleates when Γ * ≃ H 4 * .Using eq.(2.9) one finds that (S 3 /T DS ) | t * ≃ log A/H 4 * .During radiation domination the Hubble parameter is of the order H 2