Introduction.-The direct detection of gravitational waves (GWs) by the LIGO and Virgo collaborations [1] has opened up a fascinating era of astronomy and cosmology that looks at our Universe through entirely new eyes. Since GWs propagate without interaction, their detection enables us to probe physics in the early Universe. While ground-based interferometers such as LIGO and Virgo have the best sensitivity at frequencies of around 100 Hz, searches for GWs with lower frequencies have been also conducted or planed in various types of experiments. Future space-based GW observers such as LISA [2] , DECIGO [3] , and BBO [4, 5] have the best sensitivity at frequencies of order mHz. GWs with even lower frequencies of O(10 −9 ) Hz are searched for by pulsar timing array (PTA) experiments such as EPTA [6] , PPTA [7] , and NANOGrav [8, 9] . Now, the discovery of such low frequency GWs may be right around the corner.
Recently, the NANOGrav collaboration of a PTA experiment has analyzed their 12.5 years of data and reported a signal that may be interpreted as a GW background [10] . One possible source of the signal is an astrophysical GW background generated by mergers of supermassive black-hole binaries [11] [12] [13] [14] . Another possibility is a stochastic GW background emitted by a cosmicstring network in the early Universe, which can give a favored flat spectrum of frequencies in the GW energy density [15, 16] (see also, e.g., Refs. [17] [18] [19] [20] [21] [22] [23] [24] [25] for earlier works). A stochastic GW signal associated to primordial black hole formation has also been investigated [26, 27] .
In this paper, we explore an interpretation of the reported NANOGrav signal in terms of a stochastic GW background generated by a new strongly first-order phase transition which occurred in a dark sector. Generation of GWs from a strongly first-order phase transition has been actively discussed. A possible detection of GWs from the QCD phase transition through PTAs was pointed out in Ref. [28] . GWs from a supercooled electroweak phase transition and their detection with PTAs were discussed in Ref. [29] (see also Ref. [30] ). Such a supercooled phase transition has been known to be realized in warped extra dimension models or their holographic duals [31] [32] [33] [34] [35] [36] [37] [38] [39] . GWs from dark sector phase transitions were explored in Refs. [40] [41] [42] [43] [44] though they considered the case in which the dark sector is not completely decoupled from the visible sector. In particular, Ref. [40] discussed a range of GW frequencies covered by PTAs. In this paper, we focus on the first order phase transition in the dark sector decoupled from the visible sector. The production of GW in such a decoupled dark sector was studied in Refs. [45, 46] . As we shall see below, however, the estimated GW amplitude in Ref. [45] was underestimated, and we will also use the updated GW spectrum taking account of the effective lifetime of the GW sources. We derive a consistency relation between the GW energy density to explain the reported NANOGrav signal and the allowed abundance of dark radiation components in a dark sector. If the NANOGrav signal is confirmed, it predicts a future discovery of dark radiation components in our Universe. The amount of dark radiation will provide us with information on details of the phase transition. Interestingly, in some (rather realistic) parameter space our scenario predicts ∆N eff ∼ 0.4 -0.5 which can ameliorate the socalled Hubble tension [47] [48] [49] .
The rest of the paper is organized as follows. In the next section, we derive the relation between the amount of dark radiation and the amplitude of stochastic GWs arXiv:2009.09754v2 [astro-ph.CO] 27 Sep 2020 generated by a first-order phase transition in a dark sector. Then we compare the NANOGrav data with the GWs generated by the phase transition and discuss its implications for the amount of dark radiation. Finally, we conclude and discuss a possible model of the dark sector. In Appendix, we quote and summarize the GW spectrum generated from the first order phase transition that is used in our numerical calculations. GW and dark radiation.-We consider the case in which the GW signal reported by NANOGrav comes from the first-order phase transition in the dark sector. As we stated in the introduction, we assume that the dark sector is coupled to the visible sector only via the gravitational interaction and study its implications for dark radiation. We do not specify how the first-order phase transition occurs in the dark sector, but describe nature of the phase transition in terms of phenomenological parameters such as the wall velocity v w , duration of the phase transition β −1 , and an efficiency factor κ i defined later and/or in Appendix.
As the dark sector is decoupled from the visible sector, it is reasonable to assume that the entropy is separately conserved in each sector until the beginning of the phase transition and after the end of the phase transition. We denote the entropy ratio between the two sectors before the phase transition as R i :
where s D (s vis ) is the entropy density in the dark (visible) sector and T * is the temperature of the visible sector at the time of the first order phase transition. As we will see, we are interested in the case in which the dark sector never dominates the energy density of the Universe, so that the entropy production does not dilute particles in the visible sector much, including the baryon asymmetry. The ratio R i is considered to be determined by the reheating process after inflation. We take it as a free parameter, but expect it to be neither very small nor very large, as in the case of the universal reheating. When the phase transition is a strong first-order phase transition, a large entropy is generated from the latent heat in the dark sector. Denoting the entropy production factor by ∆(≥ 1), we can express the entropy ratio after the phase transition as
If one assumes that the phase transition is not a strongly supercooled one, R ∼ R i and is neither very small nor very large.
After the phase transition, some amount of energy (and entropy) should remain in the dark sector that is decoupled from the visible sector. Then, the lightest particle(s) in the dark sector should be (almost) massless and behave as dark radiation, since otherwise the abundance of the remnant in the dark sector would easily exceed the observed dark matter density and overclose the Universe. 1 The amount of dark radiation at the recombination epoch, ρ DR,0 , is calculated as
where g * 0 (g * s0 ) and g (D) * 0 (g (D) * s0 ) are the effective numbers of relativistic degrees of freedom for the energy (entropy) densities in the visible sector and the dark sector, respectively, and ρ rad,0 is the energy density of photons and three neutrinos in the visible sector. The subscript 0 represents the value at the recombination epoch. This gives the extra effective neutrino number of
The Planck data combined with the BAO observation and the local Hubble measurement gives the constraint [48, 50] 
The prediction in the standard cosmology is N (std) eff = 3.046. Note that there is a tension between the local measurement of the Hubble parameter and the Hubble parameter inferred by the Planck and BAO with ∆N eff = 0. The tension is known to be relaxed if ∆N eff 0.4 -0.5 [47] [48] [49] , which roughly amounts to R 0.1 -0.2. Now we shall relate the entropy ratio and the density parameter of the GW. We define
where ρ vac is the false vacuum energy in the dark sector at the phase transition and ρ rad,
is the total radiation energy density just before the phase transition. Here, we denote T (D) * i and T (D) * f as the temperatures of the dark radiation before and after the phase transition, respectively. After the phase transition, the vacuum energy is converted into the radiation energy in the dark sector. Using the entropy production ∆, we can relate α and R as
where we assume the instantaneous reheating after the phase transition and use ρ vac + ρ DR (T
We are interested in the case in which the GW is efficiently emitted. From this consideration, we assume that the phase transition is a strong first-order phase transition, which implies ∆ O(1). In this case, we can neglect the second parenthesis in Eq. (7). We also note that R 0.1, which leads to α r 1. Combining Eqs. (4) and (7), we obtain
where g * 's are evaluated at T = T * or T (D) * f like Eq. (8). We also define κ i that generically represents the fraction of latent heat converted to the GW source labeled by i, which we specify later and in Appendix. Since the GW comes from the quadrupole moment, the resulting GW amplitude is proportional to the energy density of the source squared. We also note that the amplitude of stochastic GW is involved with two time integrals. Motivated by this observation and taking into account the redshift factor, we factorize the density parameter of the GW at present as (see, e.g., Refs. [51] [52] [53] )
where H * ≡ H(T * ) and Ω rad,0 h 2 = 4.16 × 10 −5 [54] with h being the reduced Hubble parameter. Here, β −1 represents the duration of the phase transition, defined by
with Γ being the bubble nucleation rate. The remaining factorΩ GW,i is determined by numerical simulations and/or (semi-)analytic calculations (see Appendix). In particular, it depends only on k/β and the bubble wall velocity v w for the GW emission from the bubble collision under a certain assumption [53] . Using Eq. (10), we can see that the density parameter of the GW is proportional to ∆N 2 eff . These are the formula that relates the amplitude of the GW to the amount of dark radiation ∆N eff . This can be applied to any models where the GWs are emitted in a dark sector that is completely decoupled from the visible sector except via the gravitational interaction. Contrary to the ordinary scenario where the phase transition occurs in the visible sector, the amplitude of the GW is determined by α rather than α. Here α cannot be much larger than the order 0.1 due to the constraint on the dark radiation (see Eq. (5)). GW and NANOGrav Pulsar Timing Data.-There are three possible sources of the GWs from the first-order phase transition; the collision of true-vacuum bubbles [51, 53, [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] , the sound waves [67] [68] [69] [70] , and the turbulence [51, [71] [72] [73] [74] [75] . For each source, we define κ i with i = bubble, SW, turb. We summarize the GW spectra produced from these sources in Appendix.
Note that κ i are determined by the dynamics in the dark sector, so that they are related to parameters in the dark sector such as α rather than α . 2 Since α can be larger than of order unity, some κ i can be as large as of order unity [51, 76] . To compare the NANOGrav data, we provide the GW spectra for two cases using the formula written in Appendix. The first case is the GW spectrum only from the sound wave and turbulence while the second case is the one from all three sources. We take κ turb = 0.1κ SW for both cases, as suggested by numerical simulations [77] . We assume the maximal efficiency, such as κ SW = 1/1.1 in the former case and κ bubble = κ SW = 1/2.1 for the latter case.
In our numerical calculations, we take v w 1 and (α/(1 + α)) 1 and set g (D) * = g (D) * 0 and g (D) * s = g (D) * s0 for simplicity. The duration of the phase transition, β/H * , is typically about 100 for T * ∼ 1−100 MeV though it can be as small as of order unity depending on models [78] . We take it to be a free parameter within (1 -100) to show examples. Figure 1 shows the GW spectra produced by the firstorder phase transition and the NANOGrav 12.5 year results with older experimental constraints. In the upper figure, we show the GW signals from the sound wave and turbulence with β/H * = 5 and T * = 5.6 MeV. In the lower figure, we include all three sources with β/H * = 2 and T * = 75 MeV. The red lines represent the GW spectrum for several values of ∆N eff . The double-headed arrow shows the amplitude favored by the NANOGrav 12.5 year signal within 2-σ posterior contour when the spectrum is assumed to be flat [10] . The blue, yellow and green shaded regions are the previous bounds from EPTA [6] , NANOGrav(11yr) [79] , and PPTA [80] . While the results of the NANOGrav 12.5yr are not compatible with the older constraints, the tension is understood by the improvement of the pulsar red noise treatment in the NANOGrav 12.5yr [10] .
In Ref. [10] , they fitted the NANOGrav pulsar timing data by a power law spectrum and a broken power law spectrum and concluded that the data favor spectrum with a power of −0.5 -1.5 around f = 3.7 × 10 −9 . However, the most relevant data are two bins at the first and second lowest frequencies, which can be fitted by a spectrum with a peak at f 3.7 × 10 −9 , just like the ones in the upper figure in Fig. 1 . The spectrum shown in the lower figure has a slightly positive power in the relevant frequency scale, which is also favored by the NANOGrav data. We also plot β/H * on the T * -∆N eff plane in Fig. 2 . The upper figure corresponds to the case for the GW produced by the sound wave and turbulence while the lower one corresponds to the case for all three sources. The yellow, blue, and green lines correspond to the case with the peak amplitude of Ω GW (f )h 2 = 2 × 10 −9 , 6 × 10 −10 , 2 × 10 −10 for the peak frequency of f = 3.7 × 10 −9 Hz, respectively, which are implied by the NANOGrav data. In more details, in the lower figure, we require the peak of the bubble collision corresponds to the NANOGrav data implied values. However, we note that the NANOGrav data can also be explained by the highest peak of the spectrum that comes from the sound wave (rather than the one comes from the bubble collision). In this case, the results are almost the same as the one in the upper figure even if we include all three sources.
The numbers shown near the dots on the lines represent β/H * . We note that β/H * is typically as large as about 100 for realistic cases but can be as small as O(1) for the case of, e.g., supercooling phase transition. The gray-shaded region is excluded by the Planck and BAO observations while the red-shaded region is favored to ameliorate the Hubble tension. The black-dotted line in the upper figure is the prospected 2σ sensitivity of the stage-IV ground-based detector, CMB-S4 [81] . Our scenario predicts ∆N eff that can be measured in the near future, by explaining the NANOGrav data by the phase transition in the dark sector. In fact, some parameter space is also favored to ameliorate the Hubble tension.
Discussions and conclusions.-We have discussed the implications of NANOGrav data for the case in which the stochastic GWs are produced from the first order phase transition in the dark sector, assuming that the dark sector is completely decoupled from the visible sector except via the gravitational interaction. Since there must be a remnant in the dark sector after the phase transition, the GW amplitude is related to the amount of dark radiation. Interestingly, the predicted abundance of dark radiation is within the reach of CMB-S4 in the near future. In some parameter space, the amount of dark radiation is as large as the one preferred to ameliorate the Hubble tension. In other words, the signal observed by NANOGrav is closely related to the Hubble tension.
If the dark radiation is not observed in the near future, it may imply that the dark sector is coupled to the visible sector. In this case, T * should be larger than of order 1 -10 MeV so that the phase transition and the subsequent energy injection to the visible sector should not spoil the success of the Big Bang nucleosynthesis. Note that the Big Bang nucleosynthesis bound is significantly weaker in the case of the dark sector completely decoupled from the visible sector.
The peak frequency of the GW is related to the energy scale of the phase transition, which we found is T * = O(1 -100) MeV. As the energy scale is close to the QCD scale, one may think of a possibility that the GW signal comes from the phase transition of a dark QCD. This is the case, e.g., of a parallel world, where the dark sector has a similar (but not exactly the same) structure to the visible sector [82] . If the parameters in the dark sector are slightly different from the ones in the visible sector, it is possible that the dark-QCD phase transition is the strong first-order phase transition and its energy scale is T * = O(1 -100) MeV. This scenario also provides a dark-matter candidate, which is the dark neutron. Since the dark neutron can have a sizable self-interaction cross section, this scenario might be confirmed by astrophysical observations for the dark-matter density profile. Also the dark neutrinos behave as hot dark matter, which will mitigate the σ 8 tension when dark radiation is introduced to solve the Hubble tension.
So far we have focused on the GWs generated from a first order phase transition in the dark sector, but it is also possible to generate similar GW spectra by considering decays of topological defects such as domain walls and/or cosmic strings (see e.g. Ref. [83, 84] ). Our consistency relation can easily be extended to these cases as well. Also, the clockwork QCD axion [85] is known to have an extremely complicated network of strings and domain walls which store a large amount of energy [86, 87] . The string/wall network annihilate around the QCD phase transition, which may produce sizable GWs with frequencies covered by the PTA experiments [86] . Appendix: GW spectrum.-In this appendix, we quote and summarize explicit values ofΩ GW (see, e.g., Refs. [77, 88] for more detailed discussion).
There are three possible sources of the GWs from the first-order phase transition; the collision of true-vacuum bubbles [51, 53, [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] , the sound waves [67] [68] [69] [70] , and the turbulence [51, [71] [72] [73] [74] [75] . The parameter κ i are given by
where ρ bubble , ρ SW , and ρ turb are the energy densities of a thin shell around the bubble wall, bulk motion of the fluid, and turbulence, respectively. Here, ρ bubble should be evaluated just before the end of the phase transition while ρ SW and ρ turb should be evaluated just after the phase transition. As we stated in the main part of this paper,Ω GW is a function of k/β and v w for the GW emission from the bubble collision. On the other hand, it includes a factor of β/H * and other time scales for the GW emission from the sound wave and magnetohydrodynamics turbulence because the duration of the GW emission from these sources is relatively long [72] . It is enhanced by a factor of (β/H * )(1 − 1/ √ 1 + 2t sw H * ) for the GW emission from the sound wave, where t sw (8π) 1/3 v w /(βU f ) is the sound-wave period [89] [90] [91] [92] and U 2 f (3/4)κ sw α/(1 + α) is the root-mean-square fourvelocity of the plasma [69, 70] . In our numerical calculation, we take U 2 f = κ sw . On another hand, it is enhanced by a factor of (β/H * )(ρ DR /ρ turb ) 1/2 for the emission from the turbulence [75] .
The GW spectra produced by the bubble collision [60] , sound waves [69, 93] , and turbulence [72, 94] 3 are given bỹ
respectively, where
F sw (f ) = f f sw
F turb (f ) = (f /f turb ) 3 (1 + (f /f turb )) 11/3 (1 + 8πf /h * )
.
Here, f is the frequency today and h * is the Hubble parameter at T = T * that is redshifted today: 
The peak frequencies f bubble , f sw , and f turb are given by 
3 SinceΩ GW,turb is the one for the magnetohydrodynamics turbulence, we implicitly assume that the dark sector has a gauge interaction. As discussed in Ref. [69] , however, this contribution is subdominant and this assumption does not change our result.
If the vacuum bubble interacts with the thermal plasma strongly enough and if the energy density of the thermal plasma in the dark sector is large enough, the accelerating bubble wall receives a large friction from the interaction with the thermal plasma [95] . As a result, the bubble wall velocity may reach a terminal velocity to balance between the friction effect and the pressure due to the false-vacuum energy. Then most of the kinetic energy of the accelerating bubble wall is injected into the thermal bath. In this case, sound wave and turbulence of the plasma are the main sources of GWs. On the other hand, if the interaction between the vacuum bubble and the thermal plasma is not strong enough or the energy density of the thermal plasma in the dark sector is not large enough, the bulk energy is negligible and the bubble collision is the dominant source of GWs [77] . This is the case for e.g., a significant supercooling (i.e., the case of ∆ 1). In the intermediate case, all three sources are relevant.