Strong gravitational radiation from a simple dark matter model

A rather minimal possibility is that dark matter consists of the gauge bosons of a spontaneously broken symmetry. Here we explore the possibility of detecting the gravitational waves produced by the phase transition associated with such breaking. Concretely, we focus on the scenario based on an SU(2)D group and argue that it is a case study for the sensitivity of future gravitational wave observatories to phase transitions associated with dark matter. This is because there are few parameters and those fixing the relic density also determine the effective potential establishing the strength of the phase transition. Particularly promising for LISA and even the Einstein Telescope is the super-cool dark matter regime, with DM masses above O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(100) TeV, for which we find that the gravitational wave signal is notably strong. In our analysis, we include the effect of astrophysical foregrounds, which are often ignored in the context of phase transitions.


Introduction
Cosmological and astrophysical observations strongly suggest that, in contrast to the ordinary substances found on Earth, baryons are not the dominant constituent of the matter in the Universe [1]. Such non-baryonic matter is called dark because its interactions with the Standard Model (SM) particles -particularly with photons -are constrained to be very weak. This, along with the obvious fact that dark matter (DM) must be stable on cosmological timescales, are the two most important properties of any DM candidate.
The first property is often invoked as an argument for the electroweak (EW) nature of DM interactions. In fact, models where DM is directly coupled to the W or Z bosons naturally explain the DM relic density by means of the thermal freeze-out of DM annihilations in the Early Universe. Nevertheless, these scenarios have been dramatically constrained JHEP05(2019)190 in the past couple of decades by direct and indirect detection experiments, together with colliders, most recently the LHC [2][3][4]. In contrast, models where DM is directly coupled to the Higgs and not to the W or Z bosons are much less constrained by the aforementioned experiments, especially in regimes where the DM is heavier than the Higgs boson. Interestingly, gravitational waves (GWs) offer a new complementary way to probe the latter scenarios. This is because they typically require the existence of additional scalar fields, which can potentially trigger a first-order phase transition (PT) in the Early Universe and therefore the emission of GWs [5][6][7][8]. Of course, DM can be probed in this way only if its properties are closely related to the PT [9][10][11][12]. This is the subject of the present work.
In order to motivate a concrete choice for the DM model, we will invoke the second DM property mentioned above, i.e. its stability. This is often ensured by imposing a discrete symmetry in the DM sector. The most common examples being Z 2 symmetries or the socalled R-parity in supersymmetric theories. Nevertheless, such symmetries are not known to exist in nature. 1 Better-motivated scenarios are those where DM is stable as a result of its own dynamics. In fact, this is exactly what happens with the stable particles of the SM. For instance, proton stability follows from baryon number conservation, which is an accidental symmetry due to the SU(3) C × SU(2) L × U(1) Y charges of the matter fields. An incomplete list of examples of this type of scenarios include Minimal DM [13], spin-one DM models [14][15][16][17] and QCD-like models of DM [18].
The previous observations motivate us to study the spin-one DM model proposed in [14], in which the DM portal to the SM is the Higgs boson. Concretely, we extend the SM with a dark SU(2) D local symmetry, under which all the SM particles are assumed to be singlets. In addition, we postulate a dark scalar doublet which carries no SM charges and whose vacuum expectation value (VEV) breaks the SU(2) D symmetry via a Higgs mechanism in the dark sector, ensuring the theoretical consistency of the model containing massive spin-one fields. After symmetry breaking, the particle content includes -besides the SM -three mass-degenerate particles of spin-one and one dark Higgs boson. In this model there is a custodial SO(3) symmetry remaining in the broken phase, under which the gauge bosons transform, ensuring their stability. Collectively, these comprise our DM candidate, which only couples to itself, to the SM Higgs h, and to the dark Higgs h D . The Higgs portal interaction allows h D to decay to light SM particles, thus avoiding it becoming a DM component.
We will consider two production regimes for the DM relic density. First, the standard thermal freeze-out of DM annihilations into dark Higgs bosons. Second, super-cool DM [19], a more exotic possibility in which we assume a classically scale invariant potential for our model [20][21][22][23]. As pointed out recently, this can result in a period of late-time inflation which sets the relic density in a completely novel way. In both cases, a PT takes place in the early Universe from a SU(2) D symmetric vacuum in which the would-be-DM is massless, to a vacuum in which the dark gauge symmetry is broken and the DM is massive. The key point of our analysis is that the parameters setting the relic density also enter the effective potential determining the PT. As we will see, this allows us to find correlations between the GW signal and the DM properties.

JHEP05(2019)190
This study is timely, as much work is being done on understanding GWs from cosmological PTs in anticipation of LISA [8], and follow-up proposals such as BBO [24]. Our analysis differs from recent similar works in at least three aspects. First, as already mentioned, our scenario is rather minimal, with only four parameters in the general case and two for super-cool DM. This allows us to establish a close connection between the emission of GWs and the relic density or direct detection. Second, in our analysis astrophysical foregrounds will be taken into account. These are mostly due to binaries of white dwarfs and are crucial for estimating the signal-to-noise ratio at future GW observatories [6,7]. Finally, we discuss for the first time the GW signatures of the super-cool DM regime. The paper is organized as follows. In section 2, we present our DM model and its phenomenology. In section 3, we calculate the GW signal arising from the PT for the standard and the classically scale invariant cases. We conclude in section 4 by presenting a summary and outlook for this work. Appendix A is devoted to details concerning the effective potential, which determines the nature of the PT, appendix B summarises the contributions to the GW spectra, appendix C includes some additional material regarding the classically scale invariant potential, and appendix D discusses bubble percolation in the vacuum dominated regime.

The model
In this section we will describe the model and define notation. As mentioned in the introduction, we consider an extension of the SM with a dark SU(2) D gauge symmetry [14], under which all the SM particles are singlets. In addition to the dark gauge bosons A i Dµ (i = 1, 2, 3), the model has a dark scalar doublet, H D , which carries no SM charges. Hence, the Lagrangian of the model is 2 and H is the SM scalar doublet. Here, F D is the field strength tensor of the SU(2) D gauge symmetry and D = ∂ + ig D τ i · A i D /2 is the corresponding covariant derivative. We write scalar doublets as where φ and η are the classical field values breaking the EW and the SU(2) D symmetries, respectively. In addition, h, h D , G i and G i D (i = 1, 2, 3) are the corresponding Higgs and Goldstone boson fields.
Symmetry breaking at tree level. In this case, the minimum of the potential associated with eq. (2.1) is located at (φ, η) = (v φ , v η ), where v φ = 246 GeV and

JHEP05(2019)190
The mixing of the real scalars is captured by the usual angle This is constrained by the Higgs signal strength measurements, |θ| O(0.1) [25,26], with the precise limit depending on which combination of measurements is taken. For convenience we commit a small abuse of notation, and from now on also label the mass eigenstates with h and h D , where m h = 125 GeV. The mass eigenvalues are given by All the dark gauge bosons obtain the mass, m A = g D v η /2. In fact, they transform as a triplet under a custodial SO(3) symmetry. Notice the presence of light fermionic fields transforming under SU(2) D would spoil the stability of the vector DM, allowing the gauge bosons to decay, as occurs in the SM [14]. The absence of such fermions allows the model to remain rather minimal with only four parameters in the DM sector, which we take as m A , g D , θ and m h D .
Radiatively-induced symmetry breaking. An alternative possibility is to consider a classically scale invariant realisation of this model [19][20][21][22][23], where the mass terms in eq. (2.1) are forbidden and symmetry breaking is achieved through radiative effects. This is known as the Coleman-Weinberg mechanism [27,28]. A systematic analysis of radiative symmetry breaking with the above field content can be found, e.g. in [23]. In the present analysis, the parameter regime of interest corresponds to what has been termed sequential symmetry breaking [23]. The running of λ 2 results in it turning negative in the IR, breaking the SU(2) D symmetry via the Coleman-Weinberg mechanism. 2 If λ 3 < 0, the breaking of EW symmetry follows sequentially from the induced tachyonic mass λ 3 v 2 η /4, which leads to v φ = v η −λ 3 /(2λ 1 ). Since we are interested in DM above the EW scale, i.e. v φ v η , the magnitude of the portal must be very small, |λ 3 | 1. This implies that the approximation of ignoring the φ direction in studying the SU(2) D symmetry breaking is consistent. Under these assumptions, we can study SU(2) D symmetry breaking by focusing on the term L ⊃ −λ 2 η 4 /4, where the coupling λ 2 is evaluated at a sliding scale given by the value of the η field, giving with η 0 being the scale at which λ 2 flips sign. Here, we neglect the contributions of λ 2 and λ 3 to R.H.S. of eq. (2.7), which is a valid approximation provided λ 2 g 4 D and λ 2 3 g 4 D [23]. We also ignore the running of g D . Note the dark and visible sectors are close to decoupled not only because the portal coupling is small but also because the corresponding beta Figure 1. The dominant DM annihilation channels for m A m h D and θ 1.

JHEP05(2019)190
function is proportional to λ 3 . In fact, the running of the latter between v η and v φ is not so large as to affect our analysis. 3 As alluded to above, λ 2 < 0 signals the breaking of the SU(2) D gauge symmetry. In fact, the minimization conditions to leading order in λ 3 , give v η = η 0 e −1/4 together with As in the previous case, the dark gauge bosons obtain a mass m A = g D v η /2. Notice also that, after accounting for m h = 125 GeV and v φ = 246 GeV, there are only two free parameters, which we choose as m A and g D . Before discussing DM production, we would like to emphasize that this scenario is not simply a limit of the previous case when µ 1 and µ 2 approach zero because here the breaking of the symmetry does not occur at tree level.
(For a detailed discussion on such a limit, see [27].)

Relic density
We will consider two production regimes for the DM relic density: the standard freeze-out scenario and super-cool DM. 4 The latter only takes place for the classically scale invariant case, i.e. when the gauge symmetry is broken radiatively. Details are given in section 3.3.
For the former case, we make the mild assumption that λ 3 and g D are large enough so that DM was in thermal equilibrium with the SM fields in the Early Universe. Freeze-out leads to the observed dark matter abundance, Ωh 2 0.12, when the corresponding cross section is of the order 2.3 × 10 −26 cm 3 /s. This means that for given m A , m h D , and θ, the relic density fixes the dark coupling g D . We are interested in the regime in which m A > 2m h D so that DM (semi-)annihilates into dark Higgs bosons. We make the further simplifying assumption, m A m h D and θ 1, so that the annihilations into SM particles by means of a scalar exchange in the s-channel are negligible and the dominant annihilation channels are those shown in figure 1. In this regime the correct relic density is achieved for, A more accurate determination can be achieved by numerically solving the Boltzmann equations. Given the uncertainties of the gravitational wave spectrum, however, the use 3 One can make this statement more precise by considering a scalar potential improved with Renormalization-Group effects, as recently suggested in refs. [23,29]. This can cure any potential pitfall associated to the disparity of VEVs in the scalar potential. However, such analysis lies beyond the scope of this work. 4 Other production mechanism for this model have been discussed in [30]. of eq. (2.9) is sufficient for our purposes. The coupling g D is fixed by the relic abundance, which effectively collapses the higher dimensional parameter space to three (one) dimensions in the standard (classically scale invariant) case, which would otherwise have to be scanned over.

Direct detection
The spin independent scattering cross-section of dark matter off nucleons is [14] where m N denotes the nucleon mass and f 0.3 is a constant that depends on the nucleon matrix element. Thus, current experiments such as Xenon1T [31,33,34] constrain the mixing angle between the scalars. This is shown in figure 2 assuming m A 0.1 TeV together with the corresponding future sensitivity from LUX-Zepelin [32,[35][36][37]. In the left panel, the gauge coupling is fixed by freeze out. In the right panel, we consider radiativelyinduced symmetry breaking and therefore we use eq. (2.8) to calculate the cross section. No production mechanism is assumed in the latter case.
The left panel of the figure can be understood as follows. Working in the limit m A m h D m h and substituting eq. (2.9) into σ SI , one finds σ SI ∝ m A sin 2 2θ, i.e. with no further dependence on the other DM parameters. Then comparing to direct detection constraints -which also scale as σ limit SI ∝ m DM for m DM 100 GeV -we find current experiments such as Xenon1T demand θ O(0.1). Future experiments such as LUX-Zepelin, JHEP05(2019)190 will improve σ limit SI by two orders of magnitude and therefore probe down to θ ∼ O(0.01). On the other panel, eq. (2.8) indicates that σ SI ∝ 1/m 6 A with no further dependence on the other dark sector parameters. Hence, σ limit SI ∝ m DM excludes contours of constant DM mass. In fact, for the classically invariant case, the current direct detection constraint demands m A 0.9 TeV [19].

Calculation of the spectrum
A first-order PT takes place by means of the nucleation of true-vacuum bubbles in a falsevacuum background. At a given temperature T , the rate per unit of volume at which this occurs scales as T 4 e −S , where S is the Euclidean action evaluated at the solution describing one bubble. In practice, S can be approximated by the smallest of S 3 /T or S 4 , where S n is the O(n) symmetric action. For our model, we have checked 5 that S 3 /T < S 4 and therefore Here, r is the radial coordinate of the bubble, V (φ, η, T ) is the finite-temperature effective potential, which we calculate using the well-known techniques of thermal field theory, and (φ 0 , η 0 ) are the field values of the false vacuum. The thermal functions which enter V (φ, η, T ) are evaluated numerically. Further details are given in appendix A. To calculate S 3 , we use the over/under shooting method implemented in our code to find the bubble profile by solving the equations of motion, with the boundary conditions dφ/dr r=0 = dη/dr r=0 = 0 and φ| r→∞ = φ 0 and η| r→∞ = η 0 . Physically, these conditions correspond to demanding a smooth profile at the centre of the bubble, r = 0, together with the Universe being in the false vacuum well outside of the bubble, r → ∞. In our discussion here we have remained general by including both fields, however, in the PTs studied below only the η field value will be changing which simplifies our calculation of the action.
Nucleation occurs at a temperature T n , when the bubble nucleation rate in the horizon volume becomes comparable to the Hubble parameter, from now on denoted H. Hence,

JHEP05(2019)190
we find the nucleation temperature by solving [38] where ρ vac is the vacuum energy density, and g * counts the effective radiation degrees of freedom. With the field content of the present model, g * = 116.75 when all species are relativistic and thermalised. For the parameter regions of interest in this work, we find that the nucleation temperature is significantly different from the one associated with EW symmetry breaking. As a result, the phase transition -and therefore the GW emission -is mostly determined by the one-dimensional η direction.
An isolated spherical bubble does not radiate gravitationally because its quadrapole moment is zero. In contrast, the collision of several bubbles generates GWs by means of at least three different processes: collision of bubble walls (mostly determined by the dynamics of the scalar fields), bubble percolation producing sound waves, and magnetohydrodynamic turbulence in the plasma. These give rise to where ρ c is the critical density, and dρ GW /df is the differential GW energy density. Its determination is an active area of research with a large number of ongoing investigations. This is illustrated by the fact that the study of GWs from spin-one DM that was briefly discussed in ref. [20] only accounted for the collision of bubble walls. Nevertheless, as discussed at length below, recent developments [39] when applied to our model indicate that sounds waves and turbulence give the dominant contribution, at least for the case when symmetry breaking occurs at tree level.
In the light of this, here we use the compendium of results presented in ref. [8] and summarized in appendix B in finding the GW spectrum. The spectrum depends on four parameters: the Hubble parameter at nucleation, the wall velocity, v w , the latent heat (here normalised to the radiation density ρ rad ), where (φ n , η n ) is the true vacuum at T n , and the timescale of the transition, (3.7) Figure 3. An example of the gravitational wave spectrum when the symmetry breaking occurs at tree level, together with the white-dwarf white-dwarf binary foreground, and LISA and BBO sensitivity curves. Here we assume v w = 1.
For the strongest of PTs, we expect the wall velocity to be close to luminal, v w 1. This is because the mean field potential typically satisfies, This is the Bodeker-Moore (BM) criterion [40]. We shall make clear on our plots where the BM criterion holds and what we assume regarding v w when it does not. Even when the BM criterion holds, however, the wall is not expected to runaway with γ → ∞, due to the transition radiation effect from the gauge bosons [41]. Therefore, provided the gauge boson population has not been overly diluted by false vacuum inflation, energy released in the transition is transfered to the radiation bath, in which sound waves [39,42,43] and magnetohydrodynamic turbulence [44,45], rather than the bubble wall collisions directly [46][47][48][49][50][51], lead to a gravitational wave signal.
In summary, determining T n , α, and β from eqs. (3.1), (3.3), (3.6) and (3.7) allows us to find h 2 Ω GW (f ) by means of the spectra summarized in appendix B. An example of the spectrum is shown in figure 3. We use estimated sensitivity of the gravitational wave detectors to stochastic backgrounds, h 2 Ω sens (f ), for LISA [8], BBO [24], and when applicable the Einstein Telescope (ET) [52][53][54] (for which we use the updated sensitivity curve from [55]). The signal-to-noise ratio can be estimated using [8] where t obs is the time of observation in years. We assume t obs = 5 throughout.
Confusion noise from astrophysical foregrounds may be an issue at these frequencies. We shall compare to some estimates of the unresolvable components given in the literature.

JHEP05(2019)190
The ensemble of white dwarf -white dwarf (WD-WD) binaries are thought to be the dominant source of this foreground, exceeding the unresolvable neutron star -neutron star (NS-NS) foreground [56,57]. In this work we restrict ourselves to the foreground from the extragalactic WD-WD ensemble and use the central value given in [56]. We make this choice because, in contrast to the extragalactic ensemble, it is thought the WD-WD Galactic foreground [58][59][60] can be subtracted [61,62]. The continuous extragalactic NS-NS foreground extends to higher frequencies in the BBO band, however, it is thought that this can also be subtracted [63,64]. We also adopt an alternative, foreground-limited, estimate of the signal-to-noise ratio in which we attempt to naively capture the degradation of the sensitivity once the foreground, h 2 Ω FG (f ), is taken into account. The aim of introducing eq. (3.10) is to be able to roughly capture, in a single number, whether the signal extends above the sensitivity and foreground estimate. Whether such a PT signal could actually be separated from the astrophysical foreground depends, of course, on a myriad of factors, e.g. the robustness of the estimates of the amplitudes and spectral shapes of the signal and foreground, together with the confidence in our knowledge of the instrumental noise. These are topics worthy of further study, but we will not attempt to do them justice here. Nevertheless, we would like to remark that the LISA SNR FGL value associated to the spectrum of figure 3 clearly illustrates the importance of astrophysical foregrounds, even though they are often ignored in similar studies. Furthermore, we wish to emphasize that our sensitivity analysis in terms of SNR significantly improves that from ref. [20], where GWs from spin-one DM were briefly discussed.

Symmetry breaking at tree level
As discussed above, in this case the DM production proceeds via the standard freezeout mechanism. Interesting for us is the regime with a large g D , as this will lead to a strong phase transition. This pushes us to large m A and v η , see eq. (2.9), and the dark phase transition will generally occur prior to the EW one. Thus the task of studying the phase transition reduces to one dimension in field space. We have seen an example of the gravitational wave spectrum, together with the dominant foreground, in figure 3. We have also fixed θ to 0.1 and 0.01 (motivated by present and future direct detection constraints as shown in the left panel of figure 2), scanned over the parameters m A and m h D , 6 and calculated the GW signal. The result is shown in figure 4. Likewise, using the expressions of appendix B, we calculate the peak frequency and the peak GW energy density for each point of the parameter space and show the results in figure 5. These plots can be understood as follows. A larger DM mass, m A , requires a larger gauge coupling in order to return the observed DM density. This results in a stronger phase JHEP05(2019)190   transition from the one-loop effects of the gauge bosons. Similarly, in analogy with the SM, a lighter dark Higgs -corresponding to a smaller quartic λ 2 -also leads to a stronger transition because the broken phase minimum is shallower. Nevertheless, for particularly large values of m A /m h D , the one-loop effects can raise the broken phase minimum too far, resulting in the Universe becoming stuck in the symmetric phase. The latter can either be a false or true minimum, corresponding to the orange and red shaded regions of the figures respectively. We expect the allowed parameter space to be increased somewhat, into the orange region, if S 4 nucleation at lower temperatures were to be taken into account.

Radiatively-induced symmetry breaking: standard freeze-out and supercool DM
Another possibility is to impose classical scale invariance on the theory, as explained above. This scenario with our field content has been studied in [19,20]. Such a potential typically exhibits a large amount of supercooling [65][66][67][68][69][70][71][72][73][74][75][76][77][78][79][80][81][82][83]. This is because, lacking a mass term, the T = 0 potential is very flat in field space. Furthermore, the positive thermal corrections from the gauge bosons will lead to a barrier being present for any finite T . Somewhat counter intuitively, a smaller g D actually leads to more supercooling because the T = 0 potential becomes shallower, as shown in figure 6. The thermal barrier also becomes smaller, but the shallower potential ends up being the more important effect.
Of importance for the DM relic density in this scenario, is not just the DM annihilation cross section, but also the details of the phase transition. In particular the nucleation temperature, T n , the temperature when inflation starts, T infl , and the reheating temperature, T RH . The latter two quantities are calculated following the methods in [19].
Furthermore, due to the large amount of supercooling, the PT may actually not take place before the temperature falls to T ∼ Λ QCD . In this particular case, the SU(2) D PT is induced by QCD effects [65,67,79,82]. Our calculation of the nucleation temperature, ignoring the QCD trigger for now, is shown in figure 6.

JHEP05(2019)190
As a result, in the classical invariance scenario two distinct possibilities for the relic density can play out.
(ia). T n > Λ QCD . There is a large thermal abundance of massive gauge bosons after the phase transition, i.e. if T RH /m A and g D are large enough to bring (or keep) the gauge bosons in thermal equilibrium. Therefore, following the phase transition, the relic density is set through the usual freeze-out mechanism. Typically this occurs for gauge couplings g D ∼ 1 and m A 1.2 TeV.
(ib). T n < Λ QCD . This is similar to above, except the sequence of PTs is switched. Most of the parameter space corresponding to this regime has been ruled out by direct detection [19], except for the mass range 0.9 TeV m A 1.2 TeV, see figure 2.
(iia). T n > Λ QCD . There is sufficient supercooling for a period of late time inflation to take place. Before the phase transition, the gauge bosons are massless and have a large abundance. This abundance is diluted away by the period of late time inflation. The relic density in principle consists of the diluted, now super-cool, population of gauge bosons, together with an additional sub-thermal component created through scatterings after reheating. Numerically, however, we find the sub-thermal population is negligible in the parameter space corresponding to this regime, leaving the DM relic abundance set by the super-cool population of gauge bosons. The parameter space here corresponds to g D ∼ 1 and m A 370 TeV.
(iib). T n < Λ QCD . This is again conceptually similar to above except the PTs are switched. The sub-thermal DM population is now important for a large range of the parameter space, which corresponds to g D 1 and m A 370 TeV.
In all regimes, once the relic density constraint is used, we are left with one free parameter which we take to be m A . Here we wish to point out, supported by our calculations, that large portions of the parameter space of the classically scale invariant scenario can be probed through GW observatories. We shall now in turn discuss the GW signal in regimes (ia) and (iia), which both exhibit promising GW signals. Regimes (ib) and (iib) on the other hand, which are less promising and include larger uncertainties, are relegated to appendix C.
Regime (ia). The GW signal in this regime has previously been discussed in [20]. Here we provide our own -updated and expanded -calculation for completeness. For simplicity we assume the spectrum is given by the sum of the sound wave and turbulent contributions over the entire regime, although H becomes vacuum dominated in the lower m A range. For justification of this choice, together with details of the GW spectrum used, see appendix B. The key phase transition parameters are shown in figures 7 and 8, together with the foreground-free and foreground-limited signal-to-noise ratios. Note reheating is efficient in this regime: there is no period of matter domination immediately following the PT, as the decay rate of the inflaton is sufficiently large, Γ > H.  As can be seen from figure 7, LISA can probe DM masses in this regime up to m A ∼ 4 TeV, even in the presence of the WD-WD foreground. This is more than competitive with projections for future direct detection experiments [32,[35][36][37], which can probe up to m A ∼ 2 TeV [19]. (The current direct detection constraint demands m A 0.9 TeV [19,31,33,34].) The BBO proposal could test the entire parameter space shown here, well into what corresponds to the neutrino floor for direct detection experiments. Note for m A 1.2 TeV we find ourselves in regime (ib), which is discussed below.

JHEP05(2019)190
Regime (iia). Following the methods in [19], we find this regime corresponds to parameters g D ≈ 1 and m A 370 TeV. Notice that these DM masses are well above the usual unitarity constraint from the thermal freeze-out of DM [84,85], which does not take place here. Numerically the required g D grows slowly, from g D = 0.95 for m A = 370 TeV, to g D = 1.02 for m A = 10000 TeV. Our calculation of T RH and T infl is shown in figure 9. In this regime, reheating is inefficient following the PT, thus T RH = T infl . Indeed, there is a period of matter domination following the PT, as η oscillates about the minimum of its potential. More precisely, the ratio of scale factors between the PT and the end of JHEP05(2019)190 Figure 9. Left: the temperature when inflation starts, T infl , the reheating temperature, T RH , and the nucleation temperature, T n , in regime (iia). The ratio (T infl /T RH ) determines the amount of additional redshifting of the signal due to the matter dominated reheating period following the PT. Right: the detectability of the GW signal in regime (iia). Here the BM criterion holds over the entire range. reheating is given by This leads to greater expansion of the Universe between the PT and today, suppressing the signal, and redshifting the frequency further than would otherwise be the case. The GW spectrum is determined in the following way. First of all, due to the large amount of supercooling the scalar field configuration -and not sound waves or turbulence -is the source of the signal (see appendix B.2 for further discussion). It has been suggested that the oscillations of the scalar field after the PT may increase both that peak frequency and energy density of the GW signal by an order of magnitude [86]. We choose to remain conservative, however, and base our spectrum on the non-oscillating scalar field contribution, as indicated in appendix B.1. Once the Universe enters the late inflationary stage at T infl , the energy density remains constant until the plasma temperature reaches T n , and so the Hubble scales at both temperature are the same H(T n ) = H(T infl ). Taken

JHEP05(2019)190
together, β and H ∼ T 2 infl /M Pl set the initial frequency of the GW signal. We then redshift this value to T RH when the Universe once again enters a radiation dominated phase. The redshifting from T RH to today then follows the standard calculation [8]. Taking all this, together with eq. (3.11) into account, the peak frequency is given by where g * counts the effective degrees of freedom contributing to the radiation density, and f * /β = 0.62/(1.8−0.1v w +v 2 w ) is taken from simulations [51]. Due to cancellations between the various factors, we find the peak frequency here is ∼ 10 −2 Hz and almost independent of m A , as shown in figure 10. The amplitude of the spectrum is also suppressed with respect to the case with no early period of matter domination, because Ω GW = ρ GW /ρ c ∝ a 0 , (a −1 ) in a radiation (matter) dominated Universe. Accounting for these factors, we find the GW spectra and summarise their detectability in figure 9. Examples of the spectra are shown in figure 11. Notice that for these large masses, the frequency of the gravitational waves extends well above 1 Hz, motivating us to compare our signal against sensitivity curves from current and future LIGO configurations O1 and O5 [64], and ET [52][53][54][55]. Finally note we have explicitly checked the phase transition completes even though we are in the vacuum dominated regime. Details are presented in appendix D.

Discussion and conclusion
We have explored the possibility of spin-one DM from a hidden SU(2) D gauge group. The stability of DM is elegantly assured through a custodial symmetry. Given the massive vector bosons, unitarity demands that the SU(2) D be broken through the Higgs mechanism. This implies a phase transition or crossover occurred in the dark sector, i.e. the symmetry was initially unbroken at high temperatures. A strong phase transition will result in gravitational waves possibly detectable at future gravitational wave observatories. In this scenario the SU(2) D gauge coupling plays a crucial role in determining the relic abundance through freeze-out or late-time inflation. The same gauge coupling controls both the scattering cross-section and the thermal effects of the gauge bosons relevant for the phase transition. The model is therefore well suited as a case study for the sensitivity of future gravitational wave observatories to phase transitions in DM sectors.
We studied both tree level and radiatively-induced symmetry breaking. After finding the resulting gravitational wave spectra we identified parameter space which can be probed by LISA and BBO. As is known from previous studies, only limited parameter space of standard polynomial type potentials can be tested by LISA. The prospects improve for the classically scale invariant scenario. In this case, LISA is competitive with future direct detection experiments in the freeze-out regime and can probe the new regime of super-cool JHEP05(2019)190 Figure 11. Examples of GW spectra in regime (iia). Although α 1, and β/H is similar for both phase transitions, the period of matter domination after the PT is longer for larger m A , leading to a suppressed signal. For purposes of illustration we also include the unresolvable foreground from black hole binaries with masses 10 2 M − 10 10 M (MBH-MBH) [57]. DM, which is inaccessible to direct and indirect detection. Nevertheless, a conclusive test could only be performed by a more powerful observatory such as BBO.
We saw how foregrounds, which have so far been largely ignored in phase transition studies, apart from in [6,7], can be taken into account in the estimates of the signalto-noise ratio. Our results should be taken as indicative; we expect updated estimates of foregrounds to become available as our knowledge of the binary populations improves. More sophisticated studies, taking into account the precise capability of the LISA and eventually BBO spacecraft are required. Simulations of sound waves in the plasma for α > 0.1 should also be performed. Only then will it be possible to conclusively rule out models from their implied gravitational wave signals using future LISA and BBO data. A positive signal at LISA -which requires a very strong phase transition -would most likely point toward exotic new physics at the TeV scale such as the close-to-conformal potential studied here.
Note added. After our paper was released on arXiv, ref. [87] appeared, which confirmed the validity of our one-dimensional field space approximation for the classically scale invariant potential.

A.1 Symmetry breaking at tree level
The full effective potential is composed of four pieces The tree-level piece. This directly follows from eq. (2.1) and it is given by Coleman-Weinberg potential at zero temperature. Knowing the field dependent masses, m i (φ, η), in the Laudau gauge the one-loop T = 0 contribution is given by where µ is the MS renormalization scale and g i = {1, 3, 6, 12, 1, 9, 3, 3} for the h, Z, W ± , t, η, A, G, G D . In addition, C i = 5/2 for gauge bosons and C i = 3/2 otherwise. Finally, F = 0 (1) for bosons (fermions). The masses as a function of the scalar field values for the fermions and gauge bosons of the SM are Similarly, for the dark gauge bosons m 2 A (φ, η) = g 2 D η 2 /4. Due to the coupling λ 3 in eq. (2.1), the scalar sectors mix with each other and the masses for the real scalar fields entering in eq. (A.3) are the eigenvalues of the matrix In spite of this, the Goldstone bosons do not mix at tree level. In fact, in the Landau gauge, their masses are given by which vanish at (φ, η) = (v φ , v η ), as follows from eqs. (2.3).

JHEP05(2019)190
The counter-term potential. The counter terms to the potential in eq. (A.3) are By demanding no changes to the masses and VEVs of the scalars from their tree level values, that is, by imposing we calculate the couplings in eq. (A.8). Moreover, we find 7 In this equation, the prescription for the scalars m 2 i (φ, η) is the following. They are the eigenvalues of the mass matrix in eq. (A.5), ordered in such way that Notice that Σ ± F ± (a, b) = a + b and that, when a and b are both positive (negative), F + (a, b) is the maximum (minimum) of them.
Finite-temperature potential. The one-loop finite T contribution is given by We evaluate these integrals numerically. In order to take into account the resummation of the Matsubara zero modes one includes the daisy term There is a subtle issue for the contribution of the Goldstone bosons to eq. (A.3). As explained above, their tree-level masses vanish at (v φ , vη) leading to an infrared divergence in eq. (A.3). Such a divergence is spurious [88] and disappears after accounting for the one-loop contributions to the Goldstone-boson self energies. Neglecting any possible mixing effect due to non-vanishing λ3, the latter can be calculated by means of δm 2

JHEP05(2019)190
where the sum runs only over scalars and the longitudinal degrees of freedom of the vector bosons, i.eḡ i ≡ {1, 1, 1, 2, 1, 3, 3, 3} for h, Z, γ, W ± , η, A, G, G D . Here the thermal masses are given by [89] Π Higgs = where n f = 3 is the number of fermionic families with SU(2) × U(1) charge. Note for the scalars and the Z/γ, the prescription here is that m 2 i (φ, η) represents the relevant eigenvalue of the zero temperature mass matrix and m 2 i (φ, η)+Π i (T ) the relevant eigenvalue of the zero temperature mass matrix with the thermal masses added along the diagonal. This means the Z and γ mix at finite temperature. To avoid spurious contributions to the thermal masses from the SU(2) D gauge bosons at large field values, we cut off the g D contributions with a factor (m A /T ) 2 K 2 (m A /T )/2, where K 2 (x) is the modified Bessel function of the second kind of order two.

A.2 Classically scale invariant potential
As explained in the text, in this case we have radiative symmetry breaking and the potential at one loop becomes [19,20,27,28] V 0 1 (η) where the φ direction plays a completely negligible role in the area of parameter space in which we shall be interested. (The EW symmetry is broken by the induced mass term, If directly after the PT the Universe becomes radiation-dominated, the stochastic GW background receives a number of contributions, summarised in [8]. First, if no significant

JHEP05(2019)190
plasma is present, the scalar field contribution [51] dominates Alternatively, if a significant plasma is present, the following contributions dominate Note we do not sum eq. (B.1) with (B.2), following the updated recommendations in [41].
Here H * is the Hubble parameter when the GWs are emitted and κ η , κ v , and κ turb are the fractions of the latent heat that is converted into energy density in the scalar field, bulk motion of the fluid, and turbulence respectively. These can be calculated in terms of the wall velocity and α. For this, we use the expressions reported in refs. [8,90]. In addition, the spectral shapes are given by 3) The corresponding peak frequencies are where f * /β = 0.62/(1.8 − 0.1v w + v 2 w ) is taken from simulations [51]. If the after the PT there is a period of matter domination, the previous expressions must be rescaled as explained in the main text. See eqs. (3.12) and (3.13).
In this work, we are interested in phase transitions with significant supercooling, which lead to signals possibily observable by LISA. We caution that for strong phase transitions, α > 0.1, considerable uncertainty enters into the use of eqs. (B.1)-(B.4), despite the fact that they become insensitive to α when the latter takes values much greater than one. This is partly due to the uncertainties in determining the relevant contribution, which have only been started to be explored following the updated results of [41] (see appendix B.2 below). Furthermore, numerical simulations of sound waves in the plasma have also only been performed up to α = 0.1. Until further simulations have been performed, we are therefore resigned to extrapolating the results to large α, as in [8].

B.2 Determination of the relevant contribution
We provide an estimate to justify our use of the sound wave plus magnetohydrodynamic turbulence in regime (ia), and the scalar field contribution in regime (iia), of the classically scale invariant potential. This is crucial as in the vacuum dominated regime, depending on the amount of supercooling, the energy transferred to the plasma and hence sound waves and magnetohydrodynamics may become negligible. Bodeker and Moore [41] find a nextto-leading-order contribution to the pressure difference due to transition radiation across the wall which, using the parameters relevant to our model, scales as

C Super-cool DM regimes (ib) and (iib)
In this appendix we comment briefly on these regimes, in which the SU(2) D phase transition occurs after QCD confinement, although they are less promising from the point of view of GWs. In these regimes the QCD phase transition occurs with six massless quarks and is first order [92,93]. There is a chance this could lead to an observable GW signal [9,79] (although note the vacuum domination continues for some time after the QCD PT, diluting the signal). We cannot, however, use our techniques above to accurately calculate the phase transition parameters, α and β for the QCD phase transition [79]. Most likely a lattice study is required in order to more carefully explore this possibility. Alternative techniques have been pursued in [77,80].
We now turn to the details of the SU(2) D phase transition following after the QCD one. The quark condensate formed after chiral symmetry breaking leads to a tadpole term and hence a VEV for the EW Higgs. This in turn leads to a mass term for η through the cross-quartic. Provided 3m 2 A 2m 2 h , which corresponds to our regime of interest, the thermal barrier from the gauge bosons is still large enough to prevent immediate SU(2) D breaking. Instead, a first order phase transition occurs just before the barrier disappears at T ∼ m h Λ QCD /m A . As can be checked numerically, the non-zero mass term for η means this phase transition now occurs with a very large β, and does not lead to an observable gravitational wave signal. An example calculation showing these conditions are met is shown in figure 12. The underlying reason this can occur is because the nucleation rate continues to grow exponentially as the temperature falls, allowing sufficient bubble formation to overcome the Hubble driven expansion of the false vacuum. In contrast, for very strong phase transtions in standard polynomial type potentials, S 3 /T may remain constant or even grow as the temperature drops. Meaning the phase transition may not complete even if the nucleation condition Γ ∼ H 4 is achieved [91,[94][95][96][97].

D Completion of the phase transition
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.