LIGO-VIRGO constraints on dark matter and leptogenesis triggered by a first order phase transition at high scale

We study the possibility of constraining a scenario with high scale first order phase transition (FOPT) responsible for the cogenesis of baryon and dark matter using gravitational wave (GW) (non)-observations. While the FOPT at high scale is responsible for generating baryon asymmetry through leptogenesis and dark matter via the \textit{mass-gain} mechanism, the resulting GW spectrum falls within the ongoing LIGO-VIRGO experimental sensitivity. The dark matter is preferred to be in the non-thermal ballpark with sub-GeV masses and the criteria of successful dark matter relic rules out a large portion of the parameter space consistent with high scale FOPT and successful leptogenesis. Some part of the parameter space allowed from dark matter and leptogenesis criteria also gives rise to a large signal-to-noise ratio at ongoing experiments and hence can be disfavoured in a conservative way from the non-observation of such stochastic GW background. Future data from ongoing and planned experiments will offer a complementary and indirect probe of the remaining parameter space which is typically outside the reach of any direct experimental probe.


I. INTRODUCTION
Presence of dark matter (DM) and baryon asymmetry in the universe (BAU) has been suggested by several astrophysical and cosmological observations [1,2].These have been two longstanding problems in particle physics and cosmology given the fact that the standard model (SM) of particle physics fails to provide any explanations for the same.While several beyond standard model (BSM) proposals have been put forward to explain these observed phenomena, none of them have been experimentally tested yet.With popular production mechanisms for BAU namely, baryogenesis [3,4] or leptogenesis [5] typically remain a high scale phenomena out of direct reach of terrestrial experiments, particle DM has not been discovered yet at direct detection experiments [6].This has motivated alternative and indirect ways of probing such mechanisms behind the origin of DM and BAU.One such avenue is the detection of stochastic gravitational wave (GW) background, which has been utilised in several baryogenesis or leptogenesis scenarios [7][8][9][10][11][12][13][14][15][16][17][18][19] as well as particle DM models [20][21][22][23][24][25][26][27][28][29][30][31][32].
While most of the previous works focused on future detection of stochastic GW arising in DM and baryogenesis setups, here we consider the possibility of constraining such models with existing data from LIGO-VIRGO-KAGRA (LVK) experiments taken during their first three observing runs (O1, O2, O3).There have already been tight constraints on isotropic GW background from LVK observations [33].The same constraints have been used in the context of particle physics models with first order phase transition (FOPT) capable of generating stochastic GW in the LVK ballpark [34][35][36][37].Motivated by this, here we consider a high scale leptogenesis and DM triggered by a FOPT.Similar to the baryogenesis and leptogenesis scenarios proposed in [13,15,16], we consider a minimal setup where both DM and right handed neutrino (RHN) responsible for leptogenesis acquire masses in a FOPT by crossing the relativistic bubble walls.In our previous work [19], we focused on a low scale FOPT or low scale leptogenesis such that the resulting GW spectrum remains within the sensitivities of future experiments.In the present work, we consider a high scale version of this setup which can already be constrained by existing GW experiments like LIGO-VIRGO.
We show the parameter space consistent with successful leptogenesis which is disfavoured by LVK data.The minimal version of this model also predicts non-thermal fermion singlet DM with mass in the sub-GeV ballpark.This paper is organised as follows.In section II we briefly discuss our model followed by the details of our results in section III.We finally conclude in section IV.

II. THE MODEL
In order to show the key results, we consider a minimal setup where the type-I seesaw for light neutrino masses [38][39][40][41] is extended by a scalar doublet η with an additional Z 2 symmetry under which η and one of the RHNs are odd.We also impose a classical conformal invariance such that the mass terms arise only after a singlet scalar S acquires a non-zero vacuum expectation value (VEV) while also driving a FOPT.The relevant part of the Yukawa Lagrangian is given by where Φ 1 is the SM Higgs doublet.Thus, two of the RHNs even under Z 2 take part in type-I seesaw while the Z 2 -odd sector contributes radiatively to one of the light neutrino masses in scotogenic fashion [42].The scalar potential of the model can be written as Since we are interested in the singlet scalar induced FOPT at high scale, we denote the singlet scalar as S = (ϕ + M )/ √ 2 with M denoting the singlet scalar VEV as well as the scale of renormalisation.Since we are assuming a classical conformal invariance, the singlet scalar VEV not only decides the physical masses of RHNs and η, but also generates the scale of electroweak symmetry breaking dynamically.This constraints the parameter λ 6 < 0 to a small value for high scale FOPT.For simplicity, we consider the two Z 2 -even RHNs to be quasi-degenerate while the Z 2 -odd RHN to be much lighter, playing the role of DM.

III. RESULTS AND DISCUSSION
In order to study the details of the FOPT, we consider the tree level potential V tree mentioned above, one-loop Coleman-Weinberg potential V CW [43] along with the finitetemperature potential V th [44,45] such that the full potential is V tot = V tree +V CW +V th .While calculating the thermal potential we also include the Daisy corrections [46][47][48] which improve the perturbative expansion during the FOPT.Out of the two popularly used schemes namely, Parwani method and Arnold-Espinosa method, we use the latter.The details of the finite-temperature potential can be found in appendix A. The FOPT proceeds via tunneling, the rate for which is estimated by calculating the bounce solution.The bounce action S 3 (T ) determines the tunneling rate per unit volume defined as where and S 3 (T ) are respectively determined by the dimensional analysis and given by the classical configuration, called bounce.The bounce solution can then be obtained by following the prescription given in [49].The details of this prescription and our calculation are given in appendix B. The nucleation temperature T n of the FOPT is then calculated by comparing the tunnelling rate to the Hubble expansion rate as where t = log(ϕ/µ) with µ = M being the scale of renormalisation and the function G(t) is given by The relevant Yukawa and scalar potential couplings as a function of energy scale are calculated by solving the corresponding renormalisation group evolution (RGE) equations [19], the details of which can be found in appendix C. Finally, the temperature at which the FOPT is completed, known as the percolation temperature T p is calculated by following the prescription given in [52,53].According to this prescription, T p is obtained from the probability of finding a point which is still in the false vacuum, given by P(T ) = e −I(T ) . Here, The percolation temperature is then calculated by using I(T p ) = 0.34 [52] which implies that at least 34% of the comoving volume is occupied by the true vacuum.GeV, In Fig. 1, we show the parameter space in terms of physical masses of the heavier quasidegenerate RHNs denoted by M N , neutral real component of scalar doublet η and scalar singlet denoted by M N , M η , M S respectively.The parameter space is consistent with a high scale FOPT with the colour code showing the nucleation temperature T n .We also identify the points with M N < T n where one has to consider thermal production of heavy neutrinos into account while estimating the lepton asymmetry produced.
In order to calculate the stochastic GW spectrum as a result of the FOPT, we take all relevant contributions into account from the bubble collisions [54][55][56][57][58], the sound wave of the plasma [59][60][61][62] and the turbulence of the plasma [63][64][65][66][67][68].The two important quantities namely, the duration of the phase transition and the latent heat released are calculated and and α * respectively.The action is evaluated numerically by fitting our potential using the procedure laid out in [69] and utilised in our earlier work [19].The bubble wall velocity v w is estimated by first calculating the Jouguet 70,71] which then leads to the bubble wall velocity v w as Here, ρ rad = g * π 2 T 4 /30 denotes the radiation energy density and ∆V tot denotes the energy difference between the true and the false vacua, given by The latent heat released during the phase transition can be estimated as with which is also related to the change in the trace of the energy-momentum tensor across the bubble wall [74,75].
The stochastic GW spectrum, considering all the sources, can be written as [76,77] where Ω ϕ , Ω sw , Ω turb correspond to the individual contributions from bubble collisions, sound wave of the plasma and turbulence in the plasma respectively.In general, each of these contributions has a peak type feature with peak frequency f peak .The spectrum can be parametrised as where the pre-factor R takes into account the red-shift of the GW energy density, S(f /f peak ) parametrises the shape of the spectrum and ∆(v w ) is the normalization factor which depends on the bubble wall velocity v w .κ is the efficiency parameter which denotes the fraction of latent heat driving that particular source of GW in a first order phase transition.The numerical values of the parameters p, q depend upon the source.For bubble collision as the source, the spectrum has been re-estimated in several recent works [52,53,[78][79][80] and can be written as [77,78,81] T n 100 GeV 0.64 2π The efficiency factor κ ϕ for bubble collision is given by [63] The GW spectrum generated from the sound wave in the plasma has been studied through large hydrodynamical simulations [62] which has also been updated in several recent works [74,81,82].The corresponding spectrum can be written as [81] The corresponding peak frequency is given by f sw peak = 8.9 × 10 −6 Hz g * 100 The efficiency factor for sound waves, applicable for relativistic bubble wall velocity v w ∼ 1 in our model, is [71] Here, 1+2τswH * is a suppression factor which depends on the lifetime of sound wave τ sw [82] and it can be written as τ sw ∼ R * / Ūf with mean bubble separation, R * = (8π) 1/3 v w β and rms fluid velocity, Ūf = 3κ sw α * /4(1 + α * ); z p ∼ 10.Finally, the spectrum generated by the turbulence in the plasma is given by [76,77,81] with the peak frequency being [76] f turb peak = 2.7 × 10 −5 Hz g * 100 The efficiency factor for turbulence is κ turb ≃ 0.1κ sw [76] and the inverse Hubble time at the epoch of GW production, redshifted to today is The total contribution to the stochastic GW spectrum is shown in Fig. 2 for a few benchmark points shown in table I. Since we are focusing on high scale FOPT, we are showing the relevant sensitivities of future experiments like DECIGO [83], ET [84] and ongoing LIGO-VIRGO (HVO3) [33,85] as shaded regions of different colours.
We then implement the leptogenesis via relativistic bubble wall or mass-gain mechanism [13,15,16,19] for the high scale FOPT scenario.The right handed neutrinos N 1,2,3 and scalar doublet η acquire masses after entering the bubble formed due to the FOPT induced by the singlet scalar discussed above.The quasi-degenerate heavier RHNs namely, N 2,3 then decay into leptons and Higgs to generate the lepton asymmetry while N 1 , being the lightest Z 2 -odd particle emerges as the DM candidate.The CP asymmetry parameter corresponding to the CP violating decay of RHN N i (summing over all lepton flavours) is given by [86] where, only self-energy correction is considered, which dominate for quasi-degenerate RHNs.
Here M i denotes the mass of RHN N i .Since we consider quasi-degenerate N 2,3 , we denote their mass as M N hereafter.The relativistic nature of the bubble walls arising out of the supercooled phase transition ensures the penetration of RHNs to maintain a large abundance inside the bubble, same as the equilibrium abundance without a Boltzmann suppression.For this, we need to ensure that the Lorentz boost of the bubble wall should be more than the Lorentz factor of the particle in the plasma frame where T n is the nucleation temperature.However, due to T n < M N , RHNs can not be thermally produced but yet having a large comoving abundance Y N = n N /s inside the bubble without any Boltzmann suppression, with n N , s being equilibrium number density of N and entropy density of the universe respectively.The comoving abundance of N inside the bubble is then evaluated as where g N and g * are the degrees of freedom of RHN N and the total relativistic degrees of freedom of the universe, respectively.For the parameter space with M S > 2M N , we also take the additional contribution to Y N from singlet scalar S decay.The final baryonic asymmetry can then be approximated as where ϵ N ≡ ϵ 2,3 ≃ sin(2δ)/(16π) [86,87] is the CP-asymmetry and δ is the relative CP phase between the quasi-degenerate RHNs, κ sph = 8/23 is the sphaleron conversion factor for our model, and T RH is the reheating temperature after the FOPT.T RH is defined as where T inf can be found by comparing radiation energy density to the energy released from the FOPT or equivalently ∆V tot namely, We also check the feasibility of RHN decay into L, Φ 1 at the reheating temperature by considering the thermal masses of daughter particles.Since thermal masses of L, Φ 1 scale as T , the requirement of Table II shows the required CP asymmetry along with other relevant parameters involved in calculating baryon asymmetry for the benchmark points discussed before.
The validity of baryon asymmetry estimate given by Eq. ( 27) also depends upon the assumption that the washout processes are absent or negligible.In particular, the dominant wash-out coming from inverse decay should be suppressed to validate Y B estimated in Eq. (27).While this is true due to M N > T n , presence of L, Φ 1 with energy more than the mean thermal energy can induce such washout.As the RHNs N 2,3 enter the bubble, they receive a boost in the plasma frame with energy γ N M N , where γ N ∼ M N /T n .Consequently, the decay products L, Φ 1 and L c , Φ * 1 also appear to be boosted, with the energy of the order M 2

N
Tn .These boosted particles can interact with SM particles in the plasma of energy T n and go through cascade scattering with particles in the plasma while redistributing their energies.
At the same time, these energetic particles can also lead to washout asymmetry through on-shell processes such as L Φ 1 → L c Φ * 1 via N 2,3 .With the simplest assumption of energy distribution, the energy of boosted particles in n-th step cascade scattering is 2 n Tn .Following the procedure of [13,15], the rate of the washout in the above process can be estimated is total decay width of RHN.Similarly, the thermalization rate can be estimated [13,15] as Γ th ≈ ln( where g EW = 46, α W = 1/137.We check that for a typical benchmark point considered in our analysis, Γ th ≫ Γ wash such that all the boosted particles get thermalised quickly via cascade interactions with the SM bath, before inducing any washout effects via inverse decay.Therefore, for our scenario, washout processes are negligible and the baryon asymmetry estimate given by Eq. ( 27) remains valid.
Adopting a conservative approach to apply the GW constraints on the parameter space of our model consistent with a high scale FOPT and leptogenesis, we define the signal-to-noise FIG.2: GW spectrum for different benchmark points of our model shown in table I.
with τ being the observation time for a particular detector, which we consider to be 1 yr.To register a detection SNR needs to be more than a threshold value ρ > ρ th .In this work, we consider ρ th = 1, such that the parameter space giving rise to SNR ρ > 1 at LIGO-VIRGO latest run namely O3 is disfavoured due to non-detection of such stochastic GW, as per our conservative selection criteria.
In Fig. 3, we show the parameter space of our model in terms of physical masses of heavier quasi-degenerate RHNs, singlet and doublet scalar masses in a way similar to the ones in Fig. 1  be calculated in a way similar to thermal leptogenesis while solving the relevant Boltzmann equations explicitly.We also calculate the peak amplitudes of GW for the parameter space consistent with FOPT requirements and show them against LIGO-VIRGO sensitivities on the left panel of Fig. 4. We show the sensitivity curves of run 2 (HVO2) [90] as well as run 3 (HVO3) [33] for comparisons.While the peak frequencies of the points with SNR > 1 lie outside the sensitivity curves, they can still be detected if we consider the entire GW spectrum.On the right panel of Fig. 4, we show the heavy quasi-degenerate RHN massM N versus SNR with colour code indicating the nucleation temperature.Clearly, some of the points with T n ∈ (10 6 − 10 7 ) GeV get disfavoured.
Dark matter or the Z 2 -odd RHN can, in principle, be produced either thermally or nonthermally.For TeV scale FOPT in this model, studied earlier [19], it was found that thermal DM parameter space is almost entirely ruled out due to stringent direct detection bounds.This is precisely due to the constraint on quartic coupling between singlet scalar and the SM Higgs from the requirement of successful electroweak symmetry breaking.This coupling, namely λ 6 in our notation then decides Higgs portal coupling through which fermion singlet thermal DM typically annihilate.While DM has Yukawa coupling with leptons and η, Yukawa portal annihilations typically remain sub-dominant especially due to large η mass required by FOPT criteria and often require large Yukawa couplings to satisfy relic.Such Yukawa coupling faces tight constraints from neutrino mass as well as charged lepton flavour violation [91].If we go to high scale where some part of the parameter space gets constrained by LIGO-VIRGO data, thermal DM is likely to hit the unitarity limit [92].While DM mass need not be same as the scale of FOPT, making it much lighter suppresses the annihilation rates due to heavy mediators in the form of η as well as singlet scalar.In addition to mediator suppression, DM coupling with singlet scalar also gets smaller for lighter DM masses.Therefore, thermal DM gets overproduced in such a scenario due to insufficient annihilation rates.While entropy dilution at the end of FOPT could lower this thermal abundance, we find such entropy dilution to be negligible as discussed above, in the context of leptogenesis.Therefore, we consider non-thermal DM which remains feasible for some part of the parameter space under study.
Such non-thermal or feebly interacting massive particle (FIMP) type DM can be produced dominantly from singlet scalar (S) decay inside the bubble.The corresponding Boltzmann equations for comoving densities of DM and singlet scalar S can be written as with z = M S /T and assuming the relativistic dof to be constant, which is valid at temperatures above electroweak scale.In the second equation, Γ sDM , Γ sh , Γ sN 2 , Γ sN 3 denote the corresponding decay width of singlet scalar into a pair of DM, SM Higgs, N 2 , N 3 respectively.
Similar to N 2,3 , the singlet scalar with mass M S > T n is also out-of-equilibrium inside the bubble having a large initial abundance, given by Since singlet scalar coupling with the SM Higgs is very small, one requires a dominant decay channel of S such that it decays at least before BBN without overproducing DM.Due to M S > T n applicable to the entire parameter space, dilution of singlet scalar abundance via annihilation is also suppressed.In the minimal model we are studying, DM overproduction from singlet scalar decay can be avoided only when M S > 2M N 2,3 ∼ 2M N keeping the branching ratio BR(S → DM DM) small.We have chosen our benchmark points in table I, II such that this condition is satisfied.The corresponding DM masses, consistent with correct relic abundance, are shown in the last column of table II.On the left panel of Fig.
The thermal contributions to the effective potential are given by where n B i and n F i denote the dof of the bosonic and fermionic particles, respectively.
The J B and J F functions are defined by following functions: (A6) In thermal potential, we also consider the contribution from daisy diagrams [46][47][48].
Considering Arnold-Espinosa method [48], the thermal potential with the daisy correction can be written as V T (ϕ, T ) = V th + V daisy (ϕ, T ), (A7) Denoting m 2 i (ϕ, T ) = m 2 i (ϕ) + Π i (T ), the relevant thermal masses can be written as [93] Π η (T ) = ( g where, the coefficients A and B depend on the potential.In Euclidean space or imaginary time formalism, the coefficient B is S 3 /T with S 3 being the Euclidean action at finite temperature in three dimensions.As described in [49], the dimensional estimation of preexponential factor A is T 4 S Here, we used a fitting method to calculate the action.We fitted the effective potential to a generic potential for which the action calculations are done in semi-analytical way [69].The generic quartic and logarithmic potential is We show that the actual potential and the fitted generic potential overlap quite well, depicted in left panel of Fig. 6 with H * ≡ H(T = T n ).The strength of the FOPT is conventionally decided by the order parameter ϕ(T c )/T c ≡ v c /T c with ϕ(T c ) being the singlet scalar VEV at critical temperature T = T c at which the two minima of the potential are degenerate.Larger is the order parameter v c /T c > 1, stronger is the FOPT.The bounce calculation is done by rewriting the zero temperature one-loop effective potential as[50,51]

FIG. 1 :
FIG. 1: Parameter space in M η versus M N plane (left panel) and M η versus M S plane (right panel) consistent with a high scale FOPT v c /T c > 1 with the colour code showing the corresponding nucleation temperature T n .In this scan, the scale of phase transition (M) is varied from to 10 8

FIG. 3 :
FIG. 3: Parameter space in M η versus M N plane (left panel) and M η versus M S plane (right panel) consistent with a high scale FOPT showing the SNR with respect to LIGO-VIRGO O3.In this scan, the scale of phase transition (M) is varied from 10 4 to 10 8 GeV, λ 7 (0) is varied from 1 to 3,

10 7FIG. 4 :
FIG. 4: Left panel: GW spectrum peak vs its corresponding Peak frequency with the colour code showing the SNR.Right panel: Heavy RHN mass versus SNR with colour code indicating the corresponding nucleation temperature.

5 , 7 FIG. 5 :
FIG.5: Left panel: The parameter space in M N vs M S plane where ⋄ shaped points are consistent with dark matter relic while △ shaped points are not, and ⋆ shaped points satisfy dark matter relic as well as non-thermal leptogenesis criteria.Right panel: The parameter space in M S − M DM plane which satisfies DM relic criteria.In both the plots, the points labelled as M N /T n > 1 are consistent with non-thermal leptogenesis criteria via mass-gain mechanism.

TABLE I :
Benchmark parameters and other details involved in the GW spectrum calculation of the model.