Gravitational Waves, Bubble Profile, and Baryon Asymmetry in the Complex 2HDM

This study explores the generation of the observed baryon asymmetry of the Universe within the complex Two Higgs Doublet Model (C2HDM) while considering theoretical and current experimental constraints. In our investigation, we analyze critical elements of the Higgs potential to understand the phase transition pattern. Specifically, we examine the formation of the barrier and the uplifting of the true vacuum state, which play crucial roles in facilitating a strong first-order phase transition. Furthermore, we explore the potential gravitational wave signals associated with this phase transition pattern and investigate the parameter space points that can be probed with LISA. Finally, we compare the impact of different approaches to describing the bubble profile on the calculation of the baryon asymmetry. We contrast the typically used kink profile approximation against the explicit solution of the tunneling profile. We find that a non-negligible range of the C2HDM parameter space results in significant discrepancies in the baryon asymmetry estimation between these two approaches. Through an examination of the parameter space, we identify a benchmark point that satisfies the observed baryon asymmetry.


I. INTRODUCTION
Understanding the origin of the matter-antimatter asymmetry of the Universe, known as the baryon asymmetry of the Universe (BAU), is a fundamental question in particle physics and cosmology.The asymmetry between baryons and antibaryons in the early Universe can be quantitatively evidenced through the baryon-to-entropy ratio measurement n B /s ≃ 8.6 × 10 −11 [1], exceeding the expected value for a symmetric scenario by several orders of magnitude.As a consequence, the majority of antibaryons underwent annihilation during the thermal history, leaving behind a significant density of baryons in the present Universe.
The essential ingredients required for generating this baryon asymmetry are theoretically well understood and encapsulated by the three Sakharov conditions [2].These conditions demand the violation of baryon number, the presence of C and CP violation, and a departure from thermal equilibrium.While the Standard Model (SM) satisfies the requirements for baryon number violation and C violation, it falls short in providing a sufficiently robust source of CP violation.Additionally, the observed Higgs mass of m h = 125 GeV precludes the necessary out-of-equilibrium conditions through a strong first-order phase transition [3,4].Thus, the quest for baryogenesis requires physics beyond the SM [5][6][7][8].Among possible extensions, the complex Two-Higgs Doublet Model (C2HDM) can potentially provide both of the missing ingredients: strong first-order electroweak phase transition and additional sources of CP-violation [9][10][11].In this work, we explore the phase transition pattern and the feasibility of generating the observed baryon asymmetry within the context of the C2HDM.Central to our investigation is the shape of the Higgs potential, which plays a crucial role in determining the nature of the phase transition.We focus on the formation of the barrier and the upliftment of the true vacuum state, as these factors are instrumental in driving the phase transition from a smooth crossover to a strong first-order transition.
Our analysis builds upon previous studies for other new physics extensions [12][13][14][15], where it was observed that the intensity of the phase transition is closely linked to the elevation of the true vacuum relative to the symmetric one at zero temperature.The prevalence of one-loop effects over thermal corrections, particularly when ξ c > 1, enhances the strength of the phase transition [15].However, it should be noted that if the one-loop correction is too large, the universe may become trapped in the electroweak symmetric vacuum, resulting in an incomplete phase transition [15,16].Consequently, as we will show, a significant portion of parameter points with large ξ c values become unphysical in this scenario.
The first-order phase transition in the early Universe can generate stochastic gravitational waves (GW) whose characteristic peak frequency is associated with the phase transition temperature.After redshifting to the present time, the GW spectrum would have a peak frequency at the mHz range for the phase transition at the electroweak scale [17,18].This presents an exciting prospect to probe electroweak phase transition (EWPT) at LISA [19], designed to be sensitive to mHz frequency signals.Hence, we also investigate the parameter space points in C2HDM that can be probed using LISA.
Through an extensive exploration of the parameter space, we note that the C2HDM can describe the observed baryon asymmetry, although only for a limited set of parameter space points.In this regard, we compare two different approaches to describe the bubble profile, a key ingredient in the BAU estimation.The commonly adopted kink profile parameterization [20][21][22][23][24] and the explicit solution for the tunneling equation are examined to assess their impact on the resulting baryon asymmetry.Our analysis reveals relevant deviations between the BAU calculation between these two approaches.While the majority of parameter space points yield similar results using both methods, a notable fraction exhibits significant differences, sometimes varying by several orders of magnitude.To understand these discrepancies, we scrutinize the behavior of the source term in front of the bubble wall, which sheds light on the distinct asymmetry values obtained from the two profile assumptions.
The paper is organized as follows.In Section II, we provide a brief overview of the complex Two Higgs Double Model.Section III discusses the one-loop finite temperature effective potential.It is followed by a discussion on electroweak phase transition and GW signals in Section IV.In Section V, we study how the shape of the Higgs potential will affect the EWPT, focusing on the barrier formation and the vacuum upliftment.In Section VI, we present the details for the baryon asymmetry calculation.The results of the BAU are presented in Section VII, where we also contrast the results depending on the bubble profile estimation.Finally, we summarize in Section VIII.Details of the parameterization for the C2HDM scan are presented in Section A.

II. COMPLEX TWO HIGGS DOUBLET MODEL
The two Higgs doublet model (2HDM) lays out a compelling extension of the SM in line with current experimental constraints [25].This work considers CP-violating 2HDM with a softly broken Z 2 symmetry.Within this framework, the tree-level potential is given by where the mass term m 2 12 and quartic coupling λ 5 are complex and all other mass terms and quartic couplings are taken to be real.However, one of the phases in m 2  12 and λ 5 can be removed by a phase redefinition of Φ 2 .In this work, we always keep m 2 12 real and λ 5 will be complex at zero temperature.Hence, overall in such setup, there is only one independent physical CP violation phase.To preclude dangerous tree-level Flavor Changing Neutral Currents (FCNC) [26,27], we impose a Z 2 symmetry softly broken by the m 2  12 term, under which Φ 1 → Φ 1 and Φ 2 → −Φ 2 .Following electroweak symmetry breaking, the neutral components of Φ 1 and Φ 2 develop non-zero vacuum expectation values (VEVs).Expanding around the VEVs ω i , the scalar doublets Φ i can be written as where at zero temperature VEVs Whereas an additional source of CP-violation should decrease at zero temperature (ω θ | T =0 → 0) to comply with the stringent electric dipole moment (EDM) constraints [28], the dynamical generation of CP-violation at high temperatures offers a potential avenue for a CP-violating mechanism crucial to the success of Electroweak baryogenesis.To account for a more comprehensive scenario, we also incorporate a possible charge-breaking at high temperature, ω CB .Since a non-zero charge-breaking VEV at zero temperature would lead to massive photons, we impose v CB = 0.
The scalar sector in the CP-violating 2HDM has five physical mass eigenstates: three CP-mixed neutral scalars H i and one charged scalar pair H ± .The correspondence between mass eigenstates and gauge eigenstates is established by the mixing angle β in CP-odd and charged sectors and another three angles α, α b and α c mixing the CP-odd and CP-even scalars: The mixing angle β is defined as We also define c x ≡ cos x and s x ≡ sin x.
At zero temperature, the physical parameters in the scalar sector include the VEVs , the masses of the scalar eigenstates (m H i and m H ± ), the mixing angles (α, α b , and α c ), and m 2 12 .Note that, as we mentioned earlier, there is only one physical CP violation phase, i.e., only one of θ, α b , and α c is independent.In this work, we keep α c as an independent input while calculating θ and α b from other parameters.
Hence, we choose the input parameters to be which match the 9 real parameters in the potential Eq. ( 1).Here, m H ↑ and m H ↓ represent the masses of heavier and lighter beyond the Standard Model (BSM) neutral scalars, respectively.
The detailed mapping between the parameters in Eq. ( 5) and those in Eq. ( 1) can be found in Section A. This parameterization for the CP-violating 2HDM is similar to the scan performed for CP-conserving 2HDM in our earlier works [15,29] in the sense that it provides the scans over all physical BSM scalar masses (m H ↑ , m H ↓ and m H ± ) and CP-violating angle (α c ).The phase transition pattern in 2HDM, to a large extent, depends on the masses of additional scalars and corresponding mass splittings [15].Hence, the numerical scan performed over three scalar masses is more suitable than one of the scalar masses written as a function of other scan variables.
Within the Yukawa sector, there are four distinct Z 2 charge assignments that effectively preclude tree-level FCNC.In this study, we focus on two specific scenarios: type-I and type-II.In the type-I scenario, all fermions exclusively couple with Φ 2 , while in the type-II scenario, only up quarks couple with Φ 2 , with down quarks and charged leptons coupling with Φ 1 .To thoroughly explore these possibilities, we conduct a random uniform scan, encompassing both type-I and type-II configurations, over the parameter space region tan β ∈ (0.8, 25) , m 2 12 ∈ (10 −3 , 5 × 10 5 ) GeV 2 , m H ↑/↓ ∈ (30, 1500) GeV , We performe the parameter space scan by implementing the parametrization detailed in Section A in ScannerS [30,31].Using ScannerS, we impose constraints from perturbative unitarity [32][33][34], boundedness from below [35], vacuum stability [36,37], electroweak precision, and flavor constraints.EDM constraints are also imposed using the stringent limits from the ACME collaboration [28].Furthermore, constraints from the 125 GeV Higgs boson measurements and additional scalar searches are carried out using HiggsBounds and HiggsSignals [38][39][40].

III. ONE-LOOP FINITE TEMPERATURE EFFECTIVE POTENTIAL
We use loop-corrected finite temperature effective potential to determine the dynamics of electroweak symmetry breaking in the early Universe.Along with the tree-level potential V 0 from Eq. (1), we also include the Coleman-Weinberg potential V CW and counterterms V CT that encode one-loop corrections at zero temperature, and finite-temperature corrections V T .
The effective potential is given by The Coleman-Weinberg potential in the Landau gauge can be written, using MS renormalization prescription as [41] where the index i runs over all particles in the thermal bath with field-dependent mass m i (Φ 1 , Φ 2 ), including Higgs bosons, massive gauge bosons, Goldstone bosons, longitudinal photon, and fermions.The parameter n i represents the number of degrees of freedom for each particle, with n i > 0 for bosons and n i < 0 for fermions.In the MS renormalization procedure, the coefficient c i takes the value of 5/6 for gauge bosons and 3/2 otherwise.
Moreover, we set the renormalization scale µ to the zero-temperature VEV, µ = v(T = 0) ≈ 246 GeV. 1   The one-loop effects of the Coleman-Weinberg potential result in shifts of the mixing angle and scalar masses from their tree-level values.To perform a consistent parameter scan, we adopt an on-shell renormalization scheme, which enforces the parameters to match their tree-level values [24,44,45], by proper counterterms determined according to where ϕ i (i = 1, ..., 8) represents scalar components from the Φ 1 and Φ 2 doublets, ω denotes the ω i values, and ω tree characterizes the minimum of the tree-level potential for the fields in Φ 1 and Φ 2 .The first and second derivatives of V CW are consistently defined with an analytical expression in Ref. [44]. 2 The first renormalization condition, given by Eq. ( 9), ensures that the minimum of the effective potential is not shifted from tree-level minimum, and the second condition, shown in Eq. ( 10), guarantees that mixing angles and scalar masses remain the same as their tree-level values.
The one-loop thermal correction V T in Eq. ( 7) is given by [48] where the sum extends over fermions f and bosons.The bosonic sector can be further divided into two categories: the transverse modes of gauge bosons, represented by V T = W T , Z T , and the longitudinal modes of gauge bosons and scalars, denoted by The resummation of the n = 0 Matsubara modes of the longitudinal components V L leads to thermal corrections in their masses [49,50].The second line in Eq. ( 11) corresponds 1 A renormalization group improved calculation can be taken into account for a further refined estimation [42].For the renormalization scale µ 2 dependence of effective potential at finite temperature, we refer to Ref. [43]. 2 The first and second derivatives of the effective potential exhibit IR divergence at zero temperature as a result of the contribution from Goldstone bosons [46,47].To address this, we use the analytical expression provided in Ref. [44] to evaluate derivatives of effective potential and compute the counterterms.
to the Daisy contributions, where m V L represents the thermal Debye mass calculated using the Arnold-Espinosa scheme [45,48].The formulas for Debye masses are provided in Appendix C. Finally, the thermal functions for fermions (J + ) and bosons (J − ) are given by Whereas the effective potential in the electroweak phase transition is subject to theoretical uncertainties stemming from gauge parameter choices [42,[51][52][53][54][55][56][57], Nielsen identities offer a way to construct gauge-independent probes [58].These identities ensure that the gauge dependence cancels out at the extrema of the potential where ξ is the gauge fixing parameter.Inspired by the gauge independence guaranteed by Nielsen identities, we employ two distinct methods for phenomenological analyses.The first approach involves calculating the finite-temperature effective potential and performing a numerical scan.The second approach focuses on determining the gauge-invariant vacuum upliftment at T = 0.In Section V, we highlight that the upliftment of the true vacuum relative to the symmetric vacuum at zero temperature serves as an effective probe of the phase transition's strength.While the first method carries uncertainties associated with gauge parameter choices, the latter approach is gauge invariant, as assured by Nielsen identities [12,51].It is worth noting that we introduce additional counterterms at one-loop order to preserve the positions of the electroweak vacuum and masses.The agreement between our numerical scan and the profile derived from the vacuum upliftment serves to confirm the reliability of the numerical scan despite its inherent uncertainties.3

IV. ELECTROWEAK PHASE TRANSITION AND GRAVITATIONAL WAVES
The finite temperature effective potential dictates the phase-transition pattern.The two Higgs doublet model displays both single and multi-step phase transitions.The first-order phase transition occurs through tunneling from false to true vacua.It results in bubbles of the broken phase that pop up and expand in the surrounding region of the symmetric phase, transitioning from the false vacuum to the true vacuum.The tunneling probability is given by [60,61] where S 3 represents the three-dimensional Euclidean action associated with the critical bubble formation Here, the scalar field ϕ corresponds to the critical bubble profile, which is determined by solving the following differential equation We utilize the publicly available code CosmoTransitions [62] to solve the differential equation and compute the Euclidean action S 3 .
The first-order phase transition is considered to be completed around the nucleation temperature T n , which corresponds to the point where one bubble nucleates per unit horizon volume [63].
This condition ensures that the bubbles percolate even in the inflating Universe.For the electroweak phase transition, with a nucleation temperature of approximately T n ≈ 100 GeV, this condition can be approximated as [50] S 3 (T ) To preserve the baryon asymmetry generated through electroweak baryogenesis, it is crucial to suppress the sphaleron process inside the bubble.This requires the electroweak symmetry breaking to undergo a strong first-order phase transition [50] where ) is the Higgs VEV at the critical temperature T c .This critical temperature corresponds to the point where the broken and unbroken vacua of the electroweak symmetry are degenerate.The approximate inequality in Eq. ( 19) indicates the theoretical uncertainty in this condition [51].
The production of stochastic gravitational waves is a significant consequence of a firstorder phase transition.These GW originate from three main sources: the collision of vacuum bubbles, fluid motion resembling sound waves in the plasma, and turbulent motion within the plasma.Each source contributes to the GW spectrum, which can be described by numerical functions dependent on two parameters that capture the dynamics of the phase transition at the nucleation temperature T n [17,[64][65][66]. 4The first parameter is α, defined as the ratio of the latent heat released during the phase transition (ϵ) to the energy density of the vacuum radiation (ρ rad ), i.e., α ≡ ϵ/ρ rad .The latent heat and the vacuum radiation energy density are expressed as where ∆ represents the difference between the true and false vacua, and g ⋆ the number of relativistic degrees of freedom in the plasma.The second important parameter is β/H n , which characterizes the inverse time duration of the phase transition.This quantity is defined as where H n is the Hubble constant at the nucleation temperature T n .Detectable GW signals are typically associated with a slow phase transition (small β/H n ) and a large latent heat release (large α).
Finally, to assess the detectability of GW signal, we employ the signal-to-noise ratio where Ω Sens represents the sensitivity curve of the considered GW detector [67] and T corresponds to the mission duration.For our analysis, we adopt the LISA gravitational wave detector as a benchmark, with T = 5 years and a detection threshold of SNR = 10 [64].
The existence of percolation temperature T p , where 29% of space is covered by bubbles, guarantees the completion of phase transition [68].Supercooling emerges when the nucleation temperature is substantially lower than the critical temperature, leading to a pronounced α ≫ 1 [68,69].In our analysis, the overwhelming majority of data points in our scan, specifically 99.99%, display α < 1.Hence, we assume there is no supercooling and the percolation temperature can be approximated to the nucleation temperature T p ≃ T n .

V. BARRIER FORMATION AND VACUUM UPLIFTMENT
Introducing a second Higgs doublet to the SM Higgs sector can alter the behavior of electroweak symmetry breaking from a smooth crossover to a strong first-order phase transition.In Ref. [15], the authors studied the key ingredients that trigger this transmutation in the EWPT by focusing on the barrier formation and upliftment of the true vacuum in the context of CP-conserving 2HDM [12][13][14][15].In this model, the barrier is driven primarily by one-loop corrections and ξ c can be correlated with ∆F 0 /|F SM 0 |, a gauge independent parameter calculated at zero temperature.The ∆F 0 /|F SM 0 | is defined as where F 0 is the zero-temperature vacuum energy density of the 2HDM defined as with It is interesting to examine whether the phase transition features of the CP-conserving 2HDM prevail in the CP-violating 2HDM.In Fig. 1 (left panel), we show that the fraction of one-loop contribution to the barrier height is correlated with the zero temperature vacuum upliftment measure ∆F 0 /|F SM 0 |.We observe that the larger the one-loop correction, the higher the value of vacuum upliftment.In particular, this correlation can be seen for ξ c ≳ 1.
As one-loop effects are the dominant contributions, we can use ∆F 0 /|F SM 0 | to shed light on the properties of the EWPT.We can approximately propose ∆F 0 /|F SM 0 | ≳ 0.2 as minimal condition for strongly first-order EWPT in the CP-violating 2HDM.
In the scenario where the vacuum upliftment measure ∆F 0 /|F SM 0 | is extremely large, the tunneling from false vacuum to true vacuum becomes challenging, translating into Eq.( 18) having no solution [15,16].Thus, the Universe is trapped in a high energetic electroweak symmetric vacuum, yielding a nonphysical vacuum.In the Fig. 1 (right panel), we denote these points with the orange color.Most of the parameter points with vacuum upliftment measure ∆F 0 /|F SM 0 | ≳ 0.87 exhibit a vacuum trapped scenario.The above constraint excludes the bulk of ξ c > 2 points, which would otherwise serve as promising candidates for successful electroweak baryogenesis.
In Fig. 2 (upper-panel), we show the scanned points in the (∆m H ↑ , ∆m H ↓ ) plane colorcoded by ξ c , where m H ↑ (m H ↓ ) represents the mass of the heaviest (lightest) BSM neutral scalar, and ∆m The gray points in the background pass all the theoretical and current experimental constraints.The black points also satisfy the first-order phase transition condition with 0 < ξ c < 1.The preference to the region with The samples with ξ c < 1 have been excluded from the middle and bottom panels.In the middle panel, orange color represents parameter points trapped in false vacuum and phase transition is incomplete.The parameter space scan is implemented using ScannerS [31], where we impose the constraints from perturbative unitarity, boundedness from below, vacuum stability, electroweak precision, flavor constraints, and EDM limits.HiggsBounds and HiggsSignals are used to incorporate the searches for additional scalars as well as the 125 GeV Higgs boson measurements [38,39].
coded by ξ n .The points marked in orange correspond to locations where vacuum trapping In Fig. 2 (lower-panel), we show parameter points that can be probed by LISA in the (∆m H ↑ , ∆m H ↓ ) plane.The color coding in this case represents the logarithm (base 10) of the signal-to-noise ratio.We focus on points above the SNR threshold, SNR > 10.Among these points, we highlight the benchmark point BP4 in Tab.I, which serves as an example of a parameter point that can be probed by LISA.For the parameter points with ξ n > 1, 6% of Type-I points show a detectable GW signal by LISA, whereas it is around 2.5% for Type-II.These differences between type-I and type-II scenarios are driven by constraints from flavor physics [15].More concretely, constraints from B-meson decays impose a lower bound on the charged scalar mass requiring m H ± ≳ 580 GeV in the type-II 2HDM.In Fig. 3, we show the correlation between SNR and zero-temperature vacuum upliftment measure ∆F 0 /F SM 0 for the Type I parameter points with ξ n > 1.The bulk of parameter points that exhibit strong GW signals are associated with large ∆F 0 /F SM 0 measure.In most cases, parameter points with ∆F 0 /F SM 0 < 0.4 do not show a promissable GW at LISA.

A. Estimation of bubble wall profile
The bubble profile in the radial coordinate can be obtained by solving the tunneling equation Eq. ( 16).The baryon asymmetry calculation is performed in the bubble wall coordinate system z, where z = 0 denotes the bubble wall.We obtain the position of the bubble wall in the radial coordinate system r 0 , where the energy density obtains the maximum value.
The energy density U E is given by and the bubble wall coordinate z can be defined as A key ingredient for the baryogenesis is the complex mass of quarks and leptons, which To illustrate these concepts, we present a graphical representation in Fig. 4. The left panel displays the energy density U E (r) as a function of the radial distance r, specifically for benchmark point BP1 as defined in Tab.I.The position of the bubble wall is identified as the barrier of the tunneling profile.On the right panel, we show the dynamic variation of the CP-violating angle of the top quark with respect to temperature in the broken phase for the BP1.As thermal effects come into play, additional CP violation is induced at higher temperatures.This effect becomes prominent, whereas at zero temperature, it is roughly seven orders of magnitude smaller.The oscillatory behavior observed between temperatures of 20 GeV and 35 GeV arises due to thermal contributions that lead to a change in sign of the CP angle, which we represent in terms of the absolute value |θ t |.
In the literature, it is a customary practice to parameterize the tunneling profile θ i (z) by kink profile [20-24] where θ i brk (θ i sym ) is the phase at the broken (symmetric) minimum.The thickness of the wall L W is given by L W = v n / √ 8V b [21] with v n representing the VEV at EWPT and V b the height of the barrier that separates the two minima (at the nucleation temperature T n ).
Remarkably, this parameterization displays a tunneling profile that is symmetric with respect to the bubble wall.In Section VII, we compare the estimation of the baryon asymmetry of the universe using two different methods: the kink profile and the explicit solution from the bubble profile.In the latter case, the bubble profile is obtained by directly solving the tunneling equation.By comparing the results obtained from these two approaches, we can evaluate the consistency and reliability of the BAU estimation.

B. Semi-classical force method
The baryon asymmetry in the Universe can be estimated using the semi-classical force method.This framework utilizes the existence of a fermion with varying complex mass as it passes through the bubble wall.The particle interaction with the bubble wall can be formalized using the WKB approximation [21,22,72] or the closed-time-path formalism of thermal field theory [73][74][75][76], where the force acting on the particle is given by The radial coordinate denotes the perpendicular distance from the wall in the rest frame of the wall, where the positive direction of z points towards the symmetric phase.E 0 is the conserved wall frame energy of the quasi-particle, E 2 0z = E 2 0 − p 2 ∥ and (..) ′ denotes the derivative with respect to the z coordinate.The first term in Eq. ( 29) conserves CP, whereas the second and third terms depend on the spin and nature of the particle, with the upper sign solution corresponding to the particle and the lower sign to the antiparticle.Thus, the presence of a non-zero value for θ ′ generally indicates the appearance of CP violation [21,22,74].Assuming that the kinetic momentum is conserved in collisions, the perturbation δf i from the equilibrium density f i of species i caused by the movement of the bubble wall is given by where w is the boost factor of the wall, and + (−) refers to fermions (bosons).
In Eq. ( 29), the CP even term is first order in derivatives, while the CP odd term is second order in derivatives; thus we can solve the CP even and odd parts separately.Following Ref. [24], we introduce the following definition The evolution of f i is described by the Boltzmann equation where L[f i ] is the Liouville operator and v g is the group velocity determined by WKB dispersion relation [22] The C[f i ] is a model-dependent collision integral associated with the interaction rate of the thermal bath [77].The terms in the fluid equation can be written as the average over-phase space of the form [21,24] where f ′ 0+ (m = 0) can be written as Plasma velocities can be defined as The second-order CP odd chemical potential is defined by the difference between the second-order chemical potential of the particle and its anti-particle, and a similar definition follows for corresponding plasma velocities, The zeroth and first momenta of the collision integral can be written in terms of inelastic rate Γ inel and total interaction rate Γ tot by [77] ⟨C For the generation of the baryon asymmetry, the first step is to produce asymmetry in left-handed quarks.We consider the effects of the strong sphaleron process, W -scattering, top Yukawa interaction, helicity flip, and Higgs number violation with the rate of Γ ss , Γ W , Γ y , Γ m , and Γ h respectively.The last two processes are relevant only in the broken phase.The transport equation for chemical potentials of the left-handed top quark, the conjugate of the right-handed bottom quark, left-handed bottom quark, Higgs bosons, and the corresponding plasma velocities are given as follows [21,22,24]: • Charge conjugation of right-handed top quarks (t c ) • Left-handed bottom quarks (b) • Higgs S t denotes the source term of the top quark that can be written as The source term for the bottom quark can be neglected due to the suppression factor m 2 b /m 2 t ∼ 10 −3 .Thermal transport coefficients are defined as with the expectation values given by and the distribution function defined as The third equation in Eq. ( 51) is a Taylor expansion; hence, it is valid only for small bubble wall velocity v W .We assume v W = 0.1.The transport equation with full dependence on the wall velocity is provided in Ref. [72].The values for the strong sphaleron rate, top Yukawa rate, Higgs number violating rate and rate for spin-helicity flipping rate for the top quark are given by [21,24,78,79] Γ ss = 4.9 where z is the distance.The W exchange rate can be approximated as the total Higgs interaction Γ W = Γ tot h .Finally, the asymmetry in left-handed quarks is converted into baryon asymmetry by electroweak sphaleron transition which can be calculated as [77] where Γ ws ≃ 1 × 10 −6 T is the weak sphaleron rate estimated by lattice calculation [80] and g ⋆ ≃ 106.75 is the effective degrees of freedom at the electroweak scale.The chemical potential for left-handed quarks µ B L is given by We solved the top transport equation and estimate the baryon asymmetry of the Universe with BSMPT v2 [23]. 6

VII. BARYON ASYMMETRY IN THE C2HDM
In this section, we estimate the baryon asymmetry generated via electroweak baryogenesis using the semi-classical force method.A key ingredient in this calculation is the estimation of the bubble profile.As discussed in Section VI A, it is usual in the literature to parametrize bubble profile by the kink profile [20][21][22]24].In this section, in addition to deriving the BAU in the C2HDM framework, we pay close attention to the viability of the kink profile by comparing it with the bubble profile obtained by explicitly solving the tunneling equation using CosmoTransitions [62].
In Fig. 5, we compare the magnitude of baryon asymmetry estimated using these two bubble profiles at the nucleation temperature with the color code representing the distribution probability.First, we observe that the C2HDM can satisfy the observed baryon asymmetry, matching the observed value η obs .However, these points are rare in our parameter space scan.We highlight one of these points in Fig. 5 as benchmark point 2 (BP2) and detailed define it in Tab.I. 7 Second, we observe in Fig. 5 that for most of the points in the parameter space, the kink profile solution leads to larger values than the profile obtained by solving the tunneling equation.In addition, there is a non-negligible fraction of points where the kink profile overestimates the asymmetries by a few orders of magnitude in comparison to the profile from the solution of the tunneling equation. 7The VEV-insertion approximation (VIA) [81] has been found to produce baryon asymmetry values that are two to three orders of magnitude larger than the semi-classical force method adopted in the current study [24], displaying a larger number of points satisfying η obs .Several works have raised criticisms about the validity of the approximations used in this alternative method.One particular argument is that the expansion utilized in deriving the source term for the top quark in the VIA approach may encounter limitations due to the substantial mass of the top quark [24,72].It is important to highlight that recent improvements have been made in treating the source term in the VIA method [82].In most cases, we can understand the difference in asymmetry using two profiles by looking at the behavior of the source term Eq. ( 47) in front of the bubble wall.Specifically, the sign of ∂ z θ t in front of the bubble wall determines the sign of the source term S t , thereby influencing the overall asymmetry.In most cases, a negative (positive) ∂ z θ t results in a positive (negative) source term S t , leading to a positive (negative) asymmetry.The kink profile typically provides a higher value for the top mass around the bubble wall, and thereby a higher magnitude for the source term in Eq. ( 47).This feature is illustrated in Fig. 6 using our benchmark point 1 (BP1) as defined in Tab.I.Even when the change in phase of the top mass θ t has a larger magnitude for the tunneling profile, the value of the top mass is higher for the kink profile, and subsequently, the kink profile has a larger asymmetry.Therefore, the behavior of the top mass is the dominant factor in estimating the magnitude of the asymmetry compared to the phase of top mass θ t .
There are instances where the top mass in the tunneling profile is smaller, but does not quickly drop to zero compared to the kink profile.This translates into the source term being active for a larger wall distance for the tunneling profile compared to the kink profile.
In this case, the magnitude of asymmetry calculated using the tunneling profile would have a larger value compared to the kink profile.The above feature is illustrated in the case of BP2 shown in Fig. 7, where the asymmetry differs by two orders.Once again, we highlight that BP2 can explain the value of the observed baryon asymmetry η obs when using the explicit solution for the tunneling equation.The above characteristic of the tunneling profile will permit significant baryon asymmetry even though the change in CP violating phase is relatively small.
Finally, there are instances where even the sign of the derivative for the CP phase ∂ z θ t near the bubble wall differs between the kink and tunneling profiles.We illustrate this scenario with the benchmark point 3 (BP3) presented in Fig. 8. Near the bubble wall, θ t exhibits an increasing trend for the kink profile, while it displays a decreasing trend for the explicit solution for the tunneling profile.Consequently, the asymmetry is positive for the tunneling profile and negative for the kink profile, despite both profiles having identical endpoints.
These findings emphasize the importance of accurately determining the bubble profile and highlight the discrepancies that can arise when relying on the kink profile approximation.

VIII. SUMMARY
In this work, we explored the phase transition pattern and the feasibility of generating the observed baryon asymmetry of the Universe within the C2HDM framework while considering the theoretical and experimental constraints.First, we carefully examined the essential elements in the shape of the Higgs potential, specifically focusing on the formation of the barrier and the upliftment of the true vacuum state.These factors are critical in facilitating the phase transition from a smooth crossover to a strong first order phase transition.We observe that the intensity of the phase transition is linked to the elevation of the true vacuum relative to the symmetric vacuum state at zero temperature [12][13][14][15].
This phenomenon occurs due to the prevalence of one-loop effects over thermal corrections, particularly when ξ c > 1 [15].However, if the vacuum upliftment measure is too large, the universe becomes trapped in the false vacuum state, rendering no solution for the nucleation temperature Eq. ( 18).This leads to parameter points with ∆F 0 /|F SM 0 | ≳ 0.87 unphysical, which excludes most of the ξ c > 2 points [15,16].Therefore, in electroweak baryogenesis studies, it is crucial to look at the nucleation temperature T n and not just at the critical temperature T c .
When it comes to gravitational wave signals, only a small fraction of the parameter points in the Strong First-Order Electroweak Phase Transition parameter space of the C2HDM can be probed by LISA.However, among the accessible points, those with a higher value of the ξ c parameter display particularly strong gravitational wave signals.Notably, the Type I parameter space points generally offer more promising gravitational wave signals compared to the Type II parameter points in the C2HDM.These differences can be traced to the more stringent flavor constraints imposed on the Type-II scenario that shape the parameter space.
We note that the C2HDM can describe the observed baryon asymmetry η obs , albeit for a limited set of parameter space points.One specific point, BP2, was highlighted as a benchmark that satisfied the observed asymmetry value.Furthermore, we contrast the impact on the baryon asymmetry calculation using two different approaches to describe the bubble profile, namely the usually adopted kink profile parameterization and the explicit solution for the tunneling equation.Our objective was to access the dependency of the resulting value of η B on these two approaches and evaluate their respective contributions to the baryon asymmetry calculation.We found that the majority of points in our parameter space scan yield similar results from both approaches.Nonetheless, a non-negligible portion of points exhibits significant discrepancies between these two methods.Specifically, the kink profile approximation often displays higher asymmetry values compared to the explicit solution obtained from the tunneling equation.In some cases, the discrepancy was by several orders of magnitude.The difference in the asymmetry value for the two profiles was scrutinized in terms of the behavior of the source term in front of the bubble wall.
Undoubtedly, the task of achieving a baryon asymmetry of the universe that aligns with the observed value poses significant challenges.The requirements of a strong first-order electroweak phase transition, substantial CP violation, and stringent theoretical and experimental constraints make the generation of a compatible BAU a formidable task.However, the discrepancies observed in calculations performed using different profile assumptions provide avenues for improving the accuracy of computing the BAU.The comparison between the kink profile and the explicit solution for the tunneling profiles provided valuable insights into the estimation of baryon asymmetry, emphasizing the importance of accurately determining the bubble profile for a more robust analysis of electroweak baryogenesis.

Bubble collisions
Nucleated bubbles undergo collisions that disrupt their spherical symmetry and give rise to gravitational waves.In the β Hn ≫ 1 limit, this can be described well by the thin wall and envelope approximation [85,86].In this approximation, GW is generated mainly from envelopes, and the overlapped region is neglected.The gravitational wave spectrum from bubble collision can be modeled by [87] Ω coll h 2 = 1.67 where k c is the efficiency factor for the bubble collision [88] k c = 0.715α + 4 27 3α 2 Taking into account the redshift of the frequency, the peak frequency of the GW spectrum from bubble collision is given by [89] f env = 0.62 × β/H ⋆ 1.8 − 0.1v w + v 2 where first term in right side represents the peak frequency of the GW spectrum from bubble collision at T n .

Sound waves
The energy released from the phase transition to the plasma can either go to heat or fluid motion.Numerical estimation shows that the energy-momentum tensor of the fluid after bubble collision is similar to that of an ensemble of sound waves [90].Remarkably, these sound waves serve as a notable source of gravitational waves.The gravitational wave spectrum stemming from these sound waves can be modeled by [90][91][92] Ω SW h 2 = 2.65 × 10 where the peak frequency of the GW spectrum, accounting for the redshift factor, can be written as [89], Hz . (B6) For the radiation-dominant universe, the lifetime of the sound waves is finite.Recent studies [93,94] show a suppression factor Υ(τ SW ) due to the finite active period τ SW of sound waves where the efficiency factor for the sound wave spectrum k s is given by [64] k s = α 0.73 + 0.083 The sound wave lasts until the turbulence starts to develop.The duration τ SW is quantified with [90] where R ⋆ ≃ (8π)

MHD turbulence
The energy injected into the plasma can induce turbulence in the fluid if the early universe plasma has an extremely high Reynolds number [88].This turbulent motion can be a source of gravitational waves [96].Additionally, a fully ionized plasma can give rise to a turbulent magnetic field under turbulent motion, leading to GWs.The GW spectrum stemming from turbulence can be parametrized as [97] Ω Hz. (B12) Based on insights drawn from numerical simulation, we adopt k turb = 0.05k s [64].Finally, the peak frequency for GW spectrum arising from MHD turbulence can be written as [97], Hz. (B13)

FIG. 1
FIG. 1: The ratio FIG. 2: The parameter space scan in terms of (∆m H ↑ , ∆m H ↓ ) for Type-I (left panel) and Type-II (right panel).The heat map tracks order parameter ξ c (upper panel), ξ n (middle panel), and log 10 (SNR) (lower panel).The gray points in the upper panel pass all theoretical and current experimental constraints.Black points also show first-order phase transition with 0 < ξ c < 1.

FIG. 3 :
FIG. 3: SNR versus ∆F 0 /F SM 0 color coded with ξ n for the Type I parameter point with ξ n > 1.The dotted line in the plot corresponds to an SNR value of 10, serving as the threshold above which LISA can probe the parameter points.

FIG. 4 :
FIG. 4: The left panel shows energy density U E (r) as a function of r for BP1 as defined in Tab.I.The position of the bubble wall r 0 is determined as maxima of U E (r).In the right panel, we show the dynamic evaluation of CP violating angle of top quark θ t with temperature in the broken phase.

FIG. 5 :FIG. 6 :
FIG. 5: Comparison of the magnitude of baryon asymmetry computed using the profile from the solution of the tunneling equation and the kink profile at T n color coded with the probability distribution.The red color denotes the benchmark points presented in Tab.I.

FIG. 7 :
FIG. 7: The evolution of phase of top mass θ t and square of top mass with the wall distance for the BP2.It takes a larger wall distance z for the top mass value to drop to zero in the tunneling profile compared to the kink profile.This translates into the source term being active for a larger distance and, subsequently, the larger value of asymmetry for the tunneling profile compared to the kink profile.

FIG. 8 :
FIG. 8: The evolution of phase of top mass θ t and square of top mass with the wall distance for the BP3.The phase of top mass near the bubble wall is an increasing function for the kink profile and decreasing function for the tunneling profile; subsequently, both asymmetries have opposite signs.

TABLE I :
Benchmark points for the complex two Higgs double model.