Higgs inflation and cosmological electroweak phase transition with N scalars in the post-Higgs era

We study inflation and cosmological electroweak phase transitions utilizing the Standard model augmented by $N$ scalars respecting a global $O(N)$ symmetry. We observe that the representation of the global symmetry is restricted by the inflationary observables and the condition of a strongly first order electroweak phase transition. Theoretical constraints including the stability, perturbativity and unitarity are used to bound the model parameter space. The Electroweak precision observables and Higgs precisions limit the representation of the symmetry. We evaluate the possibility to simultaneously address the inflation and the dark matter after considering the experimental constraints from the future leptonic colliders. When the $O (N)$ symmetry respected by the N-scalar is spontaneously broken to the $O (N-1)$ symmetry, both the one-step and two-step SFOEWPT can occur within the inflation viable parameter regions, which will be tested by the future CEPC, ILC and FCC-ee. The relation between the number of Goldstones and the SFOEWPT condition depends on phase transition patterns. The situation of Goldstone faking neutrinos and contributing to the dark radiation are investigated.


Introduction
To our knowledge, the cosmic inflation [1][2][3] solves the horizon, flatness and monopole problems of the universe successfully. The primordial density fluctuations generated during the inflation can explain the formation of large scale structure of the universe observed by CMB [4]. The inflation scenario is more attractive when the inflaton field can play some important role in the particle physics, one fascinating scenario is the Higgs inflation [5,6], where the inflaton is the SM Higgs being observed by LHC [7,8], that can leads to spontaneously breaking of the electroweak symmetry SU(2) L ×U(1) Y . The original Higgs inflation suffers from the unitarity problem at high scale around ∼ O(10 13 ) GeV [9][10][11][12][13][14][15][16][17][18][19]. Therefore, the cosmologists and particle physicists turn to investigate the singlet scalars assistant case utilizing Higgs-portal [20,21] interactions.
In the post-Higgs era, how does the electroweak symmetry breaking occurs comes to be the hot research topic of particle physicists [22]. One outstanding mechanism is the electroweak phase transition. A strong first order electroweak phase transition does not happen in the Standard Model of particle physics [23,24], to obtain a strongly first order one, the new physics is always necessary where the cubic Higgs couplings deviation might be able to be probed at high energy colliders [22,23]. The most simplest way to realize the SFOEWPT can be extending the SM with an additional real singlet scalar [25,26] or complex singlet scalar [27,28].
The existence of the dark matter is supported by overwhelming astrophysical and cosmological observations. For the most simplest DM case with the hidden sector coupled with the SM through only the Higgs-portal interaction of |H| 2 S 2 [29][30][31], the one-step type SFOEWPT can be realized with a large quartic coupling that would directly leads to the under-abundant situation [32]. The straightforward approach to ameliorate the situation of requirement of large quartic coupling can be adding to the SM N-scalars with O(N) symmetry, therefore the one-step SFOEWPT can be realized with a lower magnitude of Higgs-portal interaction of |H| 2 S i S i (i = 1, ..., N) [33][34][35][36]. The previous studies of Ref. [34] indicates that the one-step SFOEWPT cannot be addressed together with a correct relic density unless the quartic couplings of |H| 2 S i S i and masses of S i are non-universal, rather one can expect a rather large N(e.g., N>50 ) to overcome the problem. Here, the additional N hidden scalars might also alleviate the hierarchy problem through positive contributions to radiative corrections of the Higgs boson mass and therefore satisfy the Veltman conditions [39].
For the Higgs inflation in the Higgs-portal scenarios, the typical quartic scalar couplings are required to be around ∼ O(10 −1 ) [20,21,38]. For that magnitude of couplings, the previous studies of N-scalars with an exact O(N) symmetry suggest the one-step SFOEWPT can be realized for a larger N, with typical frequency of ∼ O(10 −3 − 10 −1 )Hz gravitational wave signals can be probed [34][35][36] together with a substantial triple Higgs couplings deviation to be probed by the future colliders [35,36]. The WIMP DM situation in the classical scale invariant N-scalars model with O(N) symmetry [39] are ruled out for N > 4 even for rather large quartic coupling [40]. Here, we stress that the largeness of the quartic couplings can not accommodate inflations though prefer one-step SFOEWPT.
With the increasing pressure from the direct detection experiments, the relic abundance might be partially rather than fully saturated by WIMP DM. In the one-step SFOEWPT scenario, the dimensional six operators can be used to improve the vacuum barrier. While, in the Higgs-portal scenario with the hidden sector obey Z 2 symmetry, the dimensional six operator lose the ability to lift up the vacuum barrier at finite temperature due to the cancellation of the different contributions given the Z 2 symmetry is spontaneously broken [41]. Novertheless, in the simplest scenario with real singlet scalar serves as WIMPs DM candidate, one can expect the SFOEWPT through two-step types [32,42].
In this work, we study the inflation and EWPT in the Higgs-portal N-scalars model. When the N-scalar fields respect the O(N) symmetry, then all the N-scalars can serve as DM candidates. On the other hand, when the symmetry is broken to O(N − 1) spontaneously, there will be N − 1 Goldstones that can fake the effective neutrinos from the viewpoint of cosmology. In Section. 2 we revisit the O(N) scalar extended SM, and introduce the possible case with the O(N) being spontaneously broken, where we present the relevant theoretical and experimental constraints. The effects of the N scalars on the inflation are explored in Section. 3.1 after imposing theoretical constraints. The cosmological electroweak phase transitions are studied around the Higgs inflation valid parameter regions in Section. 3.2. In Section. 3.3, the dark matter situation in the O(N) model has been explored after taking into account the bounds from future colliders. We explore the EWPO constraints on N and the ability of the future lepton colliders to probe or exclude these parameter regions in Section. 4. The effective neutrino constraints are used to set bounds on the number of Goldstones, the ability of Goldstones to contribute to dark radiations are estimated also provided they gain masses from non-renormalizable gravity effects. We conclude with the Section. 5.

The Models
In this work, we study two scenario of N singlet scalar extended standard models. In the case of the N singlet scalar with O(N) symmetry, one can anticipate the O(N) symmetry broken at finite temperature and restored at zero temperature. The zero temperature tree-level Lagrangian is given by . After the spontaneously symmetry breaking of the Electroweak symmetry, the mass term of S i is given as m 2 S i = µ 2 s + λ hs v 2 /2. The N singlet scalar S i with O(N) being spontaneously broken to O(N − 1) symmetry, which is another model setup, in this case we use "s" instead of "S" to differentiate it from the O(N) scenario. The minimization conditions of the potential can be obtained when EW symmetry is broken and the O(N) being broken along the direction of s, with other directions being s i (i = 1, ..., N − 1) 1 , The mass matrix is given by In order to diagonalize the mass matrix, we introduce the rotation matrix R = ((cos θ, sin θ), (− sin θ, cos θ)) with tan 2θ = −λ hs vv s /(λ h v 2 − λ s v 2 s ) to relate the mass basis and field basis, The mass squared eigenvalues are

Theoretical constraints
Firstly, due to the additive property of the scalar contributions to the beta functions of quartic couplings, especially when the N is larger, one need to aware the possible perturbativity problem of the model at high scale when one perform the inflation analysis. We impose the following conditions to preserve the perturbativity, Since the perturbative unitarity rely on the quartic coupling of the tree-level potential, one can expect the bound is the same for O(N) and O(N → N − 1) scenarios, being given by [36], with the couplings are related with masses and VEVs through Eq. 2.8 for O(N → N − 1) case. To prevent the unbounded from bellow of the scalar potential, the vacuum stability conditions should be satisfied, A simple analysis of these theoretical limits on the quartic couplings at the Electroweak scale are given in Fig. 1 and Fig. 2 for O(N) and O(N → N − 1) models. The perturbativity roughly bounds the upper limit of λ hs,s , the sharp of the boundary is set by the unitarity bounds, the lower limit of the quartic couplings λ hs,s is given by the stability conditions where more parameter spaces are allowed by λ sh > −2 √ λ h λ s in comparison with λ sh > 0. In the Fig. 2, it's converted to the bounds on the m h 2 and v s correspondingly. In the following sections, we require the perturbativity, unitarity, and the stability from the electroweak scale to planck scale for the safe of inflation and EWPT studies, that have been evaluated with the renormalization group equations list in Sec.5.

Higgs precisions
Integrating out the heavy scalar fields results in the dimension-six operators,  O 6 can alleviate the triple Higgs coupling and is crucial for the strong first order electroweak phase transition [43]. For the Wilson coefficients generated at tree-level, we have [44,45], Thus one obtain a universal shift of all Higgs couplings. Which therefore induce the correction to the e + e − → hZ associated production cross section [47], , (2. 18) which has been defined as the fractional change in the associated production cross section relative to the SM case. For the O(N) symmetry with m S > m h , the Higgs wavefunction renormalization shift the SMlike Higgs couplings to other SM particles by c N H v 2 /(2m 2 S ) ∼ Nλ 2 hs v 2 /(2λ s m 2 S ). Which results in the constraint on c H and therefore N, λ hs,s from the LHC [41] as well as ILC, CEPC, and FCC-ee [46].
On the other hand, the operator O H induces the operator combinations O W +O B and O T operators through RGE [44,47], which results in the S and T parameters We set bounds on m S(h 2 ) and the "N" using the electroweak fit in Ref. [48] , using the feynman diagram method in Ref. [49]. One can obtain severely constraints on the θ with increasing of m h 2 . For the O(N → N − 1) scenario, the mixing angle and the heavy Higgs masses are subjective to the bounds coming from the LHC Higgs data, which force the mixing angle θ to be |cos θ| ≥ 0.84 [26]. After including of the current LHC and High-luminosity LHC Higgs production rates together with the EWPO, a moderate of θ ∼ λ 2 hs v 2 /(4λ s m 2 h 2 ) = 0.2 can be safety [41]. Here, we stress that: though for the O(N → N − 1) scenario, when fit the Higgs data we need to take into account the invisible decay of the Higgs to the N − 1 Goldstone bosons, one can firstly perform the Higgs fit without including the effects of changing of SM Higgs decay width and then constrain the number of Goldstone bosons N − 1 with the invisible decay fit results from [50]: B BSM < 0.34 at 95% CL. For the case of O(N) symmetry being broken to O(N − 1) at zero temperature, we have the following lagrangian to describe the triple scalar interactions, with h i, j denotes h 1,2 and N = 2, ..., N, the relevant triple scalar couplings are given bellow, From which, when the m h 2 > m h 1 , the decay widths of the SM-like Higgs and the second Higgs are 2 Indeed, if one want the heavy Higgs take part in the EWPT process, it cannot be highly decoupled.
given by For the case in which m h 1 > m h 2 one need take into account the decay of h 1 → 2h 2 , with decay width being given by The invisible decay of SM Higgs can be used to set upper bounds to the number of the Goldstones and the mixing angle θ. At 95% CL, the LHC (ATLAS+CMS) set B inv < 34% [50], see Fig. 3 for the constraints. With the increase of v s , more parameter space of θ is allowed.

High scale and low scale phenomenologies
We first study the cosmic inflation with large scale fields. With the temperatures of the universe cooling down the low scale fields physics come to us: the possibility to implement the cosmological strong first order phase transitions, and the dark matter physics.

The Higgs inflation with N singlet scalars
The action in the Jordan frame is 3 where M p is the reduced Planck mass, R is the Ricci scalar, ξ h,s define the non-minimal coupling of the h, S-field. The quantum corrected effective Jordan frame Higgs potential at large field values (h) can be written as which is evaluated along the higgs axis, where the scale is µ ∼ O(h) ≈ h. We impose quantum corrections to the potential following Ref. [51,52]. After the conformal transformation, and a field redefinition we obtain In this work we consider only h serves as inflaton for the O(N) model and the O(N → N − 1) model [53][54][55]. The metric in this case is given by The inflationary action in terms of the canonically normalized field χ is therefore given as with the potential in terms of the canonically normalized field χ as where the new field χ are defined by for h− inflations [53]. Note that λ h and ξ h have a scale (h) dependence. The slow-roll parameters are then given by The field value at the end of inflation χ end is obtained when ε = 1, and the horizon exit value χ in can be calculated by assuming e-folding numbers being 60 between the two periods, Then, one can relate the inflationary observables of spectrum index n s and the tensor to scalar ratio of r with the slow-roll parameters at the χ in , Meanwhile, one can determine the non-minimal gravity couplings using the constraint coming from CMB observations [4], with the amplitude of scalar spectrum fluctuations ∆ 2 R being calculated as We note that though the above formula can describe the slow-roll inflation of O(N) and O(N → N − 1) scenarios, the inflation dynamics for the two scenarios are different due to the mass spectrum and the relation between scalar masses and quartic couplings which will be explored in detail in the following sections.

O(N) scenario
First, in Fig. 4, we show the Higgs inflation feasible parameter regions in the plane of (λ s , λ hs ) for different N after imposing the theoretical constraints up to planck scale as aforementioned in Sec. 2.1. Where, the upper and lower bounds of λ hs are mostly coming from perturbativity and unitarity, and stability conditions. The inflation feasible (λ s , λ hs ) ranges are the largest for N = 1 compared with other N cases. The feasible ranges are reduced with the increased N and are overlap for the two neighbour N expect for N = 1 and N = 2. The decreasing of the inflation valid area with the increase of N is due to the fact that: a larger N leads to more contributions of λ hs to λ h at the inflation scale through RG running, and therefore the stability, perturbativity and unitarity set the lower and upper bounds of λ hs . We plot the RG running of the quartic scalar couplings for the case of N = 7, 10, 13 in Fig. 5. The perturbativity of quartic coupling and unitarity can be violated due to RG running of couplings as shown in the right panel of Fig. 5.
We show the Higgs inflation valid parameter spaces in Fig. 6, the feasible ranges are reduced with the increased N and are overlap for the two contiguous N. The left panel indicates that the λ hs increases with the increase of λ s for the O(N → N − 1) case, which is different from the O(N) scalar model cases. A interesting triangular shape shows up due to the bounds on m h 2 and v s from perturbativity, unitarity, and stability, together with the relations among quartic couplings and the Higgses masses as well as VEVs. Here, it should be note that the lower value of m h 2 is set by the stability bounds, the quartic Higgs coupling will be negative for m h 2 <330 GeV. That set the lower boundary of inflation. And the upper boundary is set by the magnitude of v s asistanted by m h 2 , and the upper bound on m h 2 < 600 GeV is required to fulfilled the perturbatitity and unitarity conditions at high scale. For explicitly, we explain the property in the

Electroweak phase transitions
With the temperature cooling down, the universe can evolve from symmetric phase to the symmetry broken phase. The behavior can be studied with the finite temperature effective potential with particle physics models [56]. Through which one can obtain the critical classical field value and temperature being v C and T C to characterize the critical phases. Roughly speaking, a SFOEWPT can be obtained when v C /T C > 1, then the electroweak sphaleron process is quenched inside the bubble and therefore one can obtain the net number of baryon over anti-baryon in the framework of EWBG. For the uncertainty of the value and possible gauge dependent issues we refer to Ref. [57]. The finite temperature effective potential include the tree level scalar potential, the Coleman-Weinberg potential, and the finite temperature corrections [58]. The tree level scalar potential for O(N) is obtained directly from Eq. 2.1,  Here, we drop the subscript since all N directions are the same and we assume only one direction get VEV during the EWPT process. The "S" should be the "s" for the O(N → N − 1) scenario to indicate the possible symmetry breaking direction with other directions s i (i = 1, ..., N − 1) do not get VEV during the EWPT process. For the O(N) scalar model, the one-loop Coleman-Weinberg potential for the scalar parts is given by  . The finite temperature corrections to the effective potential at one-loop are given by [56], and with the upper (lower) sign corresponds to bosonic (fermionic) contributions. Here, the above integral J B,F can be expressed as a sum of them second kind modified Bessel functions K 2 (x) [60], The thermal masses/corrections are given by, for the O(N) case, and for the O(N → N − 1) case one needs to replace the "S" by "s" and replace the thermal mass of the Higgs fields by Last but not least, the resummation of ring (or daisy) diagrams are also crucial for the evaluation of v C and T C with the finite temperature effective potential [61], which is given by and for O(N) and O(N → N − 1) cases. Here, again, we list only the contributions of scalar contributions for V T and V ring , the others particle fields contributions are the same as the SM, see Ref. [59,61]. It should be note that, the counter terms can keep the VEVs of the potential from shift caused by the V CW , we add that parts follows Ref. [59]. Then the critical parameters of EWPT can be calculated when there are two degenerate vacuums with a potential barrier. Due to rich vacuum structures of the potential at finite temperatures, there can be one-step or multi-step phase transitions. A SFOEWPT can be realized in the first or the second step in the two-step scenario. We will investigate two scenarios: one-step and two-step phase transition in O(N) and the O(N → N − 1) scalar models. Before the detailed study, we firstly recall the one-step strong first order phase transition conditions on Wilson coefficients of the dimensional-six operator [43,62,63], (3.34) In the Fig. 8 we show the physics picture for the energy barrier with different dimensional six operator Wilson coefficients. A suitable c 6 = c cr 6 is needed to obtain a proper vacuum barrier to seperate two degenerate vacuum at critical temperature T c , therefore make the SFOEWPT feasible. It should be note that, with the spontaneous symmetry breaking of O(N → N − 1), the two contributions of the c 6 from sh 2 and s 3 terms cancels each other when one following the method of Ref. [45] and therefore the tree-level induced dimensional six operator disappears. Which is the same as in the SM+1 singlet case being studied in Ref. [64]. And in this case the dimensional six operator shows up at loop level which is too small to affect the EWPT dynamics. This property can explain why the SM+1 singlet scalar with Z 2 does not prefer one-step strong first order EWPT, and here one may needs to pursue the two-step types where the DM can assistant to achieve the SFOEWPT [32,42]. In this work, we confirm the same property in the O(N) symmetry and O(N → N − 1) cases. The inflation occurs for small quartic scalar couplings, that motivate us mostly focusing on the two-step types.

O(N) scenario
The corresponding critical temperature and critical field value for one-and two-step EWPT types in O(N) (see the Fig.9) can be evaluated through the following degeneracy conditions, and for one step and two step EWPTs. Here the s C is the O(N) broken direction, which is analogy to the O(N → N − 1) scenarios. The survey of the one-step EWPT in the O(N) scenario shows negligible effects of the quartic couplings λ s , and the quartic coupling between the SM Higgs and the O(N) scalars S i ( λ hs ) should be large in order to make the SFOEWPT occurs, which is not favored by the slow-roll inflation. Generally, the parameter spaces with a larger λ hs one can only have a SFOEWPT where the inflation is invalid, the perturbativity of quartic coupling and unitarity are violated due to RG running of couplings as explored in Sec. 3.1.1. The numerical survey of the two step EWPT shows that a strong first order EWPT requires N≥7 in the inflation favored parameter regions. Unfortunately, as can be seen from Fig.20, the inflation valid N should be smaller than 4 after imposing the constraints from EWPO. Which therefore shout down the window to realize the SFOEWPT.

O(N → N − 1) scenario
We demonstrate the one-step and two-step phase transition in Fig.10.
For the one-step case, one can realize strong first order EWPT with a smaller λ hs with increasing of N, as shown in Fig.11. Which means that the Goldstones contributions to the EWPT is notable at finite temperature. While this property disappear in the two-step scenario as can be seen in the bottomright panel. The two-step phase transition can be strongly first order with a relatively lower magnitude of λ hs in comparation with the one-step scenario. The analysis also shows that different from one-step scenario, the strong first order EWPT occurs more easy at relatively small λ hs for two-step case, that provide the possibility to obtain the SFOEWPT in the inflation favored regions. Moreover, our study demonstrates that the rate of s B /s A can be larger or smaller than 1 for different N, as shown in Fig.12.
As shown in Fig. 6 previously, for the O(N → N − 1) scalar model, the slow-roll inflation are almost excluded for N ≥ 5 in the (λ s , λ sh ) plane. We show the parameter regions that can valid Higgs inflation and SFOEWPT together with Higgs cubic couplings and Higgs decay widths in O(N → N − 1) model in the Fig.13 after taking into account perturbativity, unitarity, and stability conditions from Electroweak scale to Planck scale. The inflation and one(two)-step SFOEWPT are allowed by the Higgs invisible decay bounds from LHC [50] for different N, which is marked by cyan. The twostep SFOEWT shows no obvious relation with N due to s A can be higher or lower than s B as shown in Fig.12. In the SFOEWPT allowed parameter spaces with relatively larger quartic couplings, the slowroll inflation is not allowed, this property is caused by the bound on m m 2 and v s set by perturbative unitarity and stability from weak scale to planck scale imposed by us when perform the inflation analysis with the RGEs. There are larger parameter spaces of (λ s , λ hs ) allowed by two-step SFOEWPT  in comparison with the one-step SFOEWPT. In the regions, we find that the ratio of the triple Higgs couplings (r 3h 1 increases with the increase(decrease) of λ hs (λ s ). With the SM Higgs resonance searched by SM Higgs pairs production, one can estimated the cross section with respect to the SM case being , and Γ tot h 2 are shown by red, blue, and brown contours, respectively. Dashed, solid and dotted lines stand for N = 1, 2, 3, respectively. Cyan regions present the allowed regions by the invisible Higgs decay bounds from B inv < 0.34 [50], in which the light and deep colors correspond to N = 2, 3, respectively. The feasible regions of inflation are shown by green color regions, and the orange regions represent one-and twostep SFOEWPT for the left and right panels, respectively. For those colors, a deeper color corresponds to a larger N for N = 1, 2, 3.
Γ tot h 2 m h 2 , the resonance interference explored in Ref. [41] can be safely neglect here. We postpone the detailed collider search of the parameter spaces to a separate publication.

Dark matter/radiations
With the assumption of instantaneous reheating, one can estimate the reheating temperature when the decay of the inflaton starts competing with expansion H ∼ Γ h1,h2 for h/s inflation [65].
here g * ≈ 100 is the number of relativistic degrees of freedom in the Universe during the reheating epoch. Freeze out of cold dark matter requires x f ≡ m DM /T f s ≈ 20, here the thermal history can occurs as T R > T C > T f s > T BBN to account for the EWPT and reheating as well as the successful Big Bang Nucleosynthesis (BBN) (with typical temperature of a few MeV). 4 The freeze-out temperature T f s being smaller than the strong first order electroweak phase transition temperature T C , set the m DM < 20T C ∼ 2 TeV with T C being around ∼ O(10 2 ) GeV.

O(N) scenario: dark matter
When the O(N) is kept at zero temperature, the N-singlet scalars can all serve as dark matter candidates. With thermal averaged annihilation cross sections being the same as in Ref. [53] for each S i and we drop the subscript here. For SM Higgs pair final states one have the cross section for gauge boson final states are, and the fermion pair final states cross section is given by, The formula of spin independent cross section is given by [30] σ S SI = λ 2 The dark matter direct detection constrains the dark matter masses and the quartic Higgs-DM couplings after taking into account the rescale effects supposing the evaluated dark matter relic density will not oversaturated the DM relic abundance, (3.45) With Lee-Weinberg method [66], we can expect Ω S i h 2 ∼ 1/σv rel ∼ m 2 S i /λ 2 hs for large dark matter masses. Previous studies show that the mass region of Higgs-portal real 1-singlet scalar DM case is excluded up to ∼TeV scale by Xenon1T to obtain a correct relic density [30,38].
It should be note that the future Linear collider constraints would set the dark matter mass in TeV mass regions (∼ O(1 − 10) TeV) in the inflation feasible region as shown in Fig. ??. In this case, the seagull diagram dominate the contributions to the dark matter pair annihilations. We show the annihilation cross section in Fig. 14 to illustrate that. It's easy to see that the contributions of the Here we point out that m s ∼ µ s due to the Electroweak symmetry is still kept at temperature higher than T c within the framework of EWBG. Then, if the relic abundance is partially saturated by S i , one need a larger H-S quartic coupling λ hs . Though this case could be good for the SFOEWPT, the large m S in fact would decouple from the EWPT as been studied in Ref. [34]. Furthermore, the largeness of the λ hs may results in the perturbativity unitarity problem of quartic couplings of λ h,hs as shown in Fig. 5 after taking into account the RG running effects and therefore shut down the possibility to realize inflation.

O(N → N − 1) scenario: Goldstone and dark radiation
Suppose gravity violates global symmetries, then the Goldstone boson may acquire a mass through nonpertubative gravitational effects [67,68]. The non-perturbative gravity effects can break the O(N) symmetry at M P scale through the lowest high dimension operators, i.e., dim-5 operators, and induce mass terms to Goldstone bosons and make N − 1 majoron like particles,  In this mass region, we can expect the majoron here decay to diphoton through the non-minimal gravity couple term to break the O(N − 1) symmetry as in the [69]. Follow which we find the Goldstone bosons here is long-lived, with τ ∼ Γ −1 ∼ 10 46 s, they can survive until the recombination era and may contribute to the Universe radiation density at the time of recombination or BBN. Now, we consider the possibility of the Goldstone contributing to the dark radiation following Ref. [70]. The effective neutrino number can be expressed in terms of the Goldstone decoupling temperature as [71,72], with ∆N s N−1 = 4(N −1)/7 due to there is (N −1) Goldstone bosons decoupled at T > T d s N−1 and present before the recombination eras, the effective number of relativistic degrees of freedoms g * (T d ν ) = 43/4 and g * (T d s N−1 ) = 57/4 supposing that the Goldstone bosons decouple just before muon annihilation. One can constrain the number of Goldstone as in Fig. 15 using the recent 1σ experimental data N e f f = 3.36 ± 0.34 [73] .  To study how the Goldstone decouples from the thermal bath, we follow Ref. [70]. For the heavy Higgs contributions are typically small, one need to focus on the light Higgs case alternatively, c.f., m h 2 < 2m h 1 . When the decay width Γ h 2 m h 2 and in the small mass region the cross-section of thermal annihilation to µ + µ − , as shown in Fig.16, is given by, which leads to the constraints on v s and m h 2 , as seen in Fig. 17. The invisible decay of h 1 requires a smaller θ or a lower magnitude of N. In the resonance enhanced region, 2m µ < m h 2 < 4 GeV, using which set lower bounds on the mixing angle of θ, see Fig. 18. For a smaller v s , the invisible decay of the h 1 set the upper limits on θ depends on the number of N as shown in the left panel of Fig. 18. For the case of m µ < m h 2 < 2m µ , we have, θ v s =1000 GeV Figure 18: Decouple conditions valid regions for 2m µ < m h 2 < 4 GeV, regions in the dashed square frame are allowed by the B inv < 0.34 [50]. The decay width of h 2 ( Γ tot h 2 ) labeled on the contours is shown in units of KeV. θ m h 2 =3m μ /2 Figure 19: Decouple conditions valid parameter regions for m h 2 = 3m µ /2, regions in the dashed square frame are allowed by the B inv < 0.34 [50]. The decay width of h 2 ( Γ tot h 2 ) labeled on the contours is shown in units of KeV.
Requiring the Goldstone bosons annihilation process contribute to the equivalent neutrino numbers, we obtain the bounds on mixing angle and the v s , see Fig. 19. A larger N allows a smaller θ to meet the decoupling conditions due to more Goldstone contributions.

The future lepton colliders constraints
The study of Ref. [46] shows that the CEPC, ILC, and FCC-ee can probe the new physics parameter spaces ( through the e + e − → hZ process ) much better than that from LHC. The c N H v 2 /m 2 S i is bounded to be smaller than 0.0038, 0.0034, and 0.0028 by CEPC with luminosity of 5 ab −1 , ILC with all center of mass energies, and FCC-ee with luminosity of 10 ab −1 . With which, we can obtain the bounds on the hidden scalar masses of m S i ∼ O(1 − 10) TeV, as shown in Fig. 20. As been explored in the Sec. 3.3.1, in this mass regions the N scalars cannot explain the correct relic density. And in N = 1 of the O(N → N − 1) case, which is the Higgs-portal 1-singlet scalar case, the corresponding mixing angle | sin θ| was constrained to be 0.062, 0.058 and 0.052 at 95% C.L. at CEPC with luminosity of 5 ab −1 , ILC with all center of mass energies, and FCC-ee with luminosity of 10 ab −1 [46]. We map that values to constraints our model parameter spaces. Here, we point out that considering the high sensitivity of CEPC, ILC, and FCC-ee, a much smaller mixing angle θ can be probed, see Fig. 21. At that moment, one obtains a smaller λ hs , and therefore the stability problem can easily preclude the possibility of slow-roll Higgs inflation. For the case of N ≥ 2 of the O(N → N − 1) model, the ILC from the Higgs Strahlung process set B inv < 1% [74], the FCC-ee set B inv < 0.5% [75,76], and CEPC set 0.14% [77]. Fig.22 imply a larger mixing angle θ is allowed for a larger v s , which corresponds to the heavy Higgs decouple cases. Generally, to make the Higgs inflation feasible, a relatively larger mixing angle θ is required to enlarge the value of λ hs and therefore to save the Higgs quartic coupling from the vacuum unstability problem. Firstly, the increasing of θ can leads to the perturbativity problem of λ h , and thus we need a lower magnitude of m h 2 . secondly, a larger v s leads to a smaller λ hs , therefore to avoid the stability problem we need a larger N. In the narrowed down region of mixing angle θ coming from the future e + e − colliders constraints, the mixing between the heavy Higgs and the SM Higgs is negligible. In this super-weak coupe scenario one have a far smaller of λ hs to obtain a SFOEWPT, thus one need a much larger N to amplified the effects of the N − 1 Goldstone to the potential in order to obtain a SFOEWPT. On the other hand, the Goldstone fake effective neutrino situation would be changed a lot due to the decouple conditions allowed parameter spaces can be covered by the B inv from ILC, FCC-ee and CEPC as aforementioned.

Conclusions and discussions
In this work, we study inflation and cosmological electroweak phase transitions utilizing the SM augmented N scalars respecting O(N) symmetry or O(N → N − 1) symmetry. With the assistant of N scalars that couple to the SM Higgs, the stability problem can easily been remedied up to the inflation scale. Therefore, the interaction of the Higgs-portal quartic scalar couplings shouldn't be too small depending on the number of scalars. The quartic scalar coupling perturbativity and the unitarity constraints set the upper bounds on the quartic scalar couplings. The Higgs inflation valid parameter regions reduces with the increases of the number of N mostly due to the perturbativity and the unitarity bounds.
The Electroweak precision observables set severe bounds on the parameters spaces of both O(N) and O(N → N − 1) scenarios. For the O(N) case, the largest number of the scalars that can valid inflation is bounded to N ≤ 3, which make the strong first order EWPT unreachable, no matter onestep or two-step paths. The future e + e − colliders, such as ILC, FCC-ee, and CEPC can set the masses of N-scalars being ∼ O(1 − 10) TeV scale. When the O(N) symmetry is an exact symmetry, in principle, all of the N scalars can serve as WIMP DM candidates. While, no way to expect the Nscalar WIMP DM can saturate the correct relic density here. And, in that dark matter mass regions, the freeze out happens earlier that the EWPT process, thus the only relevant DM annihilation process is S i S i → hh with m h (T ≥ T f s ) ∼ 0.
When the O(N) symmetry is spontaneously broken, one obtains N − 1 Goldstones, and the rest one scalar mixed with the SM Higgs. Therefore, the invisible Higgs decay is very powerful to set the bound on the scalar numbers N and the mixing angle θ. With one moderate θ = 0.2 allowed by EWPO and Higgs precisions as well as invisible Higgs decay bounds set by LHC run-I, we explore the possibility to realize Higgs inflation and EWPT through one-step and two-step types. In the parameter regions where one can obtain successful slow-roll Higgs inflation and SFOEWPT, the triple Higgs couplings λ h 2 h 1 h 1 and λ h 1 h 1 h 1 increase with the increasing of λ hs , and the decay width of the two Higgses are not larger enough to introduce large interference effects in the resonant mass regions of Higgs pair productions. The future high precision e + e − colliders, such as ILC, FCC-ee, and CEPC, set the mixing angle and the number of Goldstones N − 1 far more stringent than that from the LHC, which could probe the inflation and SFOEWPT valid parameter regions. The N − 1 Goldstones can fake the effective number of neutrinos and contribute to dark radiations supposing they obtained ∼ O(1) eV masses induced by non-renormalizable gravity effects. The roughly estimations independent of the mixing angle indicates the N is bounded to be smaller than five at 3σ. The dark radiations calculations indicates that Goldstones decouples from the thermal bath at different mass range of a small m h 2 . The parameter spaces of this scenario can be covered by the projected ILC, FCC-ee, and CEPC in the future.