Prospects of gravitational waves in the minimal left-right symmetric model

The left-right symmetric model (LRSM) is a well-motivated framework to restore parity and implement seesaw mechanisms for the tiny neutrino masses at or above the TeV-scale, and has a very rich phenomenology at both the high-energy and high-precision frontiers. In this paper we examine the phase transition and resultant gravitational waves (GWs) in the minimal version of LRSM. Taking into account all the theoretical and experimental constraints on LRSM, we identify the parameter regions with strong first-order phase transition and detectable GWs in the future experiments. It turns out in a sizeable region of the parameter space, GWs can be generated in the phase transition with the strength of $10^{-17}$ to $10^{-12}$ at the frequency of 0.1 to 10 Hz, which can be detected by BBO and DECIGO. Furthermore, GWs in the LRSM favor a relatively light $SU(2)_R$-breaking scalar $H_3^0$, which is largely complementary to the direct searches of a long-lived neutral scalar at the high-energy colliders. It is found that the other heavy scalars and the right-handed neutrinos in the LRSM also play an important part for GW signal production in the phase transition.


Introduction
The discovery of a Higgs boson at the Large Hadron Collider (LHC) heralds the completion of the standard model (SM) [1,2] and a great hope for the discovery of new physics. Obviously, the completion of the SM naturally leads to the quest of microscopic structure to its next chapter, which will be further searched by the LHC [3]. In the long list of questions which might be the key to the next chapter, a few are interesting and crucial.
For example, what is the dynamics for the electroweak (EW) symmetry breaking, what is the origin of mass of neutrinos [4], how are the parity and CP symmetries broken, and what is nature of dark matter and dark energy [5], etc. To answer these questions has been motivating various new physics models beyond the SM (BSM) at the TeV scale.
In the history of early universe, from the Planck time to today, phase transitions might have occurred when the symmetries at different energy scales are broken. For example, the symmetry breaking of grand unified theory (GUT) and supersymmetry (SUSY) breaking can induce the corresponding phase transitions at the GUT scale and SUSY breaking scale. For new physics beyond the SM, new dynamics and a larger symmetry are usually introduced at the TeV region or a higher-energy scale. Such new physics models are of transitions in the LRSM and the resultant features of the corresponding GWs. Compared with the recent and former study [102], the new things of this paper lie in the following aspects: (i) we have implemented the correct EW vacuum conditions [103] and set α 2 = 0 (α 2 is a quartic coupling in the scalar potential Eq. (2.2)), (ii) we have taken into account more recent LHC experimental bounds, which are collected in Table 1 and Fig. 1, (iii) we have found more general parameter space where the strong FOPT can occur and detectable GWs can be produced, and (iv) we have also explored the complementarity of GW probes of LRSM and the direct searches of the heavy (or light long-lived) particles in the LRSM at the high-energy colliders, and examined how the self couplings of SM Higgs can be affected in the LRSM.
With all the theoretical and experimental limits taken into consideration, it is found that the strong FOPT at the right-handed scale v R in the LRSM favors relatively small quartic and neutrino Yukawa couplings, which corresponds to relatively light BSM scalars and right-handed neutrinos (RHNs), as seen in Figs. 2, 3 and 9. The scatter plot in Fig. 5 reveals that the phase transition in the LRSM can generate GW signals with the strength of 10 −17 to 10 −12 , with a frequency ranging from 0.1 to 10 Hz, which can be probed by the experiments BBO and DECIGO, or even by ALIA and MAGIS. The GW spectra for five benchmark points (BPs) are demonstrated in Fig. 8, which reveals that the GW signal strength and frequency are very sensitive to the value of ρ 1 . Although some other quartic and neutrino Yukawa couplings are very important for the GW production, the quartic coupling ρ 1 plays the most crucial role and it also determines the mass of SU (2) R -breaking scalar H 0 3 . In the parameter space where it does not mix with other scalars, the scalar H 0 3 couples only to the heavy scalars, gauge bosons and RHNs in the LRSM [104], which makes it effectively a singlet-like particle, and thus the experimental limits on it are very weak [105,106]. As presented in Fig. 10, the GW probe of H 0 3 is largely complementary to the direct searches of H 0 3 at the high-energy colliders [104] as well as the searches of H 0 3 as a long-lived particle (LLP) at the high-energy frontier [105,106]. In addition, in a sizeable region of parameter space, the strong FOPT and GWs are sensitive to a large quartic coupling λ hhhh of the SM-like Higgs, which is potentially accessible at a future high-energy muon collider [107].
The rest of the paper is organized as follows. In Section 2 we briefly review the minimal LRSM and summarize the main existing experimental and theoretical constraints on the BSM particles in this model. Phase transition are explored in Section 3, and the GW production is presented in Section 4. Section 5 focuses on the complementarity of the GW probes of LRSM and the collider signals of LRSM. After some discussions, we conclude in Section 6. For the sake of completeness, the masses and thermal self-energies are collected in Appendix A, and the conditions for vacuum stability and correct vacuum are put in Appendix B.

Left-right symmetric model
The basic idea of LRSMs is to extend the EW sector of SU (2) L × U (1) Y of the SM gauge group to be left-right symmetric, i.e. SU (2) L × SU (2) R × U (1) B−L . Various LRSMs have been proposed to understand the parity symmetry and CP breaking of the SM, the origin of masses of matters or even DM candidates and the matter-antimatter asymmetry of the universe. The main differences between these LRSMs could be in the gauge structure, the scalar fields, the matter contents, and/or the seesaw mechanisms.
The most popular, or say conventional, LRSM is the version with a Higgs bidoublet Φ, a left-handed triplet ∆ L and a right-handed triplet ∆ R [99][100][101] When the right-handed triplet ∆ R acquires a vacuum expectation value (VEV) v R , the gauge symmetry SU (2) L × SU (2) R × U (1) B−L in the LRSM is broken to the SM gauge group SU (2) L × U (1) Y . Two triplets ∆ L and ∆ R are introduced to give Majorana masses to the active neutrinos and RHNs, respectively, which enables the type-I [108][109][110][111][112] and type-II [113][114][115][116][117] seesaw mechanisms for the tiny neutrino masses. The SU (2) R × U (1) B−L symmetry can also be broken only by a right-handed doublet H R [118,119]. In this case, heavy vector-like fermions have to be introduced to generate the SM quark and lepton masses via seesaw mechanism (see also [120]). There are also LRSM scenarios with inverse seesaw [121,122], linear seesaw [123,124], or extended seesaw [125][126][127][128] in the literature. Cold DM is not included in the conventional LRSM (a light RHN can only be a warm DM candidate [129]), but it is easy to add a fermion or boson multiplet, where the lightest neutral component is naturally stabilized by the residual Z 2 symmetry from U (1) B−L breaking [130,131]. Alternatively, based on the gauge group SU the hypercharge in the SM and Y R the "right-handed" counter partner), heavy RHNs can be the cold DM candidate [132][133][134].
In this work, we focus on the minimal LRSM with one bidoublet Φ and two triplets ∆ L and ∆ R in the scalar sector. The most general scalar potential in the LRSM can be written as [135] 2) whereΦ = σ 2 Φ * σ 2 (with σ 2 the second Pauli matrix). Required by left-right symmetry, all the quartic couplings in the potential above are real parameters. The CP violating phase δ associated with α 2 is shown explicitly. At the zero temperature, the neutral components of the scalar fields can develop nonzero VEVs, i.e.
where θ κ and θ L are CP violating phases. The two bidoublet VEVs are related to the EW In light of the hierarchy of top and bottom quark masses m b m t in the SM, it is a reasonable assumption that κ 2 κ 1 [135]. There are three key energy scales in the LRSM, i.e. the right-handed scale v R , the EW scale v EW and the scale v L which is relevant to tiny active neutrino masses via type-II seesaw. Furthermore, from the first-order derivative of the scalar potential (2.2), v L is related to the EW and right-handed VEVs via [113,135,136] where ξ = κ 2 /κ 1 . Due to the tiny masses of active neutrinos, it is a good approximation to set v L = 0, therefore we will set β i = 0 so as to simplify our discussions below.
With v L = 0, there are only two energy scales in the LRSM, i.e. the EW scale v EW and the right-handed scale v R . In light of the hierarchy structure v EW v R , a two-step phase transition is supposed to occur in the LRSM. In the early universe, the temperature is so high T v R that the symmetry SU (2) L × SU (2) R × U (1) B−L is restored. As the universe keeps expanding, the temperature decreases. When the temperature is lower than a critical temperature but much higher than EW scale, i.e. v EW T ∼ v R , ∆ 0 R develops a non-vanishing VEV and the gauge symmetry When the temperature becomes lower than the EW scale T ∼ v EW , Φ 0 1,2 obtain their VEVs and the symmetry is further broken into the electromagnetic (EM) group U (1) EM .
After symmetry breaking at the v R scale, we can rewrite the bidoublet Φ in terms of two SU (2) L doublets, i.e. Φ = iσ 2 H * 1 |H 2 . Then the bidoublet relevant terms in the potential (2.2) can be recast in terms of H 1, 2 : where the mass terms are respectively Although the potential in Eq. (2.5) seems to be very similar to that in a general 2HDM [137], there are still some obvious differences: In presence of the scale v R , all the states predominately from the heavy doublet H 2 are at the v R scale, and their masses are degenerate at the leading-order, which is clearly distinct from the 2HDMs, where all the scalars in the 2HDMs are at the EW scale, and the BSM scalar masses depend on different quartic couplings [137].
In the LRSM, the BSM particles include the heavy W R and Z R bosons, three RHNs N i (with i = 1, 2, 3), neutral CP-even scalar H 0 1 and CP-odd A 0 1 , singly-charged scalar H ± 1 predominately from the bidoublet Φ, neutral CP-even scalar H 0 2 and CP-odd A 0 2 , singlycharged scalar H ± 2 and doubly-charged scalar H ±± 1 mostly from the left-handed triplet ∆ L , and the neutral CP-even scalar H 0 3 and doubly-charged scalar H ±± 2 mostly from the righthanded triplet ∆ R . Thorough studies of the scalar sector of LRSM at future high-energy colliders can be found e.g. in Ref. [103-106, 135, 138-157]. In this paper, we assume that the gauge coupling g R for SU (2) R can be different from the gauge coupling g L for SU (2) L , which might originate from renormalization group running effects such as in the D-parity breaking LRSM versions [158].

Theoretical Constraints
For completeness, we collect all the theoretical constraints on the gauge and scalar sectors of the LRSM in the literature, which will be taken into consideration in the calculations of phase transition and GW production below.
• Perturbativity limits: In some versions of the LRSM, the right-handed gauge coupling g R can be different from g L [158]. As the gauge couplings have the relationship (with g BL the gauge coupling for U (1) B−L ) (2.9) the gauge couplings g R and g BL can not be either too large or too small if we want them to be perturbative. Renormalization group running these gauge couplings up to a higher energy scale put more stringent limits on them. Perturbativity up to the GUT scale requires that the ratio r g ≡ g R /g L to satisfy [159] 1 0.65 < r g < 1.60 . (2.10) Furthermore, as the masses α 3 /2v R of H 0 1 , A 0 1 and H ± 1 (cf. Table 5 in Appendix A) are severely constrained by the neutral meson mixings (see Section 2.2 and Table 1), perturbativity also implies an lower bound on the v R scale [159]: v R 10 TeV . (2.11) For v R below this value, α 3 is so large that all the quartic and gauge couplings will hit the Landau pole very quickly before reaching the GUT or Planck scale [151,154,160,161].
• Unitarity conditions: The parameters in the potential (2.2) should satisfy the unitarity conditions [154] when we consider the scattering amplitudes of the scalar fields at the high-energy scale √ s µ 2 i (for simplicity we neglect here the effects of all the scalar masses). In other words, the partial wave amplitudes should not violate the bound of unitarity so as to guarantee that the probability is conserved. The tree-level unitarity conditions turn out to be [154] • Vacuum stability conditions: The vacuum stability conditions require that [154,161,162] (see also [163]) • Correct vacuum criteria: After the spontaneous symmetry breaking, all the scalar fields have to form some specific structure in the phase space such that we reside in the correct vacuum, i.e. the vacuum with the lowest VEV in the potential [103,156]. For completeness, the correct vacuum criteria have been collected in Appendix B, which are obtained with the assumption α 2 = 0. Therefore, we will set α 2 = 0 throughout this paper.
In the limit of κ 2 κ 1 v R , in Eq. (2.7), the quadratic coefficient of H 2 term is proportional to α 3 v 2 R /2, thus the heavy doublet scalars H 0 1 , A 0 1 , H ± 2 will obtain a mass of α 3 /2v R at the leading-order. To get the correct EW vacuum, a necessary condition is m 2 This yields an upper bound of ξ. Approximately, we have

Experimental constraints
All the current LHC limits on the BSM particles in the LRSM are collected in Table 1 and also depicted in Fig. 1. Here are more details: • At the LHC, the W R boson in the LRSM can be produced via the right-handed charged quark currents. After its production, it can decay predominately into two quark jets (including thetb channel) and RHNs plus a charged lepton, i.e. W R → jj,tb, N ( * ) i α (with α = e, µ, τ ). If the RHNs are lighter than the W R boson, as a result of the Majorana nature of RHNs, the same-sign dilepton plus jets W R → N → α β jj constitute a smoking-gun signal of the W R boson [164]. Assuming g R = g L , the current most stringent LHC data require that the W R mass m W R > (3.8 − 5) TeV for a RHN mass 100 GeV < m N < 1.8 TeV [165,166]. The dijet [167,168] and tb [169,170] limits are relatively weaker, which are respectively 4 TeV and 3.4 TeV. The strongest W R limit of (3.8 − 5) TeV is presented in Fig. 1.
• The most stringent limits on the Z R boson is from the dilepton data pp → Z R → + − .
The current dilepton limit on a sequential Z boson is 5.1 TeV [171]. Following e.g.
Ref. [159], one can rescale the production cross section times branching fraction σ(pp → Z → + − ) for the sequential Z model, which leads to the LHC dilepton limit of 4.82 TeV on the Z R boson in the LRSM. This is shown in Fig. 1 as the Z R limit. There are also dijet searches of the Z boson, however the corresponding limits are relatively weaker [167,168].
• At the leading order, the scalars H 0 2 , A 0 2 , H ± 2 and H ±± 1 from the left-handed triplet ∆ L have the same mass [142] (see Table 5). The doubly-charged scalar H ±± 1 can decay into either same-sign dilepton or same-sign W bosons, i.e. H ±± 1 → ± α ± β , W ± W ± , which constitute the most promising channels to probe ∆ L at the LHC, and the branching fractions BR(H ±± 1 → ± α ± β ) and BR(H ±± 1 → W ± W ± ) depend on the Yukawa coupling f L and the left-handed triplet VEV v L . Assuming H ±± 1 decays predominately into electrons and muons, the current LHC limits are around 770 to 870 GeV, depending on the flavor structure [172]. In the di-tauon channel H ±± 1 → τ ± τ ± , the LHC limit is relatively weaker, i.e. 535 GeV [173]. 2 If the doubly-charged scalar H ±± 1 decays predominately into same-sign W bosons, the LHC limits are much weaker, around 200 to 220 GeV [174]. There are also some searches of singly-charged scalar H ± 2 → τ ± ν at the LHC [175][176][177]. However these searches assume H ± 2 is produced from its interaction with top and bottom quarks, therefore these limits are not applicable to H ± 2 in the LRSM which does not couple directly to the SM quarks. The strongest same-sign dilepton limits of (530 − 870) GeV on H ±± 1 (and also on other scalars from ∆ L ) is shown in Fig. 1. 2 As the singly-charged scalar H ± 2 and doubly-charged scalar H ±± 1 are mass degenerate at the leading order in the LRSM, here we have adopted the combined LHC limit from the pair production pp → H ++ 1 H −− 1 and the associate production pp → H ±± 1 H ∓ 2 . In these two channels, the separate channels are respectively 396 GeV and 479 GeV [173].
• As the W R boson is very heavy, the TeV-scale right-handed doubly-charged scalar H ±± 2 decays only into same-sign dileptons. The couplings of H ±± 2 to the photon and Z boson have opposite signs, therefore the production cross section of H ±± 2 at the LHC is smaller than that for the left-handed doubly-charged scalar H ±± 1 . Rescaling the LHC13 cross section of H ±± 1 by a factor of 1/2.4, The same-sign dilepton limits on H ±± 2 turn out to be 271 to 760 GeV for all the six combinations ee, eµ, µµ, eτ, µτ, τ τ of lepton flavors, which is presented in Fig. 1.
• The scalars H 0 1 , A 0 1 and H ± 1 from the bidoublet Φ are degenerate in mass at the leading order. H 0 1 and A 0 1 has tree-level flavor-changing neutral-current (FCNC) couplings to the SM quarks, and contribute to K − K, B d − B d and B s − B s mixings significantly. As a result, their masses are required to be at least (10 − 25) TeV, depending on the nature of left-right symmetry (either generalized parity or generalized charge conjugation), the hadronic uncertainties [142,[178][179][180] and the potentially large QCD corrections [181]. The stringent FCNC limits on the heavy bidoublet scalars is shown in Fig. 1.
• The neutral scalar H 0 3 from the right-handed triplet ∆ R is hadrophobic, i.e. it does not couples directly to the SM quarks in the Lagrangian. It can be produced at the LHC and future higher energy colliders either in the scalar portal through coupling to the SM Higgs (and the heavy scalars H 0 1 and A 0 1 ), or in the gauge portal via coupling to the W R and Z R bosons. Therefore the direct LHC limits are very weak [105,106]. However, when it is sufficiently light, say at the GeV-scale, H 0 3 can be produced from (invisible) decay of the SM Higgs or even from the meson decays [105,106]. More details can be found in Section 5.2.
• The RHNs in the LRSM can be either very light, e.g. at the keV scale to be a warm DM [129] candidate, or very heavy at the v R scale, and there are almost no laboratory limits on their masses, although their mixings with the active neutrinos are tightly constrained in some regions of the parameter space [182]. For simplicity, in the following sections we will set the masses of RHNs to be free parameters and neglect their mixings with the active neutrinos.
To be complete, the masses of 100 GeV scale SM particles, i.e. the SM Higgs h, the top quark t and the W and Z bosons, are depicted in Fig. 1 as horizontal black lines. See Fig. 7 for complementarity of GW prospects of the BSM particle masses and the current experimental limit.

One-loop effective potential
To study phase transitions in the LRSM, we consider the effective potential at finite temperature, which includes contributions of the one-loop corrections and daisy resummations.   Table 1, indicated by the blue and pink arrows, with the heights of the horizontal lines denoting the ranges of experimental limits. The horizontal black lines are the masses of SM Higgs h, top quark t, and W , Z bosons.
Renormalized in the MS scheme, the effective potential can be cast into the following form [183] is the Coleman-Weinberg one-loop effective potential [184], and V T =0 1 and V D are the thermal contributions at finite temperature. The V T =0 1 term includes only the one-loop contributions, and V D denotes the high-order contributions from daisy diagrams. In Eq. (3.1) the sum runs over all the particles in the model. The scalar mass matrices m 2 i (κ i , v R ) in the LRSM can be found in Ref. [135], and the corresponding thermal self-energies Π i (T ) are provided in Appendix A. As for the fermions, we consider only the third generation quarks and three RHNs. In the LRSM their masses are respectively with y t, b the Yukawa couplings for top and bottom quarks in the SM, M N the RHN masses and y N the corresponding Yukawa coupling. In the following study, for the sake of simplicity, we will assume three RHNs are mass degenerate and does not have any mixings among them. The degrees of freedom g i and constants C i in Eq. (3.1) are given by 3 2 ), for fermions , (3, 5 6 ), for gauge bosons , with λ = 1 (2) for Weyl (Dirac) fermions, and the functions J − (J + ) for bosons (fermions) are defined as In the limit of small x 2 = m 2 /T 2 , we can use the approximations [183]: In this paper we focus on the phase transition at the v R scale, thus as an approximation all the effects of SM components on the symmetry breaking SU can be neglected. Neglecting the daisy contributions, the effective potential V eff can be written down explicitly in the following form [6]: where D, T 0 , E and ρ T can be expressed by the model parameters as where M X is the mass for the particle X, and µ is the renormalization scale. Since there are lots of scalars in the LRSM, we deliberately separate their contributions from the vector bosons and RHNs. The contributions of scalars for each of the terms in Eq. (3.9) to (3.12) can be written in terms of the scalar masses via (3.13) (3.14) It should be pointed out that all the masses in Eqs. (3.13) to (3.16) depend upon the right-handed VEV v R instead of v. It is observed that the RHNs can also contribute to the symmetry breaking SU (2) R × U (1) B−L → U (1) Y via affecting the parameters D, T 0 and ρ T , while the parameter E receives only contributions from the scalars and gauge bosons.
As seen in Eqs. (3.12) and (3.16), the parameter ρ T receive not only tree-level contribution from the quartic coupling ρ 1 which corresponds to the H 0 Table 5), but also loop-level contributions from the heavy scalars, gauge bosons and RHNs in the LRSM. In particular, when the quartic coupling ρ 1 is small, or equivalently the scalar H 0 3 is much smaller than the v R scale, which is the parameter space of interest for phase transition and GW production in the LRSM (cf. Figs. 2, 3 and 8), the loop-level contributions in Eq. (3.12) might dominate ρ T . Furthermore, ρ T depends also on the gauge coupling g R via the heavy gauge boson masses M W R and M Z R .
To have strong FOPT, the cubic terms proportional to −ET v 3 are crucial. In the limit of E → 0, the phase transition is of second-order. In the SM, the effective coefficient E of φ 3 term is dominated by the gauge boson contributions, while in the LRSM, it receives contributions from both the scalars and gauge bosons, As a result of the large degree of freedom in the scalar sector of LRSM, it is remarkable that the scalar contributions to E can even be much larger. The order parameter describing the FOPT is given by v c /T c , where v c is the non-vanishing location of the minimum at the critical temperature T c at which the effective potential V eff has two degenerate minima. In the EW baryogenesis [185][186][187], to avoid the washout effects in the broken phase within the bubble wall, a strong FOPT is typically required to satisfy the following condition Therefore, it is justified to neglect the contributions of SM particles to the phase transition at the right-handed scale v R , since their masses m SM are at most close to v EW and their contributions are suppressed due to their tiny couplings to the right-handed triplet.
For given v R and heavy particle masses in the LRSM, the two key parameters T c and v c can be obtained from the effective potential (3.1) by requiring the two conditions V ef f (T c ; v c ) = V ef f (T c ; 0) and v c = 0. In the numerical evaluations, we change the temperature from a sufficiently highly energy scale, say v R , toward lower values around the EW scale. A reasonable critical temperature T c for the phase transition SU (2) R × U (1) B−L → ×U (1) Y is assumed to be within this range. The dependence of v c /T c on the parameters in the LRSM is exemplified in Fig. 2, where in the numerical calculations we have included all the contributions in Eq. (3.1).
Taking into account all the theoretical and experimental constraints in Section 2, we first consider scenarios with the simplifications λ 2 = λ 3 = λ 4 = α 1 = α 2 = 0. In order to identify the parameter space where the phase transition is of first-order, we calculate v c /T c at the critical temperature T c with different values of the quartic couplings ρ 1 , ρ 2 , ρ 3 − 2ρ 1 , α 3 . When we calculate the dependence of v c /T c on two of the quartic couplings, all others are fixed in the way that their corresponding scalar masses equal the W R mass, and the gauge coupling g R = g L . To be concrete, we have set the renormalization scale µ to be the v R scale in Eq. (3.1). The corresponding results are shown in the first three panels of Fig. 2. The dependence of v c /T c on the couplings ρ 1 and α 3 , ρ 3 − 2ρ 1 and ρ 1 , and ρ 2 and  Table 5), the dependence of v c /T c on the quartic couplings in Fig. 2 can also be understood effectively as the dependence of v c /T c on the mass−v R ratios Through the gauge boson masses M W R and M Z R , the parameter v c /T c depends also on the gauge coupling ratio r g , or equivalently on the right-handed gauge coupling g R . This is shown in the lower right panel in Fig. 2; as seen in this figure, the v c /T c limit on ρ 1 has a moderate or weak dependence on r g , depending on the value of ρ 1 .
Given the information on v c /T c in Fig. 2, a few more comments are now in order: • As seen in Fig. 2, a strong FOPT in the LRSM require a relatively small quartic coupling ρ 1 0.07 for the parameter space we are considering, which is qualitatively similar to the SM case where a light Higgs boson (say M h < 80 GeV) is needed in order to have a first-order EW phase transition [189]. It turns out that a small ρ 1 (and resultantly light H 0 3 ) is not only crucial for the prospects of GWs in future experiments (cf. Fig. 8), but also triggers rich phenomenology for the searches of LLPs at the high-energy colliders and dedicated detectors [105,106].
• The phase transition at the v R scale occurs when the neutral component ∆ 0 R of the right-handed triplet ∆ R develops a non-vanishing VEV v R . As a result, the strong FOPT is more sensitive to the mass of H 0 3 , or equivalently to the value of ρ 1 , than other heavy scalar masses. This is also clearly demonstrated in the plots of Fig. 2. As seen in the upper left, upper right and lower left panels, the quartic coupling α 3 , ρ 3 − 2ρ 1 and α 3 can reach up to order one, while ρ 1 0.1 in the Fig. 2.
• Although the quartic couplings α 3 , ρ 3 and ρ 3 − 2ρ 1 is less constrained by the FOPT than the critical coupling ρ 1 , as seen in the first three panels of Fig. 2, if either of these couplings is sufficiently large, it will invalidate the strong FOPT at the v R scale, no matter how small ρ 1 is. Meanwhile, the white areas in the plots of Fig. 2 indicate that in these regions the perturbation method starts to break down and theoretical predictions become more difficult.
In Fig. 2 we have fixed some parameter in the LRSM and vary two of them. To see more details of the correlation of v c /T c and the parameters in the LRSM, we take a more thorough scan of the parameter space of the LRSM. To be specific, we adopt the following ranges: and apply all the theoretical and experimental constraints in Section 2. Here follows some comments: • We have chosen ξ = κ 2 /κ 1 = 0.001 in order to satisfy the theoretical constraint in Eq. (2.15).
• We have chosen α 2 = 0 in order to meet the requirement of the correct vacuum conditions given in Eq. (B.1).
• It is known from Fig. 2 that the strongly FOPT need a small ρ 1 , therefore we have chosen ρ 1 < 0.5.
• ρ 3 − 2ρ 1 has set to be larger than zero, as it corresponds to the masses of the lefthanded triplet scalars (see Table 5).
• The quartic coupling α 1 is not a free parameter here, as it is related to λ 1 and the SM coupling λ via Eq. (5.1). As α 2 /4ρ 1 is always positive, it turn out that the quartic coupling λ 1 ≥ λ 0.13. The resultant scatter plots of v c /T c are presented in Fig. 3 as functions of the parameters ρ 1 and α 3 . The data points of strong FOPT with v c /T c > 1 are shown in red while those with v c /T c < 1 are in blue. When we set v R = 10 TeV and take the FCNC limit of M H 0 1 > 15 TeV [142], the quartic coupling α 3 should meet the condition The region shaded by the light pink in the left panel of Fig. 3 is excluded by such conditions. It is found that only a small amount of the data points can survive and have strong FOPT. When the v R scale is higher, say v R = 20 TeV, the quartic coupling α 3 is significantly smaller, i.e. α 3 > 1.13. The region denoted by the light pink shaded region in the right panel of Fig. 3 is excluded. Then there will be more points that can have a strong FOPT with v c /T c > 1, as clearly shown in the right panel of Fig. 3.

Gravitational waves
The thermal stochastic GWs can be generated by three physics processes in phase transition [190]: collisions of bubbles, sound waves (SWs) in the plasma after the bubble collision, and the MHD turbulence forming after the bubble collision. For non-runaway scenarios, GWs are dominated by the latter two sources [190], and the corresponding GW spectrum can be approximated as The SW contribution has the form of [25] where f is the frequency, g * and H * are respectively the number of relativistic degrees of freedom in the plasma and the Hubble parameter at the temperature T * , v w is the bubble wall velocity, α describes the strength of phase transition, β/H * measures the rate of the phase transition, and is the fraction of vacuum energy that is converted to bulk motion. The peak frequency f SW is approximated by The MHD turbulence contribution is [30,191] where κ MHD 0.05κ v is the fraction of vacuum energy that is transformed into the MHD turbulence, h * is the inverse Hubble time at the GW production (red-shifted to today), and is given by Hz , (4.6) and the peak frequency is (4.7) As shown in the formula above, the gravitational wave spectrum from FOPTs are generally characterized by two parameters related to the phase transition, namely α and β [31]. The parameter α is defined as the ratio of the vacuum energy density * released at the phase transition temperature T * to the energy density of the universe in the radiation era, i.e.
where * is the latent heat and can be expressed as The ∆V eff denotes the difference of potential energy between the false vacuum and true vacuum, i.e. ∆V eff = −V eff (0, T ) + V eff (v, T ), which can be simply determined by T * and the parameters of LRSM. The parameter β describes the rate of variation of the bubble nucleation rate during phase transition, and its inverse describes the duration of phase transition. To describe rate of the phase transition, a dimensionless parameter β H * is defined from the following equation where S 3 denotes the three-dimensional Euclidean action of a critical bubble. The T * denotes the temperature when the phase transition is ended and can be determined by requiring that the probability for nucleating one bubble per horizon volume equals 1, i.e.
The parameters α and β set respectively the strength and time variation of GWs during the phase transition, and their typical values in the LRSM are shown respectively in the left and right panels of Fig. 4. As demonstrated by the data points, the value of α varies roughly from 0.001 to 0.1, and β/H * can range from 10 2 to 10 4 . In the numerical calculations, all the data points in Fig. 4 have strong FOPT.
Assuming the bubble wall velocity v w ∼ 1, the corresponding GW signals of the data points in Fig. 4 are shown in Fig. 5. The correlation of the ratio v c /T c and GW signal peaks are presented in the left panel. We can read from Fig. 4 and the left panel of Fig. 5 that with large v c /T c the value α is typically larger, thus yielding stronger GW signals. The GW strength and frequency peaks are shown in the right panel of Fig. 5. The potential sensitivities of LISA [35,36], TianQin [33], Taiji [34], ALIA [37], MAGIS [38], BBO [40], DECIGO [39], ET [42], and CE [41] are also depicted in the right panel of Fig. 5. As seen in this figure, the frequency peak in the LRSM can range from 10 −1 to 10 2 Hz. Furthermore, there are some data points of the LRSM with frequencies in the range of roughly from 0.1 to 10 Hz and GW strength larger than 10 −17 , which can be detected in the future by BBO and DECIGO, or even by ALIA and MAGIS.  Figure 5. GW peaks for the data points in Fig. 4, as function of v c /T c (left) and frequency f (right). Also shown in the right panel are the prospects of LISA [35,36], TianQin [33], Taiji [34], ALIA [37], MAGIS [38], BBO [40], DECIGO [39], ET [42], and CE [41]. can reach up to few times 10 TeV, with their lower mass limits roughly round the experimental constraints in Section 2.3 (see also Table 1 and   Table 2. The resultant v c , T c , T * and the parameters α and β/H * are also shown in Table 2. The GW spectra h 2 Ω as function of the frequency f for the five BPs are presented in Fig. 8. There are a few comments on the five BPs. • It is clear in Fig. 8  and M N but different M H 0 3 can be probed in the future by BBO and DECIGO, and even by ALIA and MAGIS. It seems that the H 0 3 mass M H 0 3 , or equivalently the quartic coupling ρ 1 , is crucial for the GWs in the LRSM. The BPs (like BP5) with a heavier H 0 3 , or equivalently larger ρ 1 , tends to generate a small α and large β, and thus produce weaker GW signals with a larger frequency. This is consistent with the findings in Ref. [102]. The BPs BP1 and BP2 with H 0 3 mass below TeV-scale can produce GWs of order 10 −13 with frequency at around 0.1 Hz, far above the prospects of BBO and DECIGO. The BP4 with a 2 TeV H 0 3 can only produce GWs of order 10 −16 with frequency peaked at 1 Hz, which can be marginally detected by BBO and DECIGO.  Figure 8. The same as in the right panel of Fig. 5, but for the five BPs in Table 2. Table 2. Five BPs studied in this paper. Parameters not shown in the table are set to be v R = 10 TeV, ξ = 10 −3 , λ 1 = 0.13, α 1 = α 2 = λ 2 = λ 3 = λ 4 = 0. Their GW spectra are shown in Fig. 8. It is also noticed that all these BPs are non-runaway scenarios in term of the criteria defined in Eq. (25) of [190]. The suppression factor Υ in the last row is defined in Eq. (6.1) [196,197]. heavier in BP5 than in BP4, while all other parameters are the same. As seen in Fig. 8, the GW signal in BP5 is so weak that it can escape the detection of all the planned GW experiments in the figure. This reveals that the masses M H 0 2 , M H ±± 2 and M N , or equivalently the couplings ρ 3 − 2ρ 1 , ρ 2 and y N , are also important for GW production in the LRSM. More data points in the numerical calculations reveal that the coupling α 3 is also very important for the GW signals in the LRSM.

Complementarity of GW signal and collider searches of LRSM
In spite of the large number of BSM scalars, fermions and gauge bosons in the LRSM and the larger number of quartic couplings in the potential (2.2), it is phenomenologically meaningful to examine the role of some couplings, or equivalently the BSM particle masses, in the strong FOPT and the subsequent GW production in the early universe, as well as the potential correlations of GWs to the direct laboratory searches of these particles and the SM precision data at the high-energy colliders. In this section, we will elaborate on (i) the effects of the quartic coupling λ 1 in the scalar potential (2.2) which corresponds to the self-coupling λ in the SM, and (ii) the complementarity of GW signal, the collider searches of (light) H 0 3 and the heavy (or light) RHNs in the LRSM.

Self-couplings of SM-like Higgs boson in the LRSM
It is interesting to examine how the self-coupling λ of the SM-like Higgs boson h can be affected by the BSM scalars in the LRSM. The SM-like Higgs mass square, the trilinear coupling λ hhh and the quartic coupling λ hhhh in the SM and LRSM are collected in Table 3.
Comparing the mass square of h in the SM and LRSM, we can approximately identify the following relation among the SM and LRSM quartic couplings [104,151] As seen in the third column of Fig. 3, the trilinear coupling λ hhh of the SM-like Higgs in the LRSM only differs from the SM value by a small amount of ξ ∼ 10 −3 [104,151]. On the contrary, the quartic coupling λ hhhh in the LRSM might be significantly different from the SM prediction: as shown in the last column of In other words, at the leading-order of the approximations of v R v EW κ 1 κ 2 , the difference of quartic coupling of SM-like Higgs boson in the SM and LRSM is dominated by the α 2 1 /16ρ 1 term. As the FOPT and GW in the LRSM favor a small ρ 1 coupling, the difference in Eq. (5.2) tends to be significant for sufficiently large α 1 .
Adopting the parameter ranges in Eq. (3.19) and taking into account the theoretical and experimental limits in Section 2, the scatter plots of the quartic coupling λ hhhh and the couplings ρ 1 , α 1 and y N are shown respectively in the left, middle and right panels of  Figure 9. Scatter plots of λ 1 /λ and ρ 1 (left), α 1 (middle) and y N (right), with the blue points have v c /T c < 1 and the red ones v c /T c > 1. Fig. 9, where the data points with strong FOPT v c /T c > 1 is shown in red, while those with v c /T c < 1 are in blue. It is very clear in Fig. 9 that the deviation of the quartic scalar coupling λ 1 from the SM value λ is always positive and can be very large, even up to the order of 10, as expected in Table 3 and Eq. (5.2). We can also read from the left and middle panels of Fig. 9 that a large deviation of the quartic coupling of SM-like Higgs need a relatively small ρ 1 and/or large α 1 . As given in Eq. (3.12), a large y N tends to decrease ρ T , thus increasing the value of v c /T c . However, if y N is too large, say y N 1.5, a negative ρ T will be obtained which leads to a non-stable vacuum. Thus, the phase transition and GW in the LRSM favor a Yukawa coupling y N ∼ O(0.1) to O(1).
On the experimental side, the combined results of di-Higgs searches can be found e.g. in Refs. [198,199]. Data from LHC 13 TeV with a luminosity of 36 fb −1 only set a weak constraint λ hhh /λ SM hhh ∈ (−5, 12). The LHC 14 TeV with an integrated luminosity of 3 ab −1 can probe the trilinear coupling of SM Higgs within the range of λ hhh /λ SM hhh ∈ (0.7, 1.3) [200], while the future 100 TeV collider with a luminosity of 30 ab −1 can help to improve the sensitivity up to λ hhh /λ SM hhh ∈ (0.9, 1.1) [201]. However, this is not precise enough to see the deviation of trilinear coupling in the LRSM, which is of order 10 −3 or smaller.
Although the quartic coupling measurements can not be greatly improved at hadronic colliders [201,202], a future muon collider with the center-of-mass energy of 14 TeV and a luminosity of 33 ab −1 can probe a deviation of the quartic Higgs self-coupling at the level of 50% [107]. This can probe a sizable region of parameter space in Fig. 9. probe a mass range of roughly 100 GeV up to 3 TeV, while the searches of a long-lived H 0 3 at the high-energy colliders can cover the mass range of 10 GeV down to 100 MeV. As a new avenue to probe the phase transition in the LRSM, GWs are sensitive to a wide mass range of H 0 3 , from the 10 GeV scale up to 10 TeV, which is largely complementary to the searches of (light) H 0 3 at the high-energy colliders. Note that one of the important decay modes of H 0 3 is the RHN channel, i.e. H 0 3 → N N , which will induce the strikingly clean signal of same-sign dilepton plus jets [104,148,153]. The heavy RHNs can also be produced through their gauge couplings to the W R and Z R bosons, e.g. the smoking-gun Keung-Senjanović signal pp → W R → N ± → ± ± jj at the high-energy pp colliders [164]. If the RHNs are very light, say below 100 GeV scale, the decay widths of RHNs will be highly suppressed by W R mass, which makes the RHNs longlived [205,206]. The light long-lived RHNs can be searched directly at the high-energy colliders via displaced vertex, or even from meson decays [207][208][209][210]. The prospects of RHNs at the high-energy colliders and in meson decays depend largely on the heavy scalar or gauge boson masses (see also [211,212]). However, it is worth pointing out that, as seen in Fig. 6, GWs are sensitive to the RHN masses in the range of 200 GeV up to 40 TeV, which is largely complementary to the direct searches of (light) RHNs at the high-energy frontier.

Discussions and Conclusion
Before the conclusion we would like to comment on some open questions in the phase transition and GW production in the LRSM: • In the calculations we have assumed that at the epoch of phase transition the bubbles expanding in the plasma can reach a relativistic terminal velocity, i.e. the nonrunaway scenarios, where the velocity of bubble wall is taken to be v w 1 in our analysis, which corresponds to the denotation case [213]. A recent numerical analysis [214] has revealed that the SW contribution might be suppressed by a factor of 10 −3 in the deflagration case when α > 0.1 where the reheated droplet can suppress the formation of GW signals. While there is no such a huge suppression for the denotation case with α < 0.1, our results could still be valid, although the GW signals might be suppressed by a factor two or three. The bubble wave velocity, in principle, can be computed from the parameters of a given model, as demonstrated in [215][216][217]. Furthermore, according to the recent calculations in Ref. [196], it is found that the finite lifetime of SWs can lead to a suppression factor Υ, which can be parameterized in the following form [197] We have calculated the Υ factors for the five BPs in Table 2, and listed it in the last row of this table. It is observed that the GW signals in these BPs might be suppressed by up to a factor of 6 to 10. It might be interesting to explore how the trilinear couplings expressions model parameters of LRSM can affect the bubble wall velocity and the effects of the suppression factor Υ, which will be a topic for our future study.
• It is remarkable that for the scalar H 0 3 , which is mainly the CP-even neutral component of the right-handed triplet ∆ R , both the theoretical and experimental constraints on it are very weak. As a result, its mass could span a wide range, say from below GeV-scale up to tens of TeV. In the case that all other new particles in the LRSM are heavier than 5 TeV but with a relatively light H 0 3 below the TeV-scale (for instance the BPs BP1 and BP2 in Table 2), at the scale below 1 TeV, the scalar potential of LRSM given in Eq. (2.2) can be reduced to the effective model with the SM extended by a real singlet S, where the scalar potential has the following form: The trilinear and quartic couplings in Eq. (6.2) can be written as functions of the right-handed VEV v R and the quartic couplings in the LRSM, which are collected in Table 4. Obviously, when α 1 is switched off, H 0 3 will not affect the EW phase transition directly, and the EW phase transition should be of second-order as in the SM. When α 1 is switched on, it might be interesting to examine whether the light H 0 3 can affect the phase transitions at both the v R scale and the EW scale. When it is possible, a multi-step strong FOPT could be expected [218].
To summarize, in this paper we have studied the prospects of GW signals from phase transition in the minimal LRSM with a bidoublet Φ, a left-handed triplet ∆ L and a righthanded triplet ∆ R , which is a well-motivated framework to restore parity and accommodate the seesaw mechanisms for tiny neutrino masses at the TeV-scale. We have considered the theoretical limits on the LRSM from perturbativity, unitarity, vacuum stability and correct vacuum criteria, as well as the experimental constraints on the heavy gauge bosons and the BSM scalars. The experimental limits are collected in Table 1 and Fig. 1.
With these theoretical and experimental constraints taken into account, we have analyzed the parameter space of strong FOPT and the resultant GWs in the LRSM. As demonstrated in Figs. 2, 3 and 9, the strong FOPT at the v R scale favors relatively small quartic and Yukawa couplings, which corresponds to relatively light BSM scalars and RHNs. The GWs for some BPs in the LRSM in Fig. 5 reveal that the phase transition in the LRSM can generate the GW signals of 10 −17 to 10 −12 , with a frequency ranging from 0.1 to 10 Hz, which can be probed by the experiments BBO and DECIGO, or even by ALIA and MAGIS. Setting v R = 10 TeV, as seen in Fig. 6, the GWs are sensitive to the following mass ranges: • The heavy bidoublet scalars H 0 1 , A 0 1 , H ± 1 , the scalars H 0 2 , A 0 2 , H ± 2 and H ±± 1 from the left-handed triplet ∆ L , and the doubly-charged scalar H ±± 2 from the right-handed triplet ∆ R , with masses up to tens of TeV, with the lower bounds of their masses roughly set by the experimental limits in Fig. 1.
• The scalar H 0 3 with mass in the range of roughly from 20 GeV up to 10 TeV. As presented in Fig. 10, the GW prospects of H 0 3 are largely complementary to the direct searches of heavy H 0 3 at the LHC and future 100 TeV colliders, and the searches of light H 0 3 from displaced vertex signals at the LHC, future higher energy colliders, and the LLP experiments such as MATHUSLA.
• The RHNs with masses from roughly 300 GeV up to 40 TeV. The GW sensitivity of M N is also largely complementary to the direct searches of prompt signals and displaced vertices from RHNs at the high-energy colliders, as well as the production of RHNs from meson decays.
The GW spectra in Fig. 8 for the BPs in Table 2 shows that the quartic coupling ρ 1 is crucially important for both the frequency and strength of the GW signals in the LRSM, while other couplings such as ρ 2 , ρ 3 − 2ρ 1 , α 3 and y N are also important. In addition, the precision measurement of the quartic coupling of the SM Higgs at a future muon collider can probe a sizable region of the parameter space in LRSM, which can have strong FOPT and observable GW signals, as exemplified in Fig. 9. Table 5. Physical Higgs states and their masses when v L κ 2 κ 1 v R [142]. Here ξ = κ 2 /κ 1 = v EW /v R κ 1 /v R . h is the SM Higgs field. physical states mass squared The corresponding mass matrix elements can be found e.g. in Ref. [135]. In the basis of For the singly charged fields, in the basis of {φ ± 1 , φ ± 2 , δ ± L , δ ± R }, the thermal self-energy is the same as that for the real neutral components, i.e. Π H ± = Π H 0 . For the doubly-charged scalars, in the basis of {∆ ±± L , ∆ ±± R }, the corresponding self-energy is given by For the neutral gauge bosons, in the basis of {W 3 L , W 3 R , B}, the self-energy matrix reads while for the singly-charged gauge bosons, in the basis of {W ± L , W ± R }, the self-energy matrix is

B Conditions for vacuum stability and correct vacuum
The sufficient but not necessary conditions for vacuum stability and correct vacuum in the LRSM are worked out in [103] and listed below (simple analytic formula can only be obtained in the condition α 2 = 0): where the condition structure "p ⇐ q" means p needs to be checked if and only if the condition q is true. In this paper, we have chosen λ 2,3,4 = 0, which corresponds to the case of η = σ = 0.