Single-step first order phase transition and gravitational waves in a SIMP dark matter scenario

We investigate the non-zero temperature dynamics of a sub-GeV dark matter scenario freezing-out via self-interactions. As a prototype, we take up the case of a scalar dark matter species undergoing $3 \to 2$ number changing annihilations catalysed by another scalar. We study the shape of the thermal potential of this scenario in a parameter region accounting for the observed relic abundance. An analysis reveals the possibility of a first order phase transition with bubble nucleation occurring at sub-GeV temperatures. This finding can be correlated with the typical sub-GeV masses in the framework. The gravitational wave spectra associated with such a phase transition is subsequently computed.


I. INTRODUCTION
It has been concurred by numerous observations [1][2][3][4] that around a fourth of the total energy budget of the universe is some non-luminous matter known as dark matter (DM).The amount of such matter in the universe is also well measured from various observations, e.g., cosmic microwave background (CMB) [5] and large scale structure surveys [6].If DM is thought to be an elementary particle, then there is no information available from the current experiments on its mass and quantum numbers.More importantly, that there is no such particle candidate for DM within the Standard Model (SM), advocates the presence of additional dynamics.While DM can be hypothesized to stem from some beyond-the-Standard Model (BSM) framework, a possibility in this direction is that DM is a weakly interacting massive particle (WIMP) [2,[7][8][9][10][11].That is, the DM-SM interactions are in the weak ballpark.Moreover, the DM relic abundance in the WIMP paradigm is generated typically through 2 DM → 2 SM annihilation processes followed by a thermal freeze out.
Search for DM through its scattering with nuclei at terrestrial detectors is known as direct detection [12,13].However, non-observation of such scattering at have put stringent upper limits on the corresponding cross sections, and ultimately, on the DM-SM interactions.Such small interactions are also consistent with the null results obtained in DM searches at the energy frontier, the Large Hadron Collider (LHC) being the prime example here [14].These results have pushed a plethora of WIMP scenarios to a corner [10,11].An interesting alternative therefore could a strongly interacting massive particle (SIMP) [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29].In this case, the freeze-out dynamics is driven by DM matter number changing processes, i.e., processes that feature an imbalance of DM particles in the initial and final states.The main advantage of the SIMP paradigm is that it allows such number changing processes to be predominantly driven by DM self interactions.The DM-SM interaction stengths are thus permitted to be appropriately small so as to conform with the limits from direct detection and colliders.A SIMP therefore typically has sub-GeV mass and a large self-scattering cross section.
With the advent of gravitational wave (GW) interferometry [30][31][32][33][34][35][36][37][38], it has become possible to probe a dark sector through its GW imprints [39][40][41][42][43][44][45][46].In this study, we aim to study the possibility of a strong first order phase transition (SFOPT) [47][48][49][50] and the GW spectrum from associated with a sub-GeV SIMP framework.The model we take up is where the SM is extended by the scalars φ and δ that are singlets under the SM gauge symmetry [20].The scalar φ, stabilised by a governing Z 2 symmetry, becomes the DM candidate here.The number changing process for this case that is consistent with the Z 2 is φφφ → φh 2 , where h 2 is an admixture of δ and the Higgs doublet H.The scalar δ and its interaction with the DM φ thus play crucial roles in generating the observed relic density through the aforesaid 3 → 2 dynamics.Given the importance of δ in this setup, a pertinent question to ask is whether the same can trigger a SFOPT and how strong is the resulting GW amplitude.We plan to explore this possibility here by incorporating finite temperature corrections to the scalar potential along the direction of δ.
The study is organised as follows.We introduce the setup and outline the 3 → 2 dynamics in section II.Thermal corrections to the scalar potential are computed in section III.The same section presents the GW amplitude for representative benchmark points.We finally conclude in section IV.Some important formulae can be found in the Appendix V. We extend the SM with the real gauge singlet scalars φ and δ [20].A Z 2 symmetry is imposed under which φ → −φ while all other fields transform trivially.This negative Z 2 charge for φ stabilises the same and makes it a potential DM candidate.The scalar potential of the setup reads In the above, all parameteres are taken real and H denotes the SM Higgs doublet.The doublet H and the singlet δ receive vacuum expectation values (VEVs) v and v δ respectively.One then writes Here, v = 246 GeV.It is noted that φ does not receive a VEV and thus the Z 2 remains intact.In all, such a configuration of the VEVs introduces a mixing between h and δ .The mixing of φ is however forbidden by the Z 2 symmetry.The mass matrix in the (h δ ) basis is We have made use of the tadpole conditions ∂V ∂h = 0 = ∂V ∂δ while deriving Eq.(3).One derives the tadpole conditions to be Eq.( 3) is diagonalised through a 2 × 2 rotation parameterised by a mixing angle θ.That is, which leads to the eigenstates (h 1 , h 2 ).We identify h 1 with the discovered Higgs boson with m h1 = 125 GeV.
The DM mass is given by Some of the parameters in the scalar potential can be expressed in terms of the physical quantities of the theory such as masses and mixing angles using the relations The independent variables of the present framework thus are {m h2 , m φ , µ 21 , v h , v δ , sinθ, λ 11 , λ 12 , λ 13 , λ 22 , λ 23 }.
The relic abundance in a strongly interacting DM setup is generated through number changing processes.In this setup, the leading number changing processes consistent with the Z 2 is φφφ → φS where S = h 1 , h 2 .We list below the Feynman diagrams for illustration.
Figs.1, 2 and 3 show the 3 → 2 annihilation diagrams.We take the limit where sinθ and λ 13 are tiny.The tiny sinθ conforms with the constraints on δ − H mixing from direct search and Higgs signal strengths.The φ − φ − h 1 trilinear coupling becomes negligibly small in this limit while the φ − φ − h 2 coupling becomes (µ 21 + λ 12 v δ ).Moreover, the limit entails that the φ-mediated amplitude in Fig. 1 is the dominant one.The following is its expression: In the last step above, we have neglected the DM velocity, a justifiable assumption for non-relativistic DM.One can approximate p i (m φ , 0) in that case.The evolution of the DM comoving number density Y φ is described by the Boltzmann equation [7] where s = 2π 2 45 g * s T 3 is the entropy density and the Hubble parameter is H = 8π 3 Gg * T 4 90 1 .Besides, Y eq P denotes the equilibrium comoving number density of a particle P and σv 2 φφφ→φh2 is the thermally averaged 3 → 2 cross section [51].It can be related to σv 2 φφφ→φh2 using modified Bessel functions K n (x) as Further, Eq.( 9) becomes the following upon taking g * s g * .
We present an approximate analytical solution to Eq.( 11) to gain insight on the dynamics.We follow here the approach adopted in [21].
We take B ≡ 0.116 g φ σv 2 φφφ→φh 2 and rewrite Eq.( 11) using ∆ = Y φ − Y eq φ as Before freeze-out, i.e. for 1 x ≤ x f (x f denotes freeze out of DM), ∆ Y eq φ and d∆ dx → 0 and near freeze-out, i.e. for x x f , so ∆(x f ) = cY eq φ (x f ).One then writes Using the equilibrium distributions Y eq h 2 (x) = 0.145 Eq.( 13) takes the form 1 g * is the effective relativistic degrees of freedom.
Eq.15 can be iteratively solved for x f .Now Y eq φ ∆ for x x f .Eq.12 can then be integrated from the freeze-out epoch to finally yield . The final relic abundance is related to the comoving density Y φ (x → ∞) as Ωh 2 = 2.752 × 10 8 m φ Y φ (x → ∞) [7,52].

III. FIRST ORDER PHASE TRANSITIONS AND NUMERICAL ANALYSIS
We are interested in cosmological phase transitions in the direction of the field δ given its crucial role in the generation of the DM relic through the 3 → 2 transitions.The treelevel potential in the δ-direction looks like V0(δ cl ) = Here δ cl denotes the classical rolling field and is distinct from the degree of freedom δ .Next, the one-loop Coleman-Weinberg correction to the tree-level potential reads The sum in Eq.( 17) runs over the scalars of the theory since these are the fields that develop δ cl -dependent masses.Also, na denotes number of degrees of freedom of the ath scalar.One finds n φ = n δ = n h = n G 0 = 1 while n G + = 2.The field-dependent masses are The one-loop correction to the scalar potential induced in presence of T = 0 reads Here, the function JB( It can be approximated in the high temperature limit as where aB ≡ 16π 2 e 3 2 −2γ E .Further, the infrared effects are included using the daisy remmuation technique [53][54][55][56][57][58][59][60][61].More precisely, the Arnold-Espinosa prescription [60] in adopted in this study.That is The total thermal scalar potential as a function of δ cl and T is the sum of the individual components.That is We attempt to understand V total (δ cl , T ) more intuitively with the aid of appropriate approximations.The thermal potential can be cast as a polynomial in δ cl using the high temperature approximation for JB( m 2 T 2 ) given in Eq.( 20) and discarding terms that do not contain δ cl .One then straightforwardly derives where, λ(T ) = . The form in Eq.(III) resembles the SM formula in [62] and therefore permits an analytical treatment of the thermal evolution.We identify certain temperature thresholds in Eq.(III).An extrema at δ cl = 0 is readily identified.Moreover, this extrema is a minima (maxima) for T > T0 (T < T0) with T0 = −24µ 2 δ /(λ12 + λ22 + λ23).An inflection point is present for T = T1 is encountered.The temperature T1 can be solved from λ For T0 < T < T1, a maxima and a minima appear at δ cl = δmax(T ), δmin(T ) in addition to the minima at δ cl = 0. Thus, the T0 < T < T1 temperature band is the most interesting from the perspective of first order phase transitions on account of the coexisting minima.The critical temperature Tc for this case is the temperature at which the minima at δ cl = 0, δmin are degenerate, i.e., V total (0, Tc) = V total (δmin(Tc), Tc).A strong first order phase transition is identified by δ min (Tc) Tc ≥ 1.We put forth two benchmark points in Table I.The independent model parameters and the corresponding critical temperatures are listed therein.These benchmarks also predict relic densities within the Planck band.The shape of V total (δ cl , T ) around the crtical temperatures is shown in Fig. 4.Both BP1 and BP2 correspond to sub-GeV critical temperatures.One also inspects in Fig. 4 that δ(Tc) 1 GeV for both thereby signalling a SFOPT.
A combined analysis of the relic density and SFOPT is in order.We first fix λ13 = λ23 = 10 −5 , λ11 = 0.1, 0.5, 1 and vary the rest of the parameters as | sin θ| < 0.001, |λ12| < 6, |λ22| < 6, |µ21| < 1 GeV, 1 MeV < m φ < 1 GeV, 1 MeV < m h 2 < 2m φ .It can be noted that the small values taken for FIG.4: The extrema of the thermal potential around T = Tc for BP1 and BP2.The color coding is shown in the legends.T  sinθ and λ13 render the φ − φ − h 1(2) interactions appropriately small so as to conform with possible stringent direct detection bounds in the near future.The small sinθ also implies safety from the Higgs signal strength constraints and the direct search constraints involving h2.The full thermal potential is analysed using the publicly available tool PhaseTracer [63].We select the parameter points that lead to (a) 0.1118 < Ωh 2 < 0.1280 [5] and (b) δ min (Tc) Tc > 1.The results are shown as plots in the m φ − λ12 plane in Fig. 5.
Fig. 5 shows a O(10)-O(100)(MeV) mass range for the DM thereby corroborating the results of [20].The correlation among m φ , λ11 and λ12 can be understood from in Eq. (10).The maximum allowed value of m φ increases with increasing λ11.And this is concurred by Fig. 5 where one inspects that for the aforesaid variation of λ12, changing λ11 = 0.1 to λ11 = 1 takes one from m φ 70 MeV to m φ 190 MeV.Demanding a SFOPT restricts the parameter space further.We now discuss the GW spectrum of such a scenario and a few pertinent quantities are therefore introduced.The parameter β is defined as Here, Tn is the nucleation temperature and S3 is the Euclidean action in three dimensions.Next, we define ∆Vtot(T ) ≡ Vtot(0, T ) − Vtot(δmin(T ), T ) which measures the difference in the depths of the potential at the two vacua at a temperature T .The energy budget of the phase transition during the bubble nucleation is then given by = ∆Vtot(Tn) − T d∆Vtot(T ) dT Tn .
The α and β parameters are crucial in deciding the strength of the GW signals stemming from a SFOPT.Table II displays Tn, α, β for the chosen benchmarks.Bubble nucleation at the sub-GeV temperatures shown can be understood as an artefact of the sub-GeV DM masses involved.This is an important observation of our analysis.We also inspect in table II that α and β/H at nucleation are respectively O(10 −3 ) and O(10 3 ) for the BPs.Such ball-park values entail that the magnitude of the total GW amplitude peaks at O(10 −18 ), as is seen in Fig. 6.And the corresponding frequencies are ∼ O(10 −4 ) Hz.This is expected since the total GW amplitude is dominated by the sound wave component that has fsw ∼ O(10 −4 ) Hz. Overall, an ΩGWh 2 ∼ O(10 −18 ) with peak at f ∼ O(10 −4 ) Hz is beyond the reach of the proposed GW detectors.However, u-DECIGO [80] is expected to detect GW signals of a similar magnitude, albeit at a different frequency range.In all, a future detector with higher sensitivity in the ∼ O(0.1) mHz regime can potentially detect GW

IV. CONCLUSIONS
In this study, we look at possible GW signals originating from a SIMP.The SIMP species freezes-out via 3 → 2 number changing processes that are in turn driven by another scalar δ.The observed relic density is attained for sub-GeV DM masses.We have studied the shape of the thermal potential along δ and identified the region of the parameter space permitting a SFOPT.For a few representative benchmarks, we have then computed the GW amplitude as a function of frequency.Our analysis reveals that a futuristic detector with a sensitivity ΩGWh 2 ∼ 10 −18 in the mHz regime can detect GWs stemming from such a SIMP framework.

V. APPENDIX
The peak frequencies corresponding to bubble collisions, sound waves and turbulence are expressed as arXiv:2212.09659v1 [hep-ph] 19 Dec 2022 II.THE SUB-GEV DARK MATTER FRAMEWORK AND 3 → 2 DYNAMICS Benchmark model parameters along with the predicted relic densities and the crtical temperatures.

TABLE II :
The nucleation temperature and other GW parameters corresponding to BP1 and BP2.