Probing a Dark Sector with Collider Physics, Direct Detection, and Gravitational Waves

We assess the complementarity between colliders, direct detection searches, and gravitational wave interferometry in probing a scenario of dark matter in the early universe. The model under consideration contains a $B-L$ gauge symmetry and a vector-like fermion which acts as the dark matter candidate. The fermion induces significant a large dark matter-nucleon scattering rate, and the $Z^\prime$ field produces clear dilepton events at colliders. Thus, direct detection experiments and colliders severely constrain the parameter space in which the correct relic density is found in agreement with the data. Nevertheless, little is known about the new scalar responsible for breaking the $B-L$ symmetry. If this breaking occurs via a first-order phase transition at a TeV scale, it could lead to gravitational waves in the mHz frequency range detectable by LISA, DECIGO, and BBO instruments. The spectrum is highly sensitive to properties of the scalar sector and gauge coupling. We show that a possible GW detection, together with information from colliders and direct detection experiments, can simultaneously pinpoint the scalar self-coupling, and narrow down the dark matter mass where a thermal relic is viable.


I. INTRODUCTION
The evidence for the presence of dark matter (DM) in the universe is compelling, and thermal relics stand out among the possible explanations for it.They typically experience interactions that could be probed by current and near-future experiments and easily yield the correct relic density, thus attracting much attention from the community.The nature of these particles allowed us to explore the so-called dark matter complementarity, which refers to the use of data from various sources, such as direct and indirect detection experiments, as well as colliders, as a way to narrow down the viable properties of a dark matter particle in a given model (see [1] for a review).With the detection of gravitational waves (GWs) by LIGO/Virgo/KAGRA [2][3][4][5], together with the NANOGrav evidence for a stochastic GW background [6], and the near-future launch of new-generation space-based interferometers [7], a new era has surfaced.If the scalar in the dark sector was responsible for a first-order phase transition, there could be a resulting a stochastic gravitational wave spectrum which could indirectly constrain the dark sector.For a DM model at TeV scale, the spectrum would be in the range of optimal detectability at LISA/DECIGO/BBO [8,9].In other words, gravitational waves have become an interesting laboratory for dark sectors that feature a scalar particle inducing a first-order phase transition.
In this paper, we illustrate the complementarity between colliders, direct detection experiments, and gravitational wave detectors for testing models with a dark sector.Our main argument is that any near future collider will likely be insensitive to the self-couplings of a dark scalar, but this parameter could be constrained from possible detections of cosmological GWs induced from this scalar sector in a first-order phase transition.On the other hand, the GW spectrum depends also on the gauge and fermionic sectors, so GW detectors alone would not be able to disentangle the information on the selfcoupling.More concretely, we will show, in a specific model, that none of these experiments alone could probe all dimensions of the parameter space, but a full picture could be obtained by considering the symbiosis of colliders, detectors, and GW experiments.
To be concrete, we work on a model with an additional U (1) B−L gauge boson Z ′ , together with a DM candidate, which is a Dirac vector-like fermion coupled to the visible sector via a Z ′ portal.By computing the spin-independent nucleon-DM scattering cross-section, we can impose constraints on the model due to direct detection limits from XENONnT and LUX-ZEPLIN (LZ) [10,11].Moreover, collider searches by ATLAS and LEP II impose bounds on the mass and coupling of the Z ′ [12][13][14].But these experiments can say little about the scalar sector responsible for breaking the U (1) B−L symmetry: the scalar is not relevant for the DM phenomenology and may decay only into invisible particles, making it difficult to detect at colliders.However, the U (1) B−L symmetry-breaking process may be a first-order phase transition in the early Universe, resulting in GWs testable at the aforementioned detectors, as shown in Ref. [15,16].The GW spectrum is highly sensitive to the scalar sector but depends also on the Z ′ and dark matter couplings.Thus, the parameters can only be disentangled with the aid of direct detection and collider experiments.We show that, once we have a measurement of these couplings, a de-tection of cosmological GWs could even lead to a measurement of the scalar self-coupling at good precision using GWs.This illustrates how GW detectors, collider, and direct detection searches are complementary probes for dark sectors and beyond the Standard Model (SM) physics in general.
The DM phenomenology has been explored in the context of the minimal B − L model.In [17], the authors contrasted their theoretical results with the bounds from collider and direct detection on a minimal B − L model in which the lightest right-handed neutrino is the DM candidate.In [18], they look into the DM phenomenology of a scalar singlet dark matter charged under B − L, where gauge and scalar interactions are considered.Moreover, for the latter case, semi-annihilation processes are very important regarding the relic density.
The discussion of GWs in the context of DM models is not new.Ref. [19,20] have explored GWs working in a secluded DM model and in a gauged U (1) ℓ leptonic Abelian symmetry extension, respectively.In [21], the authors studied a SU (2) 0 × SU (2) 1 × SU (2) 2 gauge extension.Ref. [22,23] investigated the DM phenomenology and GW production in two-component DM models.In [24,25] the authors worked in a similar vein for scalar sector extension models.
One could wonder that our paper is a combination of Ref. [17] and Ref. [15], however, our work differs from others because our dark matter candidate is taken to be a vectorlike fermion, and new collider bounds are derived and updated limits from direct detection experiments are employed.Moreover, we generated the GW spectra for choices of the gauge coupling and the scalar singlet self-coupling, bearing in mind the parameter space favored by the dark matter phenomenology.Thus, exploiting the interplay between dark matter and GW physics.Our emphasis is not on showing that the model could yield a detectable GW spectrum, but especially on highlighting the complementarity between GW detectors, collider searches, and direct detection experiments for probing the phenomenology of a model with a gauged B − L symmetry, opening also a portal-like to the scalar sector.This illustrates that the inclusion of GW detectors in the particle physicist's experimental arsenal will not render collider and direct detection experiments obsolete, but will only enrich their results.
This work is organized as follows: In Section II, we present the minimal B − L model describing the particle content and the properties of the scalar singlet and DM candidate that we are interested in.The generation of GWs via first-order phase transition in the context of the minimal B − L model is discussed in Section III, and the DM production in Section IV.We explain the DM model constraints in Section V. Finally, we discuss the results and present the conclusions in Sections VI and VIII, respectively.

II. THE MINIMAL B-L MODEL
The minimal B −L model is a simple extension of the Standard Model (SM).An additional U (1) B−L Abelian symmetry enlarges the SM gauge structure to SU (3 where B stands for baryon number and L for lepton number.Hence, the gauge content increases by one Table I.The matter particle content of the minimal B-L model and their respective charge assignments under each symmetry. new gauge boson, say, Z ′ .In the fermion sector, we add three Majorana right-handed neutrinos N iR (i = 1, 2, 3) to cancel out the gauge anomalies, and the DM candidate is a vectorlike Dirac fermion χ, which also carries B − L charge as the other fermions.Consequently, the DM-SM interactions happen through a Z ′ portal.Moreover, we include a scalar singlet Φ s , which breaks the B − L symmetry and gives rise to a Majorana mass term for right-handed neutrinos that is essential to realize the seesaw mechanism type I to turn active neutrinos to massive particles.In Table I, we present the matter particle content and their respective charge assignments under each symmetry, including the Z 2 , which arises after B − L breaking and ensures DM stability [26].
The Lagrangian that describes the DM phenomenology and encodes the process of B − L symmetry breaking can be written as where F ′ µν and g BL are the strength tensor and coupling of the B − L symmetry.The n j stands for the B − L quantum number of the particles j = χ, ℓ, ν ℓ , and q, with ℓ = e, µ, τ , and q = u, d, c, s, t, and d.We have neglected the kinetic mixing between the photon and the Z ′ boson in our study.This is well justified because the gauge coupling g BL in our study will be sizeable, whereas the bounds on the kinetic mixing impose it to be suppressed, say, as strong as 10 −2 , yielding no effect on our phenomenology [1].This limit was essentially due to the previous generation of 1-Ton Direct Detection facilities as XENON1T and, consequently, it is expected to be even stronger nowadays [27].Thus, the kinetic mixing can be safely neglected throughout.
We remark that the DM charge, n χ , must be other than ±1 to avoid DM decay via an additional Yukawa term involving χ R .Notice that H = iσ 2 H is the isospin transformation of the SM Higgs doublet H = ϕ + , ϕ 0 T .Moreover, in the scalar sector, we have the kinetic term of Φ s which will give mass to Z ′ , and the scalar singlet potential at tree level We leave further details about the effective potential V eff (Φ s ) that leads to the firstorder phase transition for the next section.The parameterization of the scalar field singlet is given by, where v s = √ 2⟨Φ s ⟩ is the U (1) B−L vacuum expectation value (VEV), and ρ s is the Goldstone boson which will be eaten by the Z ′ field after spontaneous symmetry breaking.Its mass arises from the kinetic term of scalar singlet, and the right-handed neutrinos N iR via the Majorana mass term in Eq. ( 1), We highlight that it happens at high-energy scales, say, v s ≫ v, where v is the VEV of the SM Higgs field, H.In the same way, the tree-level mass of the scalar singlet is given by The active neutrinos become massive via the popular type I seesaw mechanism, which nicely reproduces the neutrino data [28][29][30].
In Fig. 1, we display the relevant Feynman diagrams for DM phenomenology.The first and the second diagrams give the DM relic abundance.In the regime m χ < m Z ′ , the schannel overcomes the t-channel, and the cross-section is, where n f stands for the B − L charge of the SM fermion, and n f c is the number of colors of the final state SM fermion.However, when m χ > m Z ′ , the second diagram also contributes to the DM relic abundance whose annihilation crosssection is given by, In summary, the free parameters that govern the DM phenomenology are the DM mass m χ , the Z ′ mass m Z ′ , and the B − L coupling g BL .Concerning the third diagram in Fig. 1, it is relevant for direct detection, as we will see later.Having reviewed the key ingredients for the relic density, in the next section, we assess the GWs production in our model due to a first-order phase transition, where g BL and v s play a crucial role.

III. STOCHASTIC GW SPECTRUM FROM FIRST-ORDER PHASE TRANSITION
If the U (1) B−L breaking is a first-order phase transition, there are regions of broken phase in the plasma where the symmetry remains unbroken.These so-called bubbles expand and induce plasmatic motion in the form of sound waves and turbulence.At the end of the transition, when the bubbles collide and fill the entire space, GWs are produced due to a timevarying quadrupole moment in the kinetic energy-momentum of the plasma [31,32].
To estimate the shape of this spectrum, we need to study the dynamics of the phase transition.In the presence of a thermal plasma, the dynamics of the scalar field ϕ s is described by an effective potential which, at 1-loop order, takes the form The second line corresponds to the Coleman-Weinberg potential plus counter-terms added to ensure that the minimum of the 1-loop zero-temperature potential remains at v s , and that the mass of the scalar is still given by Eq. ( 5).We take into account the Z ′ running in the loop, and neglect the right-handed neutrinos (assuming small Yukawas1 ) and the scalar (since typically λ s ≪ g BL and it has only one degree of freedom, so its effect is subdominant against the Z ′ ).The thermal part V th is computed at 1-loop via standard methods of thermal field theory [33,34], and the last term accounts for the resummation of daisy diagrams, with the thermal mass for the Z ′ given by Notice that the value of the zero-temperature potential at the unbroken state ϕ s = 0 is parametrized by λ s , whereas at the broken vacuum ϕ s = v s it goes with g BL .Hence, for too small values of λ s the symmetric minimum may become the lowest energy state and U (1) B−L would remain unbroken, which is nonphysical and should be avoided.Explicitly, defining the energy difference between the two vacua ∆V (T, ϕ) ≡ V eff (T, ϕ) − V eff (T, 0), one finds that at zerotemperature implying λ s > 3g 4 BL /π2 .The bubble nucleation rate per unit volume is Γ nuc /V ∼ T 4 e −S3/T , where is the Euclidean action of the critical bubble, corresponding to the saddle point configuration that connects the two vacua in field space [34,35].At some temperature T , this configuration can be found by solving the bounce equation −∇2 ϕ + dV eff /dϕ = 0 using a shooting method.The nucleation temperature is found by imposing that Γ nuc /V equals the Hubble expansion rate.This corresponds roughly to the temperature at which S 3 /T ≈ 140 [32].We can then define the fractional amount of energy released by the transition [36] to be, where Q lat = T d∆V /dT − ∆V is the latent heat of the phase transition, the numerator Q lat − 3∆V is the difference of the trace of the energy-momentum tensor in both phases [36,37], ρ rad = π 2 g eff T 4 /30 is the radiation energy density at the time of the transition and g eff = 117.5 is the effective number of relativistic degrees of freedom in the plasma.The amplitude of the GW spectrum will depend on how efficiently this energy is converted into the fluid motion of the plasma.For that, we construct efficiency factors κ sw and κ turb such that the fractional energy in sound waves (respectively in plasmatic turbulence) is proportional to κ sw α (resp.κ turb α).An approximate formula for κ sw can be found in ref. [36], but for turbulence this conversion factor is unknown.Sometimes it is estimated to be 1 − 10% of κ sw [32,38], and we have explicitly checked that for these values its contribution is subdominant.Due to this indeterminacy and subdominance, we neglect here the contribution from turbulence and consider only sound waves.Moreover, we take into account the suppression factor of the sound wave spectrum due to the finite lifetime of this source and a possibly early onset of turbulence [39].Altogether, this means that our approach underestimates the spectrum, so our results are conservative.
Another important parameter for estimating the GW spectrum is the (inverse) duration of the transition, estimated as 2 Given that the bubble expands at a velocity v w , its radius at collision will be proportional to v w (H/β), and the larger the bubble, the larger the GW amplitude.Calculating the wall velocity v w is a daunting task since it involves non-equilibrium phenomena and depends on how we model the plasma away from equilibrium.There have been recent discussions in the literature on the appropriate way to achieve this description, but the debate is still unsettled [40][41][42].Here we approximate v w = 1 for simplicity, which should suffice for an adequate estimate of the spectra at the correct order of magnitude.More specifically, for the amplitude of the GW spectrum from sound waves we find [38,43,44], where is the suppression factor due to the finite lifetime τ sw = R/U of the sound wave source, where R = (8π) 1/3 v w /β is the typical bubble separation and U = 3 the plasma root mean squared velocity [39,43,44].
This spectrum has a peak at We will see that, for typical benchmark values of the parameters in our model, the peak frequency lies in the mHz band, hence possibly within reach of future interferometers such as LISA, DECIGO and BBO.We will now discuss the dark matter relic density and scattering rate to later put our findings into perspective.

IV. DARK MATTER RELIC ABUNDANCE
In this section, we assess the production of DM particles in the standard thermal freeze-out paradigm.In such a narrative, after reheating, the DM particles were in thermal equilibrium with the SM particles, which means that they were pairannihilated and pair-produced in equal rate in the early universe.However, the universe is expanding and cooling down.Although both the Hubble rate and the interaction rate are decreasing, at a given temperature, the Hubble rate overcomes the interaction rate, and freeze-out takes place fixing the DM relic abundance.
In this paradigm, the evolution of the DM number density n DM is described by the Boltzmann equation, where Y DM = n DM /s is the comoving number density, with representing the entropy of the primordial plasma, with g ⋆s being the relativistic degrees of freedom that contribute to the entropy, and x = m DM /T is a "time" variable that helps us to simplify the integration and physical interpretations.In this parametrization, the Hubble expansion rate and comoving abundance are written as, where g ⋆ accounts for the relativistic degrees of freedom, with g ⋆s ≈ g ⋆ = 106.75 at the time of freeze-out, g DM corresponds to the DM degrees of freedom, and K 2 (x) is the modified Bessel function.After solving Eq. ( 17), we obtain the DM abundance, where we use the limit of x → ∞ in the solution of the Boltzmann equation.The most current value is given by Planck collaboration, Ω DM h 2 = 0.1200 ± 0.0012 within 68% C.L. [45].
In Section II, we described how a vector-like Dirac fermion χ can be the DM candidate in a minimal B − L model.We computed the DM relic density using micrOMEGAs [46,47].In Fig. 2, we present the behavior of the relic abundance as a function of the DM mass.We have set v s = 7 TeV.The red and green solid curves are the relic abundance for g BL = 0.45 and g BL = 0.80 (corresponding respectively to m Z ′ = 6.3 TeV and m Z ′ = 11.2TeV, according to Eq. ( 4)).The gray horizontal line delimits the region that reproduces the Planck data [45].
It is remarkable that the largest part of the parameter space of the thermal relic abundance is overabundant.It occurs because, in general, the model provides small annihilation crosssections.However, at the resonance regime, when mχ ≈ m Z ′ /2, it reaches sizeable values via the s-channel in Fig. 1 (a).Such an enhancement brings the relic density down to the observed value.Because of the (inverted) resonance peak, we see that there are typically two viable DM masses: one slightly smaller than m Z ′ /2, the other slightly larger.
Thus far, we have addressed how to produce dark matter and GWs.In what follows, we will show that our model is amenable to collider and direct detection constraints.

V. THE CONSTRAINTS
Let us now turn to a discussion of the limits on the Z ′ mass from the ATLAS and LEP-II results, and direct detection bounds on the spin independent (SI) cross-section as reported by XENONnT and LUX-ZEPLIN (LZ).pairs as for example dileptons [12], dijets [13] and di-top [14].

A. Collider limits
We will adopt the most restrictive limits to our model, which stem from searches for heavy dilepton events conducted by ATLAS collaboration during Run 2 of the Large Hadron Collider, and corresponds to an integrated luminosity of 139f b −1 [12].It is well-known that the presence of Z ′ couplings to DM pairs can weaken the latter constraints [48], however, such invisible decay yields no effect on the dilepton bound, because of its small contribution to the total Z ′ width [26].ATLAS collaboration has not searched for a B-L Z ′ boson.Therefore, we had to derive our own limit by comparing the theoretical prediction with the upper limit derived at 95% CL on the fiducial cross-section times branching ratio quoted in [12].
For each combination of gauge coupling and Z ′ mass, we computed the Z ′ branching ratio using CalcHEP [49] and fed this information into Madgraph5 [50,51] where we performed the Monte Carlo simulation adopting a parton distribution function (PDF) NNPDF23LO [52].We computed the signal event pp → Z ′ → ℓ l at √ s = 13 TeV, with ℓ = e, µ and compared our result with the public result from ATLAS Collaboration [12].Following the collaboration, we required the signal events to feature two opposite charge leptons, with transverse momentum p T > 30 GeV, and pseudorapidity |η| < 2.5.In doing so, we obtained the collider limits shown in Table .II.We emphasize that these findings represent a new result in the literature.We highlight that the limit for g BL = 0.8, the production cross-section times branching ratio falls outside the sensitivity of ATLAS data.This happens because when we increase the gauge coupling, we also increase the production cross section of Z ′ bosons.This shift upward in the crosssection makes the theoretical prediction cross the exclusion limit from ATLAS at larger Z ′ masses.As the highest pole mass probed by ATLAS was 6 TeV, when we significantly increase the gauge coupling, our model can no longer be excluded by it.For this particular case, we applied the Collider Reach β tool [53] which allows us to forecast the new bounds on the Z ′ mass for a different collider configuration.Here we maintained √ s = 13 TeV, and simply ramped up the luminosity and selected to be L = 300f b −1 .Doing so, we project the ATLAS bound of m Z ′ > 7.7 TeV.This is a simple estimation of the ATLAS reach.For g BL = 0.8, we could have adopted a more solid and weaker bound, m Z ′ > 5.6 TeV, from the Large Electron Positron (LEP) collider that will address below.An old and relevant collider bound comes from LEP-II data that reads m Z ′ /g BL > 7 TeV [54,55].This limit is not af-fected by an eventual presence of an invisible branching fraction.It stems from the comparison between the SM prediction and new physics prediction for dilepton events, rather than from resonance searches.As LEP featured fantastic precision due to the leptonic nature of the process, they were able to obtain a stringent limit on the Z ′ mass.LEP bound is more relevant for g BL ∼ 1.We explore a setup with g BL = 0.8, and as we discussed previously, current data from ATLAS cannot probe this case.For this reason, for g BL = 0.8, we use the limit from LEP which reads 5.6 TeV.

B. Direct detection bounds
Since the Z ′ boson features vector coupling with both quark and DM pairs, spin-independent interactions arise between the latter and the nucleons.The corresponding cross-section can be written as with µ χN = mχm N mχ+m N being the DM-nucleon reduced mass.Such interactions are strongly constrained by xenon-based direct detection experiments.For our study, we will consider the most recent bounds as given by the LZ [11] and XENONnT [10].Notice also that the vector coupling with SM fermions implies an s-wave dominated annihilation cross-section into SM fermions, hence potentially testable via indirect detection.One could then use searches of γ-ray signals, see e.g.[56][57][58][59], to further constrain the parameter of the model.Indirect detection constraints are, however, not competitive with direct detection and collider for the model under scrutiny.Consequently, we will not show them explicitly.That said, we will now put our findings into perspective.

VI. RESULTS
Our main results concerning the DM phenomenology are displayed in Fig. 3 for g BL = 0.45 (top panel) and g BL = 0.80 (bottom panel).The red (top) and green (bottom) curves yield the correct relic abundance, in agreement with Planck's observations.The region in between the curves gives an underabundant relic, whereas the outside region is an overabundant one.Furthermore, notice that one value of Z ′ mass is associated with two values of DM mass that yield the correct relic density.This happens because the model reaches Ω χ h 2 ≈ 0.11 in the Z ′ resonance regime.
We also exhibit the direct detection constraint from the LZ collaboration that is slightly stronger than the XENONnT one.The shape of the experiment curve can be understood looking at Eq. (22).When the dark matter mass is much larger than the nucleus mass, the scattering cross-section is independent of the dark matter mass, but decreases with m 4 Z ′ .As the experimental limit linearly weakens with the dark matter, the viable parameter space in the m χ vs m Z ′ plane is mostly sensitive to the Z ′ mass.We are plotting our findings in a Log-Log scale, thus what we see is a line representing the direct detection bound weakening as the Z ′ mass increases.
As for the collider constraint from ATLAS, m Z ′ > 6 TeV.It is simply a vertical line on the Z ′ mass once we fix the gauge coupling.This has to do with the fact that the signal scales with g 2 BL × BR(Z ′ → l l).Hence, once we fix g BL the bound on the Z ′ can be directly extracted from the collaboration report.The branching of the Z ′ into dileptons will not change when m Z ′ > 2m χ , because the partial width of Z ′ → χ χ is also controlled by g BL .Therefore, when we open the invisible decay into dark matter, no meaningful change in the branching into lepton is expected.In models, where the Z ′ coupling to dark matter is different from the Z ′ coupling to leptons, a large decay into invisible can be more pronounced, however, [60][61][62].The region between the two curves corresponds to underabundance and is in principle allowed, while outside the curves the Universe would overclose.The shaded blue region represents the parameter space excluded by direct detection.In orange, we have the corresponding collider limits, which substantially constrain the model.
Notice that in Fig. 3, v s varies because the Z ′ mass changes linearly with it.We are bringing this to the reader's attention because the value of v s will be rather relevant to the GW spectra.For this reason, we drew a dashed vertical line representing v s = 7 TeV, which realizes the GW spectra to be explored later in Fig. 5.In the same way, we benchmark with blue and violet stars DM masses associated to those GW spectra.In the top panel, m χ = 2.9 or 3.4 TeV, while in the bottom one m χ has to be 5.0 or 6.15 TeV.The shift upward in the relic density curve has to do with the Z ′ width, which is larger as it grows with g BL .In Fig. 4 we fix v s = 7 TeV and obtain the relic density curve in the m χ vs g BL plane, whereas m Z ′ varies.The shaded region is ruled out by the LZ collaboration.From Eq. ( 22) we notice that the DM-nucleon SI scattering crosssection will depend only on m χ , since the dependence on g BL enters only in the ratio (g BL /m Z ′ ) 4 , which is constant for fixed v s .As we are considering m χ ≫ m N , the scattering cross-section no longer depends on the dark matter mass, and is consequently constant in the plane of Fig. 4. Although, the experimental limit from direct detection experiments linearly weakens with the dark matter mass.In particular, we found a dark matter-nucleon scattering cross section of 3.1 × 10 −10 pb, which is consistent with the experiment limit that reads 3.2 × 10 −10 pb for m χ ≃ 1.35 TeV.Therefore, for larger dark matter masses, the direct detection bound becomes too weak to constrain the parameter space in Fig. 4.This explains why the exclusion region from direct detection represents a horizontal blue curve.Collider bounds are shown following Table II.The above discussion highlights how DM searches can shed light on the properties of the DM particle χ, the Z ′ mediator, and the gauge coupling g BL .But these searches are basically insensitive to the properties of the scalar ϕ s , such as its mass or self-coupling λ s , which do not affect the DM density.Moreover, if m ϕs < ∼ 2m Z ′ (i.e.λ s < ∼ 16g 2 BL ), the scalar will only decay to the invisible right-handed neutrinos and will be hardly detectable at colliders.Fortunately, in this case, one might constrain the scalar sector using gravitational waves, and an eventual detection of GWs would allow us to determine the scalar self-coupling.This is shown in Fig. 5, where we display the GW spectra for the benchmark value of v s = 7 TeV with g BL = 0.45 (top) and g BL = 0.80 (bottom).The different spectra in each figure correspond to varying λ s , whose corresponding values are shown along each curve.Also shown are sensitivity curves for LISA, DECIGO and BBO [63].They have been constructed by taking into account the expected shape of the signal from a first order first transition, such that if the predicted spectrum lies anywhere above the curve, the expected signal-to-noise ratio is larger than 1.Notice that decreasing λ s generates larger spectra.This is because smaller λ s will lead to a smaller energy difference in the zero-temperature potential between the broken and unbroken minima (cf.Eq. ( 10)), which leads to stronger phase transitions [64].Eventually, the energy gap is so small that the transition to the true minimum never occurs: the field remains trapped in the metastable false vacuum and the symmetry is not broken, which is obviously non-physical.Hence, there is an upper bound on the spectrum that could be achieved for each pair (v s , g BL ).
Varying v s will shift the peak frequency, but since the plot is logarithmic in f , one would need large differences in v s for this shift to be noticeable.The differences in the spectra will not be significant for the range of v s values considered here (cf. the allowed region of DM production in Fig. 3).Fig. 6 shows the maximum values of λ s , for given g BL and fixed v s = 7 TeV, that would yield a detectable spectrum at BBO, DECIGO and LISA, respectively.This is computed by finding the value of λ s at which the predicted spectrum intersects the sensitivity curve at any point.In other words, the parameter space below the curves could be probed using gravitational waves.We also show the regions where the unbroken state is metastable (red region) and stable (gray region).For concreteness, considering BBO sensitivity, the region that lies between the dotted curve and the metastability yields a detectable gravitational wave signal.Thus, the detectability at LISA is achieved very close to the metastable situation.This means that, for the perturbative range of g BL shown in the plot, one can find spectra within the LISA sensitivity, but not by a huge margin.More sensitive detectors, such as DECIGO and BBO, would be able to probe a larger range of λ s values.

VII. DISCUSSIONS
Figure 5 shows that the GW spectra are highly sensitive to λ s and to g BL .There is also a dependence on v s (and hence on m Z ′ ), which is not shown explicitly in the figure, but can be understood from the fact that this parameter sets the characteristic energy scale of the transition and therefore determines the characteristic frequency of the spectral peak.Hence, from a measurement of the GW spectrum alone, one cannot determine any of these parameters separately.On the other hand, g BL , m Z ′ and m χ govern the dark matter phenomenology, and could be probed by colliders and direct detection experiments, as shown in figures 3 and 4. Therefore, one really needs to use colliders, direct detection and GW experiments to gain access to multiple directions of the parameter space independently.
It should be mentioned that the GW spectra shown in Fig. 5 are subject to theoretical uncertainties due to the 1-loop approximation used for the effective potential.The various sources of uncertainties have been studied in ref. [65] in the context of a ϕ 6 effective theory, where it was shown that the theoretical prediction for the GW amplitude at 1-loop could vary up to 2 orders of magnitude for that model.Taking this figure at face value, disregarding the striking differences between the models, we conclude from Fig. 5 that this would imply an uncertainty of O(10%) in our mapping from selfcouplings to GW spectra.In this paper, we do not aim at such level of precision.However, once such GWs are detected, theoretical accuracy will definitely be desired.For further discussion on this matter, see also [66][67][68].

VIII. CONCLUSIONS
In this work, we have argued that GW detectors have to be allied to collider and direct detection experiments in order to be able to probe the full parameter space of models with a dark sector.Dark scalars are typically extremely difficult to detect at colliders, and measuring their self-couplings are unimaginable with the machines we will build in the near-future.How- TeV, that would yield a GW spectrum detectable at BBO (dotted curve), DECIGO (dashed) and LISA (solid).The red region shows the values of λs that lead to a metastable false vacuum, whereas the gray region corresponds to stability of the unbroken state, i.e. the broken phase not being the global minimum of the potential, as per Eq.(10).
ever, the GW spectrum from a first order phase transition is highly sensitive to this scalar sector, and depend on the scalar self-coupling as well as on the gauge and fermionic sectors.
Our main point is that, in order to disentangle the parameter dependence on these observables, one must consider collider, direct detection and gravitational wave experiments.We illustrate this point concretely in a minimal B−L model with a viable DM candidate.The DM candidate is a vectorlike fermion and the U (1) B−L symmetry is broken by a scalar singlet and induces a first-order cosmological phase transition.Interestingly, we find that the expected LISA sensitivity is barely enough to probe the most extreme cases, when the phase transition is extremely strong and the unbroken vacuum is on the verge of becoming metastable.On the other hand, more sensitive interferometers such as DECIGO and BBO might be able to probe a larger range of self-couplings.There is arguably some degree of tuning for gravitational wave detection, in the sense that an exiguous change in the value of λ s might move the spectrum from within the detectability range of LISA down to non-detectability even at BBO.But this also means that a gravitational wave detection, allied with measurements of v s and g BL from other experiments, would lead to a measurement of the scalar self-coupling with good precision.
Colliders can probe the Z ′ mass and the gauge coupling g BL , but are insensitive to details of the scalar singlet.Direct detection experiments rule out a larger fraction of the parameter space, which helps us constrain the dark sector.Knowing that, if we impose a gravitational wave detection in future probes and thermal production of dark matter, we can predict which dark matter masses reproduce the correct relic density in agreement with the data.In summary, our main result is that gravitational wave detectors offer a complementary and orthogonal probe to dark sectors, allowing us to further narrow down the parameter space of dark matter models.

Figure 1 .
Figure 1.The main channels for DM-SM interactions in this B − L DM model: (a) DM annihilation into SM fermions f (mχ < m Z ′ ); (b) DM annihilation into on-shell pair of Z ′ (mχ > m Z ′ ); and (c) DM-nucleon scattering process, where q stands for quarks.

Figure 2 .
Figure2.DM relic abundance as a function of the DM mass.We set vs = 7 TeV.The red solid curve is the relic abundance for gBL = 0.45 and m Z ′ = 6.3 TeV, while the green one represents gBL = 0.80 and m Z ′ = 11.2TeV.The gray line is the abundance in agreement with Planck observations[45].Notice that there are two viable DM masses.

Figure 3 .
Figure 3. Combined results for DM phenomenology.The solid curves correspond to the correct DM relic abundance Ωχh 2 ≈ 0.12, for gBL = 0.45 (top) and gBL = 0.80 (bottom).The region between the two curves corresponds to underabundance and is in principle allowed, while outside the curves the Universe would overclose.The shaded blue region represents the parameter space excluded by direct detection.In orange, we have the corresponding collider limits, which substantially constrain the model.

Figure 4 .
Figure 4.The observed DM relic for vS = 7 TeV is represented by the black curve.The shaded orange and blue regions are the ATLAS and LUX-ZEPLIN limits, respectively.The ATLAS bound reproduces Table II.Direct detection imposes mχ > 1.35 TeV.

Figure 5 .
Figure 5. Spectra of GWs for vs = 7 TeV.The figure at the top (respectively bottom) refers to gBL = 0.45 (resp.gBL = 0.80).The various spectra displayed are obtained by varying λs at values shown in each curve.Sensitivity curves for LISA, DECIGO and BBO have been obtained from Ref. [63].

Figure 6 .
Figure 6.Maximum values of λs, for given gBL and vs = 7 TeV, that would yield a GW spectrum detectable at BBO (dotted curve), DECIGO (dashed) and LISA (solid).The red region shows the values of λs that lead to a metastable false vacuum, whereas the gray region corresponds to stability of the unbroken state, i.e. the broken phase not being the global minimum of the potential, as per Eq.(10).

Table II .
[12]a gauge bosons with unsuppressed coupling to the SM fermions, as the one considered in this work, produce strong collider signals in the form of resonances decaying into SM Lower mass bounds derived on the Z ′ mass for different gauge couplings using ATLAS results reported in[12].