LIGO as a probe of Dark Sectors

We show how current LIGO data is able to probe interesting theories beyond the Standard Model, particularly Dark Sectors where a Dark Higgs triggers symmetry breaking via a first-order phase transition. We use publicly available LIGO O2 data to illustrate how these sectors, even if disconnected from the Standard Model, can be probed by Gravitational Wave detectors. We link the LIGO measurements with the model content and mass scale of the Dark Sector, finding that current O2 data is testing a broad set of scenarios where the breaking of $SU(N)$ theories with $N_f$ fermions is triggered by a Dark Higgs at scales $\Lambda \simeq 10^8 - 10^9$ GeV with reasonable parameters for the scalar potential.


I. INTRODUCTION
Much of the Universe is dark, and many theories have been built trying to explain it. Our hopes for probing these theories often rely on their possible connection to regular matter via some form of non-gravitational interaction For example, direct searches for Dark Matter hinges on some sort of coupling to nucleons or electrons, and constraints on those couplings usually assume a mechanism of communication between the Dark Sector and the rest of the Universe which establishes some form of tracking between these two sectors.
The first observation of Gravitational Waves (GW) by the LIGO and Virgo collaborations [1][2][3] in 2015 initiated a new way to see the Universe, and since then exciting new observations have provided information about astrophysical objects like Black Holes [4,5]. However, the physics reach for LIGO is not circumscribed to detection of mergers [6][7][8][9]. The detection, or the lack of, a stochastic GW background allows us to explore interesting, nonstandard sectors. We will explain how, with the current public data from LIGO, one can probe plausible Dark Sector scenarios, regardless of their non-gravitational interaction with visible matter.
These Dark Sectors could resemble Standard-Model dynamics, with new forces, Dark Higgses and states charged under them. Influenced by thermal contributions from the degrees of freedom in the Dark Sector, the thermal history of the Dark Higgs could then lead to firstorder phase transitions. Many studies have been devoted to the prospects that future interferometers could offer to explore Dark Sectors, e.g., Ref [10]. In this paper we explore the possibilities that LIGO and its current public dataset present, and bridge the gap between generic studies of thermal parameters, e.g., Ref [11], and specific particle-physics models.
The paper is structured as follows. In Sec. II we first describe the analysis of GWs from first-order phase transitions, then discuss in Sec. III the connection between the phase transition thermal parameters with classes of particle-physics models, especially of SU (N )/SU (N − 1) Dark Sectors. We finally link these models to current exclusions set by LIGO, and in Sec. IV we conclude the discussion.

II. GRAVITATIONAL WAVES FROM PHASE TRANSITIONS
The stochastic gravitational wave background (SGWB) is often considered as an isotropic, unpolarized, stationary and Gaussian background generated by a large number of unresolved gravitational-wave sources [6][7][8][9]. Its power spectrum is characterized by the dimensionless quantity where ρ GW is the energy density of the stochastic gravitational wave background, f is the frequency of the GW  and is the critical energy density of the universe today. In principle, the total SGWB is a superposition of all possible astrophysical and cosmological sources. However, we can obtain a conservative upper limit for the SGWB from phase transitions that occur in the early universe by assuming that phase transition dynamics is the main source of SGWB.
The SGWB generated from phase transitions in the early universe consists of three parts [12][13][14][15][16][17]: in which the three terms on the right hand side correspond to the contribution from bubble collisions, sound waves in the fluid and the turbulence, respectively. For simplicity, we shall assume in this work that contributions from sound waves are always dominant. We emphasize that this is typically the case for models in which gauge bosons acquire masses during the phase transition [18]. However, the analysis we present can be easily generalized to cases in which other types of contribution become more important. The phase transition is in general characterized by just a few parameters: the velocity of the bubble wall v w , the ratio of the free energy density difference between the true and false vacuum and the total energy density, ξ, the speed of the phase transition β/H, and the nucleation temperature T N . With these parameters, the GW power spectrum can be expressed as [19] Ω GW h 2 8.5 × 10 −6 g * 100 where g * is the effective number of relativistic degrees of freedom at the time of the transition, Γ ∼ 4/3 is the adiabatic index,Ū 4 f ∼ (3/4)κ f ξ is the root-mean-square fluid velocity with the efficiency parameter given by the approximate expressions [15] The spectral shape S sw is given by with the peak frequency The shape of S sw is shown in FIG. 1, noting that S sw is equal to Ω GW (f )/Ω GW (f pk ). In this figure one observes that varying the peak frequency f sw amounts to simply shifting the spectrum horizontally. Also note that the asymptotic behavior of S sw goes as ∼ f 3 for f f sw (dotted line), whereas one expects a behaviour ∼ f −4 for f f sw (dashed line). The behavior of S sw in the LIGO frequency range (indicated by the two vertical lines) can therefore transition from a simple descending power law, to one with a peak in between, and eventually to a ascending power law as we increase f pk .
With this spectrum, we follow the procedure laid out in Refs [20,21] to compute the upper limit of Ω GW for different values of f pk using data from LIGO O2 [21].
More details can be found in Appendix A. The only difference respect to these references is related to the choice of the optimized estimator, i.e., the estimatorΩ f ref is optimized by putting Ω GW (f )/Ω GW (f ref ) instead of some power of f /f ref . Note that this means the upper limit essentially depends only on the shape of S sw within the LIGO frequency band, since all the other factors drop out when taking the ratio. The 95% confidence level upper limit is then obtained by setting In FIG. 2 we show the upper limit at f ref = 25 Hz by the curve with triangles. For f pk 10 Hz or f pk 10 3 Hz, the bound approaches a constant since the spectrum is essentially a power law with fixed exponent within the LIGO frequency band. This behaviour allows us to easily extrapolate the constraint to even smaller or bigger values of f pk . However, for intermediate f pk there is a smooth transition on the upper limit between the ascending and descending asymptotic behavior. We see that LIGO provides a stronger constraint on ascending spectra than on the descending spectra, and the difference can be as large as an order of magnitude.
In practice, LIGO could only provide reliable constraints on SGWB in the band 20 − 1726 Hz and we have chosen f ref = 25 Hz as it is the frequency where LIGO is most sensitive [20,21]. However, the thermal parameters of the phase transition are connected to the amplitude at the f pk . In order to place constraint on the thermal parameters, we notice that, for a given GW spectrum as in Eq. (6), the constraint at a particular frequency in the LIGO band can be mapped to the peak frequency f pk using Ω up. lim.
The curve made by circles in FIG. 2 shows the result of this mapping. For f pk very small or very large, the upper bound from LIGO becomes substantially weaker, even when comparing with the constraint from BBN, which is obtained by integrating out the spectrum and requiring h 2 d ln f Ω GW (f ) < 1 × 10 −6 [9,22].
In what follows, we shall discuss the LIGO constraints on the thermal parameters and how such constraint can be utilized to constrain particular models of phase transition.

III. SCENARIOS FOR PHASE TRANSITIONS AND THEIR THERMAL PARAMETERS
Typically, one represents phase transitions as driven by the dynamics of a scalar field which transitions from one vacuum to another under the influence of the evolving thermal potential. In that context, the thermal parameters β/H and ξ are defined by in which S E is the Euclidean action, V is the thermal potential of the scalar field, and the bubble nucleation temperature can be obtained by solving in which ρ N = g (T N )π 2 T 4 N /30. On the other hand, the calculation of the bubble-wall velocity v w for a particular model is highly non-trivial. Therefore, instead of directly calculating it, we shall follow the customary convention of considering a few reference values, v w = 0.5 and v w = 1.
To connect the thermal parameters, ξ and β/H to specific models we will follow the approach described in Ref. [10]. We will consider classes of potentials which consists of competing terms with alternating signs. Specifically, we will look into two types of finite-temperature potentials in which the coefficients of the scalar field h D are all positive at the time of transition. Indeed, most dark phase transitions can be mapped onto these effective scenarios. For example, in a Dark Sector where its particles acquire mass from a Dark Higgs as its gauge group SU (N ) breaks into SU (N − 1), the potential in Eq. (13) can be realized with renormalisable operators, whereas the potential in Eq. (14) can be obtained from non-renormalisable sextet interaction [10].
In what follows, we shall discuss the Dark Higgs -SU (N )/SU (N − 1) models with renormalisable and nonrenormalisable operators specifically.

Models with renormalisable operators
For the type of potential in Eq. (13), we can parametrise zero-temperature parameters as in which v D is the zero temperature vacuum expectation value, and Λ is the scale of the potential. With this parametrisation, the finite temperature potential can be expressed as In the second equality, the following high temperature expansions are used and in which the field dependent masses can be read as Note that in the second line of Eq. (16), only the massive gauge bosons are taken into account, i.e., the higher order term (∼ m 3 /T 3 ) from the Goldstone boson is neglected.Besides, the part of the expansion which gives rise to terms independent of h D is neglected, since it only amounts to a constant shift in the potential V . Finally, the mapping to the temperature dependent parameters in Eq. (13) is straightforward by matching the terms with the same powers of h D .
Eq. (16) is also often written in the following form in which the minimum of the potential can be easily obtained by minimizing the potential By identifying one finds With these, for α(T ) ∈ [0.51, 0.65], the effective action can be fitted by [23] S with the fitting parameters a = −71.06, b = 71.62, c = 0.008805, and d = 0.009263.

Models with non-renormalisable operators
For the type of potential in Eq. (14), similar to the previous case, one can perform the following parametrisation The finite temperature potential then becomes Note that, in the high-temperature expansion, the cubic term is assumed to be subdominant, i.e., we have only kept the part proportional to m 2 /T 2 [10]. Moreover, the field-dependent masses of the Goldstone bosons and the Dark Higgs which goes into the thermal correction are The terms proportional to αh 4 D /v 4 D would give rise to the thermal correction of the quartic term.
Just as we have done in the previous section, one can write Eq. (29) in terms of temperature-dependent parameters: The non-vanishing VEV v D (T ) = λ(T ) + Λ 2 (T ) − 24c 6 (T )m 2 (T ) 12c 6 (T ) is obtained by minimizing the potential. Suppose the non-vanishing VEV does exists (Λ 2 (T )−24c 6 (T )m 2 (T ) > 0), then α(T ) can be obtained by solving Therefore, Following [10], the Euclidean action can be fitted by Finally, for g (T N ), at high temperatures, i.e., T O(100) GeV m i , m h D , we shall assume that the particles in the Dark Sector are the only degrees of freedom in addition to the SM. Therefore, the effective number of relativistic degrees of freedom will be given by where we have included all degrees of freedom from the SM, as well as the massless and massive gauge bosons and the dark fermions charged under SU (N ).

Results
In FIG. 3, we present the LIGO constraints on the thermal parameters at 95% confidence level together with the constraint from BBN. In the left panels, the constraints are evaluated with v w = 0.5, whereas the bounds in the right panels are evaluated at v w = 1.
Note that, with v w fixed, the thermal parameters (β/H, T N , ξ) constitutes a 3-dimensional space. In our analysis, we find the the excluded region can be conveniently projected on the 2-dimensional plane of β/H and T N as the contours of ξ do not intersect each other. This thus enables us to use colour variation to represent the values of ξ on the exclusion contours. Indeed, within a particular contour, any value of ξ larger than the value associated with the contour is excluded. As a result, the regions in the 2D plots where the LIGO ellipsoidal contours have a colour lighter than the BBN vertical contours suggest that LIGO can set a better bound than BBN, and viceversa.
On the same figure, we also plot the values of the thermal parameters obtained from different models defined by Λ, N and N FL . We show examples with N = 2 and 5, and a fixed number of flavours N FL = 3, as representative of the behaviour expected from the choices on matter content and gauge groups. We also fix g = 1 and y = 1, although similar behaviour is found for similar O(1) values.
If a set of thermal parameters obtained from a particular model lies within a contour, and has a larger ξ (darker colour) than the corresponding contour, then this set of thermal parameter is excluded at 95% CL.
From the top figures 3, one can deduce that the constraint from LIGO is still unable to probe the renormalisable classes of models, as the markers representing the thermal parameters barely touch the excluded region.
However, in the non-renormalisable case, we indeed see how LIGO is able to probe these Dark Sectors. In particular we observe that the dark scales Λ ∼ 10 8 -10 9 GeV are best constrained by LIGO. Note that models with larger group rank N and number of flavours N FL are more constrained as they tend to produce a larger ξ.
The regions excluded by LIGO O2 correspond to Λ in the region around 10 8 -10 9 GeV. In this region we have varied v w from 0.5 and 1 and explored different options for N and N FL as illustrated by FIG. 3. The range of parameters α and x ≡ v D /Λ where LIGO currently has sensitivity is given by These regions can be translated into parameters in the scalar potential Eq. (14) by inspecting Eq. (28): Clearly, these are reasonable choices of parameter space and indicate that LIGO is testing interesting Dark-Sector theories.

IV. CONCLUSIONS
When discussing Dark Sectors and their GW signatures, we usually think on future probes like LISA, many years from now. Here we have shown that LIGO is already probing interesting scenarios for Dark Sectors.
In the context of first-order phase transitions and LIGO data, the emphasis has been placed in performing effective analyses, such as Ref. [11] where the authors explore the bounds from LIGO O3 using parametrisations of the power spectrum. In this paper, we take a complementary step and focus on examining whether concrete particle-physics models could be related to the tested regions.
In this work we answer the question whether the type of first-order phase transitions LIGO is currently probing could be represented by concrete, reasonable particlephysics models. For that reason, we have focused on classes of models which capture a broad set of features of Dark Sectors, and, at the same time, capable of producing interesting GW signatures. The renormalisable and non-renormalisable benchmarks we used had been identified in Ref. [10] as more promising for strong first-order phase transitions from Dark Sectors.
We choose to set up those models with a breaking SU (N )/SU (N − 1) which should be understood as an example in which some bosonic and fermionic degrees of freedom influence the thermal history of the Dark Higgs. Such a choice also allows a simple parametrisation in terms of the group rank N and the number of flavours N FL . We find that scales around 10 8 − 10 9 GeV are better probed by LIGO O2. We also find that the sensitive regions correspond to moderate values for N and N FL , evidencing that LIGO is not testing extreme regions in the UV parameter space. Of course, various other types of models for phase transition could also generate GW whose spectrum lies in the frequency range relevant for LIGO, e.g., models motivated by grand unification theories [24,25] and models for confinement-deconfinement phase transition [26]. It is straightforward to see whether LIGO might be able to constrain those models once the thermal parameters are computed.
Note that, when obtaining the GW peak amplitude in Eq. (4), we have used the standard formula from Ref. [19]. Recent discussions in Ref. [27] suggest the existence of an additional suppression factor due to the finite lifetime of the sound waves. Moreover, it has been shown recently that theoretical uncertainties in GW production could lead to changes in the GW spectrum as large as several orders of magnitude [28]. In addition, we have assumed in our analysis that sound waves are the dominant source for GW production. In other scenarios, for example, strongly supercooled phase transitions [29][30][31][32], or scenarios in which new heavy particles can provide sufficiently large friction [33], contribution from turbulence or bubble collision could also be important and one needs to care about their effects in the GW power spectrum. Although subject to those uncertainties, our analysis nevertheless continues to provide a concrete method to constrain Dark-Sector models with LIGO data.
Our analysis of the SGWB is based on publicly available O2 data. In Appendix A, we have provided ex-planations on how to reproduce our analysis. Recent papers from authors in the LIGO/Virgo collaboration, e.g., [11,34], make use of the O3 data. GW constraints on particle-physics models considered in this paper are expected to improve as more data becomes available.
We believe our results motivate a more systematic study of the particle-physics scenarios that the LIGO experiment is able to test. We emphasize again that traditional direct or indirect searches for dark particles assume that the Dark Sector interacts non-gravitationally with the Standard Model. On the other hand, since gravity is universal, methods for probing the Dark-Sector via its gravitational effects such as structure formation (see Ref. [35,36] and references in it for recent progress) and gravitational waves do not rely on those assumptions. These gravitational effects therefore offer unique opportunities to access Dark Sectors, which would otherwise be hidden from us if they lack a connection to the Standard Model.
-the strain time series h(t i ) and some auxiliary data such as the data quality (DQ) mask label associated with each strain measurement. The DQ mask is a 7-bit binary number each of which indicates whether a certain type of check is passed (value=1) or not (value=0). We convert this binary number into a decimal digit. For example, there is no data at time t i if DQ(t i ) = (0000000) 10 = 0. On the contrary, data is present if this value is nonzero.
Following the LIGO stochastic gravitational wave analysis [20,21], we first downsample the 16 kHz strain data to 4 kHz. We then select out the timestamps at which both detectors are taking data properly, i.e., the times t i at which both DQ H (t i ) = 0 and DQ L (t i ) = 0. After doing that we get, for each data file, a list of GPS times at which data in both detectors is available. We then combine the lists of GPS times within a file and across neighbouring files into continuous segments. The segments whose duration ≥ 600 s are further picked out to perform a stationarity cut following Ref. [20,37]. When we perform the stationarity cut, we also notch out frequencies at which the data exhibits narrowband coherent lines that are known to be instrumental or environmental artifacts [38][39][40][41]. The list of the notched frequency bands we used can be found on the public data release page: https://dcc.ligo.org/LIGO-T1900058/public.

Analysis
After selecting a clean list of strain data and GPS time which satisfies the requirements on data quality and stationarity, we then use the cross-correlation method to estimate the SGWB signal. The spectrum of the GW background is estimated with the cross-correlation statis-ticĈ(f ) [21], defined aŝ where S 0 = 3H 2 0 /(10π 2 f 3 ), H 0 is the Hubble parameter, s 1,2 (f ) are the Fourier transfroms of the strain data of both detectors, T = 192 s is the segment duration of the Fourier transforms, and ... indicates the average over all such 192s-segments. In the limit that the GW signal is negligible comparing to the instrumental noise, the variance ofĈ(f ) is given by where ∆f = 1/32 Hz, andP 1,2 (f ) are the one-sided power spectrum of each detector, which are obtained as an average over two neighbouring segments 1 ,