Phase Transitions and Gravitational Waves in a Model of Z 3 Scalar Dark Matter

: Theories with more than one scalar field often exhibit phase transitions producing potentially detectable gravitational wave (GW) signal. In this work we study the semi-annihilating Z 3 dark matter model, whose dark sector comprises an inert doublet and a complex singlet, and assess its prospects in future GW detectors. Without imposing limits from requirement of providing a viable dark matter candidate, i.e. taking into account only other experimental and theoretical constraints, we find that the first order phase transition in this model can be strong enough to lead to a detectable signal. However, direct detection and the dark matter thermal relic density constraint calculated with the state-of-the-art method including the impact of early kinetic decoupling, very strongly limit the parameter space of the model explaining all of dark matter and providing observable GW peak amplitude. Extending the analysis to underabundant dark matter thus reveals region with detectable GWs from a single-step or multi-step phase transition.


Introduction
Although dark matter (DM) constitutes 27% of the energy density of the Universe [1], its nature is still unknown.The discovery of the Higgs boson, over a decade ago, showed that scalars do exist in nature [2,3].In particular, there are two well-motivated candidates of scalar dark matter: the inert doublet [4][5][6][7] and the scalar singlet [8][9][10], usually stabilised by a Z 2 parity.Both these DM candidates suffer, however, in the face of the most recent results from the direct detection experiments XENON1T [11], PandaX-4T [12] and LZ [13] that have this far not seen any evidence of a dark-matter signal.
A possible solution lies in considering a larger scalar dark sector comprising both the inert doublet and a complex singlet with a stabilising symmetry more complex than the Z 2 parity.In this case, dominant processes that determine the DM relic density can belong to semi-annihilation [14][15][16][17] instead of the usual annihilation.Since semi-annihilation and direct-detection processes are driven by different couplings, the correct DM relic density may be obtained while the spin-independent direct-detection cross section lies below the current limits of direct-detection experiments.
In fact, semi-annihilation is possible already for the next-to-simplest stabilising symmetry, the Z 3 group.In the Z 3 complex singlet-inert doublet model, the impact of semiannihilations on DM relic density was first considered in [18] and a comprehensive study of the model was performed in [19]. 1 On the other hand, the Z 3 singlet model [33,34] 2 and the inert doublet model [4][5][6][7] can be considered as limiting cases of the full model.The Z 3 singlet model was reconsidered with improved bounds and impact of early kinetic decoupling in [41].Indirect detection signals have been considered in [42][43][44] and cosmic phase transitions were considered in [45].A fit of the Z 3 singlet model was made in [46].The Z 3 singlet-doublet model also shares many constraints with the inert doublet model (IDM) [47][48][49][50].Recent studies of DM phenomenology in the IDM are given by [51][52][53].The phase structure of the IDM was studied in [54].First-order phase transitions in the IDM were studied in [55][56][57][58][59] and the GW phenomenology of the inert doublet model was elucidated in [60].IDM, augmented by higher-order operators, can also given rise to baryogenesis [61].Notably, Z 3 symmetry and semi-annihilation have also been considered for a singlet-extended type II seesaw model [62].
In the part of the parameter space where semi-annihilation dominates, one has to take extra care in calculating the DM relic density.In the standard calculation of the thermal relic abundance [63] it is assumed that, at the freeze-out, DM is still in equilibrium with the heat bath.Indeed, usually the elastic scattering processes with the thermal bath particles are much more efficient than the processes of annihilation and production.If the latter are enhanced, however, or if the scattering processes are not related to the number changing ones -which is the case for semi-annihilation -this assumption may not be satisfied.It has been shown that kinetic decoupling can begin as early as the chemical one and that the change in the DM phase space distribution may modify the relic abundance by more than an order of magnitude [64][65][66].
Although direct detection of DM remains elusive, other observational channels are becoming available.The discovery of gravitational waves (GW) by the LIGO experiment [67,68] has begun a new era of GW astronomy.Of particular interest to particle physicists is the possibility to probe the vacuum structure of the scalar potential via cosmic phase transitions.In the SM, the phase transition is a smooth cross-over [69,70] which does not produce a GW signal.A larger scalar sector in an extension of the SM, however, can give rise to first-order phase transitions [71][72][73] which can produce a stochastic GW background that could be detected by future GW detectors such as LISA [74,75], BBO [76,77], Taiji [78,79], TianQin [80] or DECIGO [81,82] (see e.g.Ref. [83] for a pedagogical review).
The aims of this paper are therefore twofold: first of all to study gravitational-wave signals originating during phase transitions predicted by the inert doublet-singlet model with a Z 3 symmetry and, second, to update the phenomenological analysis of this model by including constraints from vacuum stability, unitarity and perturbativity, and electroweak precision parameters and calculating the thermal relic density taking into account the effects of early kinetic decoupling.
The paper is organised as follows.We describe the tree-level scalar potential of the model together with the loop and thermal corrections in Section 2. Various theoretical and experimental constraints are given in Section 3. Relic density calculation with early kinetic decoupling is discussed in Section 4. Calculations of phase transitions and the gravitationalwave signal are elucidated in Section 5. Our Monte Carlo scan is described in Section 6 and the signals in Section 7. We conclude in Section 8. Appendices list field-dependent scalar masses (A), effective potential counter-terms (B) and formulae for the sound-wave contribution to the gravitational-wave signal (C).

Scalar potential and parametrisation
Our scalar fields comprise, besides the Standard-Model Higgs boson H 1 , an inert doublet H 2 and a complex singlet S. Dark matter is made stable thanks to a Z 3 discrete symmetry under which the Higgs boson, together with all other SM fields, transforms trivially: The most general scalar potential invariant under the Z 3 is given by [18,19] where we take all the parameters real.We can parametrise the fields in the EW vacuum in the Landau gauge as where v = 246.22GeV is the vacuum expectation value (VEV) of the SM Higgs field H 1 at zero temperature.
Notice that the absence of the term -forbidden by the Z 3 symmetry -means that there is no mass splitting between the scalar component H and pseudoscalar component A of the neutral part of H 2 .The µ SH term in the above potential induces a mixing by an angle θ between S and the neutral part of H 2 .In terms of the complex mass eigenstates x 1 and x 2 , it yields Because x 2 is doublet-like and has gauge interactions with the Z boson, it cannot be a DM candidate due to its large direct detection cross section.Therefore, we always consider M x 1 < M x 2 and choose the singlet-like x 1 as our DM candidate.
We then express the parameters in Eq. (2.1) in terms of the physical masses M h , M x 1 , M x 2 , M H ± , the dark sector mixing angle θ and the Higgs VEV v: (2.4)

Tree-level potential for phase transitions
When dealing with phase transitions and gravitational-wave signals in Subsection 7.2, we consider the path to tunnel from the origin to the EW minimum to lie in the threedimensional field space of h, H and s.Therefore in the potential (2.1) we set all other scalar fields to zero.The resulting tree-level potential is thus (2.5)

Quantum corrections to the potential
The Coleman-Weinberg potential contains the quantum corrections to the tree-level potential (2.5) given by the sum of the one-particle-irreducible diagrams with external lines from amongst the classical background fields h, H and s.The one-loop Coleman-Weinberg potential is expressed in the Landau gauge as [84] ∆V (h, where the sum ranges over the fields i = W ± T , W ± L , Z T , Z L , γ L , t, h, H, s, G 0 , A, χ, G ± , H ± .We neglect the lighter fermions because their Yukawa couplings are very small as compared to the top Yukawa.We choose the renormalisation scale to be µ = v, and C i are constants peculiar to the renormalisation scheme.The subscripts T and L, respectively, denote the transverse and longitudinal components of the gauge bosons.The bosonic and fermionic degrees of freedom3 are given by n The constants C i = 3/2 for scalars, fermions and longitudinal vector bosons, as well as C i = 1/2 for transverse vector bosons in the MS subtraction scheme.The field-dependent scalar masses m 2 i ≡ m 2 i (h, H, s) correspond to the eigenvalues of the field-dependent mass matrices of the scalar fields M 2 S , of the pseudoscalar fields M 2 P and of the charged fields M 2 C , which are given in the Appendix A. The field-dependent top quark and gauge boson masses are given by ) where y t , g 1 and g 2 are the top Yukawa coupling, the U (1) Y and SU (2) L gauge couplings, respectively.
In order to fix the masses, VEV and the dark sector mixing angle at their tree-level values, we cancel contributions from the one-loop ∆V in the EW vacuum with opposite-sign contributions from the counterterm potential given by whose coefficients are given in the Appendix B. This cancellation is guaranteed by the following set of renormalisation conditions: where vev means at the EW minimum (h, H, s) = (v, 0, 0).In Eq. (2.9), we take care of the problematic IR divergence occurring in ln m 2 G (due to the vanishing Goldstone mass in the EW vacuum) by replacing m G with M h [86] as an infrared regulator.

Thermal corrections to the potential
The one-loop finite-temperature contribution to the effective potential is given by [87] V where the thermal bosonic/fermionic functions J B/F are defined as [87] J It is customary to deal with IR divergences arising from the zero bosonic Matsubara mode in the high-temperature limit by considering the so-called daisy resummation [88].This procedure amounts to adding the leading-order thermal correction or Debye mass Π i T 2 to the mass m 2 i in ∆V and V T [89].For the scalar sector it means where with the SM contribution taken from [90].There is no need to consider the thermal mass of the top quark and fermions in general, because daisy diagrams do not diverge in the IR, since there is no zero Matsubara frequency for fermions [91][92][93].As for gauge bosons, only their longitudinal modes acquire a significant thermal correction to their mass, thereby the Debye mass of transverse modes (protected by gauge symmetry) is neglected [91,94,95].Notice that thermal contributions to the longitudinal part of Z and γ should be added to their mass matrix in the gauge basis before it is diagonalised.The Debye mass of to their mass matrix and then diagonalise it [91,96].

Full thermally-corrected effective potential
The full thermally-corrected effective potential needed to compute phase-transition parameters is thus made of the classical potential (2.5), to which one adds quantum corrections summarised in (2.6) and (2.8), as well as thermal corrections (2.10): (2.17) 3 Constraints

Perturbativity
We require the Feynman rule vertex factors of quartic interactions to be less than 4π in absolute value to ensure that the one-loop contributions be smaller than tree-level ones [97].This yields the conditions (3.1)

Unitarity
The cross section of 2 → 2 scattering processes s 1 s 2 → s 3 s 4 of scalars s i can be expanded in the partial wave decomposition as where s is the Mandelstam variable and a l are the partial wave coefficients with the angular momenta l.The optical theorem imposes on a l the unitarity bound In the high energy limit, the s-wave amplitude a 0 (s) is dominated by contact terms because the s-, t-and u-channel processes are suppressed by the scattering energy.Therefore, in the high energy limit, a 0 (s) is completely determined by the quartic couplings of the scalar potential.The constraint (3.3) places an upper bound on the eigenvalues of the 2 → 2 scalar scattering matrix.The unitarity bounds of the 2HDM were first studied in [98,99].
For the Z 3 scalar quartic couplings, the unitarity bounds are given by4 ) ) ) ) ) and |Λ i | ≤ 1/2, where the three remaining eigenvalues Λ i with i = 1, 2, 3 of the scattering matrix are given by the solution of the cubic equation

Bounded-from-below constraints and globality of the electroweak vacuum
In order to be physical, the potential has to be bounded from below in the limit of large field values.In this limit, the dimensionful quadratic and cubic terms can be neglectedit suffices to consider only the quartic part of the potential (2.1).The analytical necessary and sufficient bounded-from-below constraints were derived in Ref. [100].We only sketch the derivation and refer the reader to that paper for details and ancillary Mathematica code.
The potential (2.1), using the parameterisation where 0 ≤ ρ ≤ 1 as implied by the Cauchy inequality 0 where we have minimised cos(ϕ + 2ϕ S + ϕ λ S12 ) = −1 so λ S12 = −|λ S12 | without loss of generality.We have to minimise the potential with respect to ρ, also taking into account separately the ends ρ = 0 and ρ = 1 of the interval.The conditions for s = 0 or ρ = 0 are given by Another necessary but not sufficient (since it is calculated with ρ = 1, but ignoring the λ S12 term) condition can be given by We minimise the potential with respect to h 1 , h 2 , s and ρ with the fields lying on a sphere, enforced by a Lagrange multiplier λ.The minimisation equations with respect to the fields and λ are and the minimisation equation with respect to ρ is The equations (3.19) and (3.20), solved together, have an analytical solution.The boundedfrom-below conditions for 0 < ρ < 1 are then given by where V min , the value of minimum at the solution, is proportional to the Lagrange multiplier λ.Note that p =⇒ q is equivalent to ¬p ∨ q.In addition, the equations (3.19) need also to be solved with ρ = 1, in which case the conditions are given by5 The necessary and sufficient bounded-from-below conditions are then given by Eqs.(3.17), (3.21) and (3.22).Eq. (3.18) can be used for faster elimination of invalid points.
Once the bounded-from-below conditions are satisfied, the potential is guaranteed to have a finite global minimum solution.Although it is possible that our EW vacuum, not coinciding with the global one, is metastable, the gain in the parameter space is expected to be marginal (at least in the singlet-like case [41]).For that reason, we simply demand that the EW vacuum be the global one.

Collider constraints
The LEP collider measured to high precision the decay widths of the Z and W bosons, which are compatible with the SM.Neglecting the small mixing between the singlet and doublet, we require that the particle masses satisfy [101] LEP searches for new neutral final states further exclude a range of masses [102], thereby forcing in addition to m H ± > 70 GeV (3.25) due to searches for charged scalar pair production [103].But since there is no splitting between m H and m A , the constraint (3.24) is always satisfied.
Similarly, if M x 1 < M h /2 or M x 2 < M h /2, the Higgs boson can decay into dark sector.This is constrained by measurements of the Higgs boson invisible branching ratio BR inv = Γ h→x 1 x 1 /(Γ h→SM + Γ h→x 1 x 1 ).Latest measurements by the CMS experiment at the LHC find that BR h→inv < 0.15 [104], while the ATLAS experiment measures BR h→inv < 0.11 [105], of which we will require the stricter ATLAS constraint.
We also use the micrOMEGAs package to compute the h → γγ decay rate (usually similar to the IDM one [47,48]) and require it to be within the value measured at the CMS [106]: µ h→γγ = 1.12 ± 0.09.

Electroweak precision parameters
The measurements of electroweak precision data at the LEP put strong constraints on physics beyond the SM.The latest electroweak fit by the Gfitter group [107] gives for the oblique parameters S and T (with U = 0) the central values with a correlation coefficient of +0.92.
To calculate electroweak precision parameters S and T , we use the results for general models with doublets and singlets [108,109].The usual loop functions are defined as with F (I, I) = 0 in the limit of J → I, and where The beyond-the-SM contributions to the S and T parameters are given by and where the approximations are given for a small mixing angle.The electroweak precision constraints mean that M H ± cannot be too different from M x 2 .We require the S and T parameters to fall within the 95% C.L. ellipse.

Relic density
The cosmological abundance of DM measured by the Planck satellite Ω c h 2 = 0.120 ± 0.001 [1] significantly constrains the allowed parameter space of the model.The model that we study contains a non-trivial dark sector that is comprised of several particles with Z 3 quantum numbers: x 1 , x 2 and H ± .Among them only the singlet-like x 1 can play the role of DM, since the physical states that originate from the inert doublet H 2 are too strongly coupled to the SM fields to satisfy the constraints and requirements for the DM candidate.Nevertheless, the processes that involve other dark sector particles in the early Universe can leave an imprint on the relic abundance and generally has to be included in the calculation of DM density evolution.In similar models with semi-annihilation [18,19,26,110], the relic density is obtained under the assumption that all the particles were initially in thermal equilibrium with the SM plasma and furthermore that the kinetic equilibrium is maintained throughout the whole DM production stage.In this case the density of DM can be obtained by solving the set of Boltzmann equations for the densities of different particles in the dark sector, e.g. using available numerical solvers like micrOMEGAs [111][112][113] that can work with multi-component DM models [114].If the processes of DM scattering on the SM plasma particles are not very efficient, the kinetic decoupling can occur before the freezeout.In this case semi-annihilation processes can redistribute the energy of DM particles and modify its temperature, 6 which has an impact on the rates of annihilation processes.Although this assumption is often justified or introduces very small corrections to the relic density, in some regimes the process that defines the dynamics of the freeze out can be strongly velocity-dependent, e.g.due to a resonance or s-wave suppression, so that its rate becomes sensitive to the temperature of DM.The density evolution of Z 3 scalar singlet DM with semi-annihilation beyond the assumption of kinetic equilibrium was studied in Ref. [41].In this work we use the same method to calculate the relic density, but extend it to the multi-component case of the model that we study and include additional interactions present in this model.
In order to properly take into account the departure from kinetic equilibrium of the DM particles and its effect on the relic abundances, we follow the formalism of [64] and the numerical implementation as presented in the DRAKE code release paper [66], with the addition of semi-annihilation collision term as discussed in [41].We focus on the cBE method, i.e. coupled Boltzmann equations for the 0 th and 2 nd moment of the distribution function, describing the evolution of number density and temperature, respectively.In doing so we neglect potential corrections from the departure of the thermal shape of the DM distribution. 7Below we briefly review the adopted formalism.
The formulation of the system of Boltzmann equations that we solve originates from the master equation governing the evolution of the phase-space density f i = f i (t, p) of particle i in an expanding Friedmann-Lemaître-Robertson-Walker universe: where where we introduce the moments of the collision term as ) This set of cBE (4.2) and (4.3) is not closed in terms of higher moments unless some assumption is made on the form of the distribution function.A well motivated one is to make an ansatz which enforces the Maxwell-Boltzmann shape, but with an arbitrary temperature T i .In the scenario studied in this work we have i ∈ {x 1 , x 2 , H ± }, but we will make a simplifying assumption that only the DM particle, x 1 , has a population which departs from kinetic equilibrium, while all three dark sector states can depart from chemical equilibrium.This assumption is justified as x 2 and H ± are relatively tightly coupled to the thermal bath, leading to their kinetic decoupling happening well after all the freeze-out processes take place, while in the regions of dominant semi-annihilation the x 1 has suppressed elastic scatterings and thus can undergo an early kinetic decoupling.The general expression for the collision term for i-th particle due to binary collisions reads where the sum goes over all possible j and two-particle {kl} states, with and without the DM, and M 2 is the amplitude squared summed over initial and final internal d.o.f.It is understood that {kl} = {lk} and it appears in the sum only once.The symmetry factor (1 + δ kl ) −1 is related to the exchange of identical particles in the final state.This formula incorporates annihilations, semi-annihilations and inelastic scatterings.Following [116] and [117], we treat the elastic scatterings in the Fokker-Planck approximation, effectively involving expansion in small momentum transfer.
Including also all possible decays of x 2 and H ± , and neglecting the quantum degeneracy effects in all the cases that we consider, such that (1 ± f ) ≈ 1, the cBE system for the i-th particle has the form: where H ≡ H/ [1 + g(x)] and g ≡ 1 3 T g s eff dg s eff dT with g s eff being the number of entropy degrees of freedom, the sum goes over all possible j states and {kl} two-particle states, K 1 is the modified Bessel function of order 1, Γ j→ikl stands for the decay width of the j-particle into ikl states and the rates are defined below. 8This is a set of 4 coupled ordinary differential equations that describe the evolution of comoving number densities of x 1 , x 2 and H ± , as well as the temperature of x 1 .
The evolution of the above system is governed by the interaction rates in the form of the momentum transfer rate γ(x), as defined in [64], and integrals I, J, K defined as follows.The number-changing processes are described by the following integral 9   I ij→kl = 1 n eq i n eq j dΠ i dΠ j g i g j W ij→kl f eq i (T i )f eq j (T j ), (4.11) through the use of the dimensionless invariant rate [119], that is related to the cross section via W ij→kl = 4p ij √ sσ ij .The integrals involved in the temperature evolution, describe the change in energy due to forward processes (annihilations), backward processes (inverse annihilations) and decays, respectively.Numerical solution of the above cBE system is relatively CPU expensive, in large part due to necessity of calculating a large number of rates. 10To render the task manageable, the numerical procedure we follow is to first determine for each of the parameter point if the full cBE system is needed or if the standard approach assuming kinetic equilibrium for the DM particle should suffice.For that we impose a rather conservative condition that if γ(x)/H(x)| x=100 ≥ 10, then for the whole period of chemical freeze-out the elastic scatterings are assumed to be efficient enough to keep the DM in kinetic equilibrium, even in the presence of other disruptive processes like decay of a heavier dark sector state, and the assumption of T x 1 = T is well justified allowing for calculations using only the 9 This is a generalization of the thermally averaged cross section for the case of (co-)(semi-)annihilating particles with different temperatures.Indeed, for Ti = Tj ≡ T ′ this expression is the same as in Ref. [119].Following their evaluation one finds where we used pij = s − (mi + mj) 2 s − (mi − mj) 2 /(2 √ s). 10 Note that in contrast to annihilation processes, the inelastic scatterings and semi-annihilations cannot be always included with the use of inclusive cross sections, but rather one has to perform the thermal average for every exclusive process separately.
The needed rates were calculated using FeynCalc [120] for the matrix elements and using the DRAKE thermal average routines.On top of that the decay widths were utilized from the micrOMEGAs calculation.The latter was also used to include in the number density equation typically sub-dominant, but in some cases important annihilations to final states including virtual Z or W bosons, where we followed the procedure as described in the micrOMEGAs manual [114].Heavier dark sector states were assumed to take part in co-annihilations and inelastic scattering processes if their mass did not exceed 50% of the DM state and the binary collision rates were truncated for processes that contribute less than 1% at the approximate time of freeze-out.In Fig. 1 we show the thermal history of two example benchmark points for which kinetic decoupling happens early enough to have an impact on the final relic abundance.The left plot shows the evolution of the yield and how it differs between the cBE and the usual approach assuming kinetic equilibrium throughout (which we call nBE), while right plot highlights how the temperature evolution can differ from the one of the SM plasma.Both example points have strong semi-annihilation and kinetic decoupling happening at a time when the annihilations are still active, thus leading to a modification of the evolution of the DM number density.

Phase transition and gravitational waves
The vacuum in which we live today corresponds to the minimum of the scalar potential located at the broken EW phase (h, H, s) = (v, 0, 0).In the first moments of the Universe, however, the temperature was so large that thermal corrections influenced the T = 0 effective potential significantly.As one goes back in time, the temperature raises and the increasing thermal effects push the low-temperature minimum (or true vacuum) upward: it becomes higher than the minimum at the origin and eventually disappears.Thus in the limit of T → ∞, the high-temperature, symmetric phase is usually (h, H, s) = (0, 0, 0); the symmetry is restored.A more complicated possibility is symmetry non-restoration, where in this limit the minimum of the potential does not go to the origin [121,122], or inverse symmetry breaking [123], where a symmetry, not broken at low temperature, is broken at high temperature.
In our case, the story thus begins at the origin of the potential, in the high-temperature phase.Then as the temperature decreases, the potential develops a second minimum, say along the h axis, where the potential is initially larger than at the origin.With time, the temperature continues to decrease until these two phases become degenerate: the value of the potential in both vacua is identical.The temperature at which it occurs is the critical temperature T c .At some lower temperature, the field will eventually tunnel from this false vacuum to a deeper stable minimum (or true vacuum), here along the h axis.This tunneling occurs through the nucleation of bubbles of true vacuum, which expand and collide with each others, filling the Universe little by little and thereby gradually converting the false vacuum into the true vacuum.The collision between these bubbles is responsible for the generation of a stochastic gravitational-wave background.The temperature at which one bubble nucleates per Hubble volume is called the nucleation temperature T n .For a fast firstorder phase transition without large reheating, this can be considered as the temperature at which gravitational waves are produced [124]. 11Finally, the value of the Higgs VEV smoothly evolves towards its tree-level value v as the temperature approaches zero.
A couple of important parameters capture the most information about a first-order phase transition.The strength of a first-order phase transition is characterized by the vacuum energy ∆ϵ released during the phase transition.To obtain a dimensionless parameter α, one normalizes this vacuum energy to the energy density of the plasma outside the bubbles ρ rad [126]: ∂T , T * the temperature at which GWs are generated and ρ rad = π 2 30 g * T 4 , where g * is the effective number of relativistic degrees of freedom at the temperature T .Note that a first-order phase transition is characterized as a strong one if the ratio of the VEV to the temperature, evaluated at T * , satisfies If this ratio is much smaller than unity the nature of the PT cannot be distinguished without lattice calculations, as also evidenced by the numerically very small value of the would-be GW signal in this case [69,127,128].The inverse time duration of the PT, β, is usually normalized by the Hubble parameter H * ≡ H(T * ) [129]: where S ≡ S(h, H, s, T ) is the 3-dimensional Euclidean action computed for an O(3)symmetric critical bubble.A large β/H * means that we can safely neglect the Hubble expansion of the Universe during the phase-transition process.
A single bubble cannot generate gravitational waves, because its quadrupole moment vanishes due of its spherical symmetry.The key ingredient for a stochastic GW background to be generated is the collisions between these bubbles of true vacuum, as mentioned above.Indeed, this will break the spherical symmetry of each of the collided scalar-field bubbles [130].For bubbles evolving in a thermal plasma, additional contributions to the GW signal come from the metric perturbation induced by the plasma: sound waves [131,132] and magnetohydrodynamic (MHD) turbulence [133].Sub-and supersonic bubbles, respectively, create compression waves beyond the bubble wall and rarefaction waves behind it.When acoustic waves from different bubbles overlap, the induced shear stress of the metric generate gravitational waves.Finally, bubble collision is responsible for the turbulent motion of the fully ionized plasma.These MHD turbulences will generate gravitational waves.The power spectrum h 2 Ω GW of the stochastic GW background from a cosmic phase transitions thus arises from three contributions [124]: 12Note that the treatment of cosmic phase transitions suffers from uncertainties due to gauge or renormalisation-scale dependence for instance, that could potentially lower the resulting GW signal stength [135,136].
The scalar-field contribution redshifted to today is given in the envelope approximation by [137] 13 with where κ col is the efficiency factor for the conversion of the vacuum energy into the gradient energy of the scalar field (energy stored in the shell of the scalar-field bubbles), v w is the bubble-wall speed in the rest frame of the plasma far away from the bubble [124] and f col is the peak frequency, thus the frequency at h 2 Ω peak col .
The sound-wave contribution is given by [141-143] with14 where κ sw is the efficiency factor for the conversion of the vacuum energy into the bulk motion of the plasma and f sw is the sound-wave peak frequency.The suppression factor for a radiation-dominated Universe is defined as [144,145] where τ sw is the time after which sound waves do not source the GW production and where Υ ≃ τ sw H when τ sw H ≪ 1, thus when this time is much smaller than a Hubble time.
Finally, the MHD-turbulence contribution is given by [124,146] with where κ turb is the efficiency factor for the conversion of the vacuum energy into the MHD turbulences, f turb is the MHD-turbulence peak frequency and h * is the today-redshifted value of the Hubble rate at T * .

Numerical scan
We used Markov Chain Monte Carlo (MCMC) with the Metropolis-Hastings algorithm in order to generate points, starting with a number of points with the DM relic density within the Planck 3σ range.We then imposed other constraints: perturbativity, perturbative unitarity, bounded-from-below conditions for the scalar potential, globality of the EW vacuum, bounds on electroweak precision parameters S, T and U , Higgs invisible width and Higgs to diphoton h → γγ decay.We performed three scans: 1. Points with non-zero quartic λ S12 coupling (and 2. Points with non-zero cubic µ ′′ S coupling (and 3. General points with non-zero µ ′′ S and λ S12 and without conditions on the quartic portals λ S1 , λ S2 , λ 3 .
In the first two cases, all the biquadratic portal couplings were set to zero, i.e. λ S1 = λ S2 = λ 3 = 0, except for λ 4 , because this would have entailed enforcing a special relation between M x 2 and M H ± .Since we need a non-zero mixing angle θ between x 1 and x 2 , the µ SH mixing term was non-zero in all cases.The dark sector scalar self-couplings were set to λ S = λ 2 = π in order to allow for larger values of λ S12 and µ ′′ S which are otherwise constrained by vacuum stability.
New points were generated from a multivariate normal distribution.The initial covariance matrix was calculated from a set of a few tens of points with DM relic density within the Planck 3σ bounds.The initial values for M x 2 and M H ± were chosen in an interval around 2M x 1 in order to enhance semi-annihilation processes.The covariance matrix was updated with each added point, while its norm was kept constant to ensure a reasonable acceptance ratio of points.Dimensionful quantities M x 1 , M x 2 , M H ± and µ ′′ S were generated logarithmically and only points with M x 1 < M x 2 , M H ± kept.Similarly to the calculation of the initial distribution, a few tens of points with DM relic density within the Planck 3σ bounds and varied couplings were used to start the MCMC chains from a few hundred up to thousand points in length.Since the desired probability distributions of our MCMCs were rather arbitrary, we do not consider the density of points in our plots to have a specific statistical meaning.
The desired distribution for each case consists of a product of Gaussians whose central values and standard deviations are given in Table 1.Extra suppression factors of an arbitrary small value 10 −12 were used to make the Gaussians essentially one-sided when the conditions sin θ > 0 and 0 ≤ λ S12 < 4 (the maximal value from vacuum stability bounds for λ S = λ 2 = π and portal couplings λ S1 = λ S2 = λ 3 = 0), and 1 ≤ µ ′′ S /GeV < 2 √ λ S M x 1 /GeV did not hold true.In the last condition, used in the second and third scans, we limited the cubic coupling by |µ ′′ S | < 2 √ λ S M x 1 , the bound from absolute vacuum stability in the Z 3 singlet limit.Notice that unlike in the singlet-only limit, the cubic µ ′′ S coupling can yield semi-annihilation together with the ever-present mixing term µ SH .

Signals
We now discuss the results of the scan, in particular for the direct-detection and gravitationalwave signals.

Direct detection
In Fig. 2 we show the predictions for the spin-independent direct detection cross section for the three cases.The top panel shows points with non-zero λ S12 , the middle panel points dominated by µ ′′ S and the bottom panel shows the points from the more general scan.Also shown are the sensitivity curves of the XENON1T (2018) [11], PandaX-4T (2021) [12], LZ (2022) [13] in grey and the projected sensitivity curve for XENONnT [147] in dashed grey.The color code shows the semi-annihilation fraction r defined by r = 1 2 where x denotes x 1 , x 2 , H ± and X stands for any SM particles.With dominant semiannihilation r approaches unity.The points with early kinetic decoupling, evaluated with the cBE method detailed in Section 4, are marked with crosses.

Gravitational waves
For each of the three scans mentioned in Section 6, we study how the Universe ends in the EW vacuum (v, 0, 0) at zero temperature.If at least one phase transition, to reach that state, is of the first order, then we subsequently compute the power spectrum of the stochastic GW background h 2 Ω GW , necessarily associated to these cosmic FOPTs.In addition to the previous constraints, we also take into account the constraint from directdetection experiments: only points not excluded by LZ are displayed in Fig. 3 and Fig. 4. Note also that we consider a relativistic bubble-wall speed, v w = 1.
In the left panel of Fig. 3, we show, for the three scans, the region of the parameter space, projected onto the plane made by the mixing angle sin θ and the DM mass M x 1 , that gives rise to first-order phase transitions.These are essentially obtained for large mixing and mild to large DM mass.We see that only points from the general scan lead to strong FOPTs.The right panel shows the values of either λ S12 or µ ′′ S , as a function of M x 1 , responsible for FOPTs for scans where they are non-zero.In both panels, triangles are points for which the cBE method has been used to compute the DM relic density, as the scattering rate is not sufficient to consider kinetic equilibrium for DM during the freeze-out process (see Section 4).  the power-law integrated (PLI) sensitivity curves of LISA, DECIGO and BBO, obtained for an observation time of 4 years and a threshold value for the signal-to-noise ratio to be 10 [148,149].Similarly to Fig. 2, we show the semi-annihilation fraction with a color code.The upper panel of Fig. 4 shows the GW signal for the scan with non-zero λ S12 .This signal results from single-step FOPTs.We can see that none of the points leads to strong FOPTs and that for some of them (triangles), the assumption of kinetic equilibrium does not hold and the full cBE system should be solved.The middle panel shows the GW signal for the scan with non-zero µ ′′ S .Again, we only find single-step weak FOPTs.A cubic term helps to generate a barrier in the potential and thus enhance the strength of the phase transition.However, we only find a single-step FOPT from the origin towards the EW minimum, thereby along the Higgs direction h.This is why a non-zero µ ′′ S is not helpful here to obtain a stronger PT, as it only impacts the s direction.Finally, the lower panel of Fig. 4 shows the GW signal for the general scan, where strong FOPTs are obtained and where one of the points clearly lies above the LISA PLI sensitivity curve, thus yielding a detectable GW signal.Note that we also found two-step phase transitions.For each of them, indicated by a triangle, the first step (0, 0, 0) → (0, v H , 0) is a second-order PT and is followed by a first-order PT (0, v ′ H , 0) → (v h , 0, 0).Note also that from Fig. 4, it seems that the lower the semi-annihilation fraction, the stronger the GW signal.However, if we do not take experimental constraints into account, the correlation between the strength of the GW signal and the semi-annihilation fraction is lost.
Finally, let us stress that in principle, in this model, if we relax the constraint on relic density and allow for underabundant DM, we can more easily find regions of the parameter space allowing for different kind of multi-step phase transitions in the three-dimensional field space (h, H, s) and/or leading to a strong enough GW signal to be detected by LISA, signal comes also with weak FOPTs whose GW signals would lie orders of magnitude below detection projections in a foreseeable future.On the other hand, a general scan of the model reveals points near the detection limits of LISA, DECIGO and BBO.
We then also computed the GW signal for uniformly distributed random general scan allowing for underabundant singlet-like x 1 with relic density below the Planck constraint.(There could be either an additional non-thermal production mechanism or another DM component.)In that case we found a sizable portion of the parameter space leading to a single-step or a multi-step phase transition producing strong enough GW signal to be potentially detected by LISA, DECIGO or BBO.
Therefore, we conclude that in the Z 3 symmetric model with an inert doublet and a complex singlet there exists a parameter range that explaining whole dark matter and providing potentially observable GW peak amplitude, though it is rather small and challenging to be observed in the near future.Nevertheless, when extending the analysis to underabundant dark matter the detection prospects of the GW peak amplitude become significantly more promising.

A Field-dependent scalar mass matrices
The field-dependent scalar masses m 2 i ≡ m 2 i (h, H, s), present in Eq. (2.6) and Eq.(2.10), correspond to the eigenvalues of the field-dependent mass matrices of the scalar fields M 2 S , of the pseudoscalar fields M 2 P and of the charged fields M 2 C : Evaluated in the EW vacuum (h, H, s) = (v, 0, 0), these mass matrices become (A.4) thus putting in evidence the mixing between H 2 and S.

B Counter-terms
We present expressions for counter-terms in the effective potential used to keep VEVs, masses and mixing at their tree-level values.where the mean adiabatic index Γ = 4/3 for a relativistic fluid, and the enthalpy-weighted root-mean-square of the fluid velocity Ūf is expressed as [143]: The mean bubble separation at percolation temperature16 R * is defined as with c s = 1/ √ 3 the speed of sound in the plasma [126].The parameter Ωgw , which quantifies how efficiently the kinetic energy is converted into gravitational waves, is numerically found to be Ωgw ≃ 1.2 × 10 −2 [141,143].
Finally, the factor F gw,0 needed to obtain the today value h 2 Ω peak sw is defined as [142,143] F gw,0 = Ω γ,0 g s0 g s * 4/3 g * g 0 , (C.4) where Ω γ,0 ≡ ρ γ,0 /ρ c,0 is the photon density today, g s * the effective number of entropic degrees of freedom, g * the effective number of relativistic degrees of freedom, while g s0 and g 0 give their present values, respectively.The degrees of freedom g s * and g * are equal at the EW scale, while g s0 and g 0 are different since they are considered after the neutrino decoupling.
The following values for the aforementioned parameters allows us to compute the factor F gw,0 [150,151]: where N eff = 3.044 [152] has been used in the computation of g s0 and g 0 .Note that g s0 ≃ 3.91 and g 0 = 2 are considered in Refs.[141,142].Considering these values, we can thus express F gw,0 as

Figure 3 .
Figure 3. Projection of the parameter space, showing points which satisfy all the theoretical/experimental constraints.cBE points that lead to FOPT are only found for λ S12 ̸ = 0 scan.Strong FOPT are only found for general scan.Left panel: Sine of the mixing angle sin θ vs. DM mass M x1 .Right panel: λ S12 vs. M x1 for both the λ S12 ̸ = 0 and general scan and µ ′′ S vs. M x1 for both the µ ′′ S ̸ = 0 and general scan.
ρ γ,0 = 0.26 eV cm 3 , ρ c,0 = 1.05 × 10 −5 h 2 GeV cm −3 , h ≡ H 0 100 km s −1 Mpc −1 , g s0 ≃ 3.94 and g 0 ≃ 3.38, (C.5) is the Hubble parameter and the collision term C[f i ] contains all possible binary and non-binary interactions between dark sector and/or SM particles.Taking the first two moments of the above equation, i.e. integrating it over g i d 3 p/(2π) 3 /E and g i d 3 p/(2π) 3 p 2 /E 2 , respectively, with g i being the number of spin degrees of freedom, and introducing new variables Y ≡ n/s and y ≡ Ewith x = m x 1 /T , one finds:

Table 1 .
The central values and standard deviations for the desired distribution of the MCMC.