A-0.1Two different local minima
<FIGURE>The tree-level potential has nine stationary points, among which there are four different points as shown in Fig. REF  marked by green, red, blue, and magenta if a \(\mathbb {Z}_2 \times \mathbb {Z}_2\)  symmetry is considered. The green point lies at the origin, and the red and magenta points are along \(h_1\)  and \(h_2\)  axes, respectively. From the Hesse matrix
\({\cal H} = \begin{pmatrix}-\mu _1^2 + 3 \lambda _{1} h_1^2 + \lambda _{12} h_2^2/2 & \lambda _{12} h_1 h_2 \\\lambda _{12} h_1 h_2 & -\mu _2^2 + 3 \lambda _{2} h_2^2 + \lambda _{12} h_1^2/2\end{pmatrix},\) 

one can find that the only possibility for two local minima happens when the blue and green points are saddle points. The definiteness of the Hesse matrix requires
\(\lambda _{12}\mu _1^2 > 2\lambda _1 \mu _2^2 > 0, \quad \lambda _{12}\mu _2^2 > 2\lambda _2 \mu _1^2 > 0.\) 

Combining these two inequalities one can put a constraint on \(\lambda _{12}\)  as
\(\lambda _{12} > 2 \sqrt{\lambda _1 \lambda _2}.\) 
A-0.2The global minimum
As is discussed at the beginning of this subsection, the global minimum should lie at the SM-Higgs axis, with \(h_1 = 246 {\text{GeV}}\) , which requires
\(\frac{\mu _1^4}{\lambda _1} > \frac{\mu _2^4}{\lambda _2}.\) 
A-0.3\(h_2\)  vacuum appears first when cooling down
As is indicated in Fig. REF -(b), a two-step phase transition requires that a local minimum along the \(h_2\)  axis appears first when cooling down.
Along each axis (\(h_2 = 0\)  or \(h_1 = 0\) ), it could be found that whether there is another extremum besides the origin depends on the quadratic term, \(-\mu _i^2 + D_i T^2\) . At high temperatures, all coefficients in the potential are positive, and the global minimum lies at the origin. When the temperature drops and the quadratic term becomes negative, the local minimum changes to \(h_i^2 \sim (-\mu _i^2 + D_i T^2) / \lambda _i\)  along each axis. We could define the turning temperature \(T_i\)  at which the quadratic coefficient changes sign, given by \(T_i^2 = \mu _i^2 / D_i\) .
Therefore, approximately the criterion for a two-step phase transition can be presented as
\(T_2 > T_1,~ \text{i.e. } \frac{\mu _2^2}{D_2} > \frac{\mu _1^2}{D_1},\) 

where \(D_{1,2}\)  are defined in Eq. (REF ).
<FIGURE><FIGURE>We plot the parameter region that satisfies the three conditions above in Fig. REF . The region enclosed by blue, yellow, and green dashed lines corresponds to each condition, as shown in the legends, and the red one is the overlap region. It could be found that \(\lambda _{12}\)  is allowed in a narrow band in each plot, scaling as \(m_H^2\) . The upper bounds for \(\lambda _{12}\)  are set by the `global minimum' condition. In contrast, the lower bounds from the `local minimum' and \(T_2 > T_1\)  conditions are quite close, indicating that the phase transition goes through two steps given the condition that there are two local minima.
We find that the two-step phase transition happens in a narrow region of the \(\lambda _{12}\)  couplings. This coupling should be moderately large and must be larger as the mass of the \(\Phi _2\)  increases. Furthermore, we find that as the isospin increases, the required \(\lambda _{12}\)  coupling should be larger to realize a two-step phase transition.
Fig. REF  shows how the potential evolves when the temperature drops. The numerical analysis of the effective potential includes the effect of tree-level, one-loop, and Daisy-resummed parts. The left column presents the case that \(\lambda _2 = \lambda _{12} = 2.0,~ m_H = 200 ~{\text{GeV}}\)  and \((I_2, Y_2) = (1/2, 1/2)\) . This set of parameters satisfies three conditions and does realize a two-step phase transition. By contrast, the right one only goes through one step, with the parameter \(\lambda _2 = 0.25,~ \lambda _{12} = 3.02,~ m_H = 300 ~{\text{GeV}}\)  and \((I_2, Y_2) = (1, 1)\)  lying outside the allowed region, even when there are two local minima at zero temperature.
Potential Barriers for Two-step Transition
In the following, we will discuss the possibility of two strong first-order phase transitions (SFOPTs) in Fig. REF -(b) by the criterion that \(v(T_C)/T_C \gtrsim 1\) , where \(T_C\)  denotes the critical temperature at which the potential values at two local minima degenerate. To evaluate the criterion in the first step (green \(\rightarrow \)  magenta), we could analyze the potential along the \(h_2\)  axis. However, Eq. (REF ) cannot generate a barrier separating the origin and another local minimum because there only exists a quadratic and a quartic term. It is necessary to add the effects of the next-to-leading order, which is roughly
\(V_{\text{eff}}\left(0, h_2, T\right) \sim \frac{1}{2} (-\mu _2^2 + D_2 T^2) h_2^2 - E_2 T h_2^3 + \frac{1}{4} \lambda _2 h_2^4.\) 

The cubic term \(E_2\)  comes from the thermal loops of bosons, and the contributions from gauge bosons and scalars are
\(E_2^g &= \frac{1}{16\pi }g^3 I_{W}^{3} + \frac{1}{4\pi }(g^2+g^{\prime 2})^{3/2}Y_{2}^3,\nonumber \\E_2^\Phi &\sim \frac{2I_2 + 1}{6\pi }\sum _{m} b_m \times {\cal O}(1) \sim \frac{2I_2 + 1}{6\pi } \lambda _2.\) 

Here the explicit expression of \(E_2^\Phi \)  requires a specific model, and we only give a rough result with the approximation that \(\lambda _2 \sim \sum _m b_m F(I_2, Y_2, m) \sim \sum _m b_m\) .
The criterion for an SFOPT shows that
\(\frac{v_{2}(T_{C,1})}{T_{C,1}}=\frac{2E_2}{\lambda _2} \gtrsim 1 \quad \Rightarrow \quad \left(1 - \frac{2I_2 + 1}{3\pi }\right)\lambda _2 \lesssim 2E_2^g,\) 

where \(v_2\)  is the value of the local minimum and \(T_{C,1}\)  stands for the critical temperature at which \(V_{\text{eff}}(0, 0, T_{C,1})=V_{\text{eff}}(0, v_{2}(T_{C,1}), T_{C,1})\) . This estimation sets an upper bound for \(\lambda _2\)  under low representations \(N \le 9\) , and for \(N \le 4\) , which indicates that an SFOPT in the first step may require a high representation of \(SU(2)_L\) .
In Sec.  we will give a numerical analysis of the doublet model.
<FIGURE>For the second step of the EWPT (i.e., the path of magenta \( \rightarrow \)  red) in Fig. REF -(b), we choose to reparametrize the field \(h_1\)  and \(h_2\)  by a polar coordinate \((z, \gamma )\)  as is shown in Fig. REF . The center of the polar coordinate is set to be the local minimum along \(h_2\)  axis. The relation between two coordinates are given by \(h_1 = z \cos \gamma \)  and \(h_2 = v_2(T) + z \sin \gamma \) , and the coordinates of two minima are called \((0, 0)\)  and \((\tilde{z}, \tilde{\gamma })\) . Since a barrier already exists between two minima due to condition Eq. (REF ), it is sufficient to use the leading-order potential Eq. (REF ) to describe the phase transition and therefore
\(v_1(T)^2 = \frac{ \mu _1^2 - D_1T^2}{\lambda _1},~ v_2(T)^2 = \frac{\mu _2^2 - D_2T^2}{\lambda _2}, \nonumber \\\tilde{z}(T)^2 = v_1(T)^2 + ~v_2(T)^2,~ \tan \tilde{\gamma }(T) = -\frac{v_2(T)}{v_1(T)}.\) 

Under the polar coordinates, the effective potential Eq. (REF ) is given by
\(V_{\text{eff}}\left(z , \gamma , T\right) &= C_0+ C_2 z^2 +C_3 z^3 +C_4z^4,\) 

where
\(C_0 & = -\frac{1}{4}\lambda _2v_2(T)^4, \quad C_2 = \lambda _2v_2(T)^2 s_\gamma ^2 + \frac{1}{4}(-2\lambda _1v^2 + 2D_1T^2 + \lambda _{12}v_2(T)^2) c_\gamma ^2, \nonumber \\C_3 & = \frac{1}{2} v_2(T) s_\gamma (2 \lambda _{2}s_\gamma ^2 + \lambda _{12} c_\gamma ^2), \quad C_4 = \frac{1}{4}(\lambda _2s_\gamma ^4 + \lambda _{12}s_\gamma ^2 c_\gamma ^2 + \lambda _1 c_\gamma ^4).\) 

Here \(\sin \gamma \)  and \(\cos \gamma \)  is denoted \(s_\gamma \)  and \(c_\gamma \) , respectively. At the critical temperature \(T_{C,2}\)  the values of the potential at two local minima degenerate \(V_{\text{eff}}(0, 0, T_{C,2})=V_{\text{eff}}(\tilde{z}(T_{C,2}), \tilde{\gamma }(T_{C,2}), T_{C,2})\) , which yields that
\(T_{C,2}^2 = \frac{\mu _2^2 - \sqrt{\lambda _2/\lambda _1}\mu _1^2}{D_2 - \sqrt{\lambda _2/\lambda _1}D_1}.\) 

Thus the criterion becomes
\(\frac{v_1(T_{C,2})^2}{T_{C,2}^2} = \frac{v^2}{T_{C,2}^2} - D_1 \sim \frac{\sqrt{\lambda _2/\lambda _1}D_1 - D_2}{\sqrt{\lambda _1\lambda _2} - \mu _2^2/v^2} > \frac{\sqrt{\lambda _2/\lambda _1}D_1 - D_2}{\sqrt{\lambda _1\lambda _2} - 2\lambda _1\lambda _2/\lambda _{12}},\) 

where we use the conditions Eq. (REF ) and (REF ) that \(\sqrt{\lambda _1\lambda _2} > \mu _2^2/v^2 > 2\lambda _1\lambda _2/\lambda _{12}\)  and \(\sqrt{\lambda _1} \sim \sqrt{1/8} \ll \sqrt{1/\lambda _1}\) . The third condition Eq. (REF ) guarantees that \(\sqrt{\lambda _2/\lambda _1}D_1 > D_2\) . Under the approximation \(\lambda _2 \sim \lambda _{12} \sim 1\)  as well as \(D_1 \sim g^2/4 + y_t^2/4 + \lambda _2I_2/3\)  and \(D_2 \sim g^2I_2^2 + \lambda _{12}I_2/3\) , we could estimate the value of the ratio to be
\(\frac{v_1(T_{C,2})^2}{T_{C,2}^2} \gtrsim \frac{g^2(\sqrt{2}/2 - I_2^2) + y_t^2\sqrt{2}/2 + (2\sqrt{2}\lambda _2 - \lambda _{12})I_2/3}{(\sqrt{2\lambda _2} - \lambda _2/\lambda _{12})/4} \sim \frac{-0.45I_2^2 + 0.6I_2 + 1}{0.1} > 1\) 

for at least \(I_2 \le 2\) . It indicates that the second step of the phase transition in Fig. REF -(b) is probable to be a SFOPT. The Fig. REF , which shows the parameter region satisfying the criterion, verifies the validation of our approximation. In this section we only utilize the leading-order effective potential to describe the phase transition, and a complete numerical analysis will be present in the following with specific models.
<FIGURE>
Gravitational Wave Spectrum
Before doing numerical analysis, we briefly discuss the GW produced by the first-order phase transition.
Typically, the GW spectrum from the first-order phase transition could be characterized by phase transition parameters: \(T_t\) , \(\alpha \) , \(\beta /H\) , and \(v_b\) .
The first parameter \(T_t\)  is the transition temperature that is defined by an equation that a bubble nucleation probability \(\Gamma \)  per Hubble volume per Hubble time reaches the unit [1]}:
\(\Gamma /H^4|_{T=T_t}=1,\) 

where
\(\Gamma \simeq {\rm max}\left[T^4\left(\frac{S_3}{2\pi T} \right)^{3/2} \exp (-S_3/T),\,\, R_0^{-4}\left(\frac{S_4}{2\pi } \right)^{2} \exp (-S_4)\right].\) 

Here, \(S_3\)  (\(S_4\) ) is a 3- (4-) dimensional Euclidean action for O(3) (O(4)) symmetric bounce solution.
Also, \(R_0\)  is the size of the nucleating bubble.
The Eq. (REF ) means one critical bubble nucleates into the causality.
The second parameter \(\alpha \)  is the normalized released energy density by the radiation energy density at \(T_t\) :
\(\alpha \equiv \epsilon (T_t)/ \rho _{rad}(T_t),\) 

where \(\rho _{rad}=(\pi ^2/30)g_* T^4\)  and \(\epsilon \)  is released energy density:
\(\epsilon (T)= \Delta V_{eff} -T\frac{\partial \Delta V_{eff} }{\partial T},\quad \Delta V_{eff} = V_{eff}(\varphi _-(T),T) - V_{eff}(\varphi _+(T),T),\) 

where \(\varphi _{+(-)}\)  is \(\varphi \)  at the broken (unbroken) phase.
The third parameter \(\beta /H\)  is the normalized \(\beta \)  parameter which is related to the time variation scale of bubble nucleation rate \(\Gamma (t)=\Gamma _0\exp (\beta t)\) .
This normalized \(\beta \)  parameter is defined as
\(\frac{\beta }{H}\equiv \left.\frac{dS_3}{dT}\right|_{T=T_t}.\) 

The last one is the bubble wall velocity \(v_b\) .
In our following analysis, we fix it by hand for simplicity.
By using these parameters, the GW spectrum from first-order phase transition could be estimated.
There are three sources of GW spectrum: bubble collision [2]}, plasma turbulence [3]}, and compression wave of plasma [4]}.
The spectra could be calculated by the details of bubble dynamics, and the fitting formulas of them are provided in Ref. [4]}.
In our analysis, we will focus on the GW spectrum from compression wave of plasma, which typically is the largest contribution among them.
The formula from the compression wave of plasma [4]} is given by
\(\Omega _{\rm comp} (f) h^2=\widetilde{\Omega }_{\rm comp} h^2 \times (f/\tilde{f}_{\rm comp})^3\left(\frac{7}{4+3(f/\tilde{f}_{\rm comp})^2}\right)^{7/2},\) 

where the tilde parameters are the peak of the spectrum
\(\widetilde{\Omega }_{\rm comp} h^2\simeq 2.65 \times 10^{-6} v_b\tilde{\beta }^{-1}\left(\frac{\kappa _v \alpha }{1+\alpha }\right)^2\left(\frac{100}{g_\ast }\right)^{1/3},\) 

and the peak frequency
\(\tilde{f}_{\rm comp} \simeq 1.9 \times 10^{-5}~{\rm Hz} \frac{1}{v_b}\tilde{\beta }\left(\frac{T_t}{100~{\rm GeV}}\right)\left(\frac{g_\ast }{100}\right)^{1/6}.\) 

Here, the \(g_\ast \)  is the number of degrees of freedom in the plasma.
The parameter \(\kappa _v\)  in Eq. (REF ) is an efficiency factor, which denotes the fraction of the vacuum energy transformed into the bulk motion of the plasma fluid [7]}:
\(\kappa _v(v_b, \alpha )\simeq {\left\lbrace \begin{array}{ll}& \frac{ c_s^{11/5}\kappa _A \kappa _B }{(c_s^{11/5} - v_b^{11/5} )\kappa _B+ v_b c_s^{6/5} \kappa _A}\quad (\text{for}~ v_b \lesssim c_s ) \\& \kappa _B + ( v_b - c_s) \delta \kappa + \frac{( v_b - c_s)^3}{ (v_J - c_s)^3} [ \kappa _C - \kappa _B -(v_J - c_s)\delta \kappa ]\quad (\text{for}~ c_s < v_b < v_J) \\& \frac{ (v_J - 1)^3 v_J^{5/2} v_b^{-5/2}\kappa _C \kappa _D }{[( v_J -1)^3 - ( v_b-1)^3] v_J^{5/2} \kappa _C+ ( v_b - 1)^3 \kappa _D }\quad ( \text{for}~ v_J \lesssim v_b )\end{array}\right.},\) 

where \(c_s\)  is the velocity of sound (\(c_s=0.577\) ) and
\(\kappa _A&\simeq v_b^{6/5} \frac{6.9 \alpha }{1.36 - 0.037 \sqrt{\alpha } + \alpha },\quad \kappa _B\simeq \frac{\alpha ^{2/5}}{0.017+ (0.997 + \alpha )^{2/5} },\nonumber \\\kappa _C&\simeq \frac{\sqrt{\alpha }}{0.135 + \sqrt{0.98 + \alpha }},\quad \kappa _D\simeq \frac{\alpha }{0.73 + 0.083 \sqrt{\alpha } + \alpha }.\) 
In order to discuss the testability of the parameter region for the multi-step EWPT, we estimate the signal-to-noise (S/N) ratio for measurements of GW spectrum [8]}, [4]}.
The S/N ratio is defined as
\(S/N\,\, {\rm ratio} \equiv \sqrt{\delta \times t_{obs} \int ^{f_{max}}_{f_{min}}df \left[\frac{h^2 \Omega _{GW}(f)}{h^2\Omega _{Sens}(f)}\right]^2},\) 

where \(\delta \)  is the number of independent channels for the experiment, \(t_{obs}\)  corresponds to the experimental observation period, \(h^2 \Omega _{GW}(f)\)  denotes the GW signal coming from the first-order phase transition and \(h^2\Omega _{Sens}(f)\)  denotes the sensitivity of experiments, such as LISA [10]}, DECIGO and BBO [11]}.
When the S/N ratio at the LISA experiment is larger than 10, we could detect the GW spectrum [4]}.
In our analysis, we take S/N>10 as a criterion for whether the GW spectrum could be detected or not.
 Supercooling phase transition
For the second step of the EWPT in Fig. REF -(b), a barrier between magenta and red points appears at zero temperature.
Then, the supercooling phase transition may be realized.
The condition of the supercooling phase transition is typically given as \(\alpha \gg 1\) .
In the case of the supercooling phase transition, the broken phase may not fill with causality.
Therefore we should estimate the temperature for bubble percolation for the phase transition to complete successfully.
The percolation temperature \(T_p\)  means the broken phase is filled at least \(34\%\)  of the comoving volume at this temperature [13]}, [14]}.
To evaluate the \(T_p\) , we use the probability of finding a point still in the false vacuum:
\(P(T)=e^{-I(T)},\quad I(T)=\frac{4\pi }{3}\int ^T_{T_c} dT^{\prime } \Gamma (T^{\prime }) a(T^{\prime })^3 r(T, T^{\prime })^3,\) 

where \(a(T^{\prime })\)  is the Friedmann-Robertson-Walker scale factor, and \(r(T,T^{\prime })\)  is the growing comoving size of a bubble from \(T^{\prime }\)  to \(T\) :
\(r(T, T^{\prime })=\int ^t_{t^{\prime }}\frac{d\tilde{t}v_w}{a(\tilde{t})}.\) 

The \(P(T)\)  in Eq. (REF ) denotes the amount of true vacuum volume per unit comoving volume.
By using it, the percolation temperature \(T_p\)  is defined as
\(I(T_p)=0.34.\) 

The physical volume of the false vacuum space \({\cal V}_{\rm false}\)  decreases around percolation to complete the phase transition.
The condition is given by [14]}
\( \frac{1}{{\cal V}_{\rm false}}\frac{d{\cal V}_{\rm false}}{dt}<0\quad \rightarrow \quad H(T)\left(3+ T\frac{dI(T)}{dT}\right)<0.\) 

From this condition and equation, we can obtain the \(T_p\) .
Since the bubbles for the broken phase collide with each other at this temperature, the phase transition parameters are determined at this temperature \(T_p\) .
In the S/N ratio phase transition, the GW spectrum is replaced as [14]}
\(\Omega _{\rm comp} (f) h^2=2.65 \times 10^{-6} v_b\frac{H\bar{R}}{(8\pi )^{1/3}}\left(\frac{\kappa _v \alpha }{1+\alpha }\right)^2\left(\frac{100}{g_\ast }\right)^{1/3}(f/\tilde{f}_{\rm comp})^3\left(\frac{7}{4+3(f/\tilde{f}_{\rm comp})^2}\right)^{7/2},\) 

where
\(\tilde{f}_{\rm comp} \simeq 1.9 \times 10^{-5} \frac{1}{v_b}\frac{(8\pi )^{1/3}}{H\bar{R}}\left(\frac{T_t}{100~{\rm GeV}}\right)\left(\frac{g_\ast }{100}\right)^{1/6}.\) 

where \(\bar{R}\sim (v_w-c_s)R_{\rm max}\)  and \(R_{\rm max}\)  is the size of the bubble of maximum-energy-configuration of bubble distributions at percolation.
We will evaluate the GW spectrum from the first-order phase transition in our model by the above equations.
In the following numerical results, we will obtain the phase transition parameters in the model with the EWPT in Fig. REF -(b) and will discuss the testability of the model parameters by the S/N ratio.
After that, we will show whether the two first-order EWPTs could be realized or not.

Numerical Results
 Patterns of EWPTs
According to the results in Sec. , the EWPT in Fig. REF -(a) and -(b) could be realized in the model, especially, multi-step EWPT in Fig. REF -(b) can occur by the red parameter region in Fig. REF .
We will numerically analyze the detail of the possibility of the pattern of the EWPT in the model with \(Y_{2}=1/2\)  and \(I_{2}=1/2\) .
In the numerical analysis in this section, the four parameters in the potential, \( \mu _1^2\) , \(\mu _2^2\) , \( \lambda _1\)  and \(\lambda _{12}\) , will be taken as \(v\) , \(v_2\) , \(m_h\)  and \(m_H\)  by using the stationary conditions and second derivative of the potential of Eqs. (REF ), (REF ), (REF ) and () in Appendix B.
The conditions for the multi-step EWPT, which are shown in Eqs. (REF ), (REF ) and (REF ), are given by Fig. REF .
<FIGURE>These color regions can satisfy the conditions, including one-loop effects in the cases that \(v_2=\)  100, 150, and 200 GeV.
In these cyan regions, the multi-step EWPT may be generated, however, we use the potential without finite temperature effects to describe these regions.
Therefore, we should consider the potential with finite temperature effects to get more detail on the possibility of the EWPTs, like the condition of "\(h_2\)  vacuum appears first when cooling down".
The results, including the finite temperature effect, are given in Fig. REF .
<FIGURE>To clarify the one-step and multi-step EWPT in Fig. REF -(a) and -(b), we compare the temperatures starting the phase transition along \(h_1\)  and \(h_2\)  axes.
In the red (green) square marks, such a temperature along \(h_2\)  (\(h_1\) ) is higher than \(h_1\)  (\(h_2\) ).
Furthermore, the one-step EWPT of Fig. REF -(a) (green squares) and the second step of the EWPT in Fig. REF -(b) (red squares) are first order.
Since a barrier in the second path of the multi-step EWPT at red squares (one-step EWPT at green squares) occurs by the non-thermal effects (thermal effects), the frequency of GW may be different between the path of the EWPT of Fig. REF -(a) and -(b).
According to the numerical results, the peak frequency at red squares is less than about \(10^{-1}\)  Hz, while the frequency at green squares is \({\cal O} (1)\)  Hz.
Therefore, we may distinguish the patterns of the EWPT of Fig. REF -(a) and -(b) by the GW observation.
At the black square marks, the EWPT cannot be completed at the current universe because the barrier for the second step of the EWPT in Fig. REF -(b) does not disappear at zero temperature.
At the red star mark in the right panel of Fig. REF , both steps of the EWPT in Fig. REF -(b) are first order.
It is difficult to realize such a phase transition because the first step relies on finite temperature effects.
The detail of that will be discussed in the following paragraph of "Possibility of two first-order EWPTs".
As the numerical results in Fig. REF , the multi-step EWPT in Fig. REF -(b) could be realized in the case of the large \(\lambda _2\)  value, that can follow the condition in Eq. (REF ).
Also, when we take heavy \(m_H\)  and large \(\lambda _2\)  value, the EWPT has not been completed, because the height difference between the magenta and blue points \(\Delta V\)  becomes high by heavy \(m_H\)  and the large \(\lambda _2\) : \(\Delta V\sim \frac{\lambda _2}{4 v^4}( 2v_2^2(m_H^2+\lambda _2v_2^2)-m_h^2v^2 )^2\) .
For example, the large \(\lambda _2\)  with \(v_2=\)  200 GeV is prohibited by the condition of height potential in Eq. (REF ).
These behaviors could follow the discussion in Sec. .

 Isospin and hypercharge effects
The isospin and hypercharge of the multiplet also contribute to the patterns of the EWPT.
Firstly, we note that the hypercharge of additional scalar boson \(Y_{2}\)  does not contribute much to the order of step of the multi-step EWPT.
The effect of \(Y_{2}\)  appears in denominator of \(T_2\) , and it is proportional to \((m_Z^2-m_W^2)\) .
Since the masses of these gauge bosons are almost the same, the \(Y_{2}\)  does not much affect the EWPT.
On the other hand, the small \(I_{2}\)  makes the multi-step EWPT because the effect of \(I_{2}^2\)  in the denominator of \(T_2\) , which does not appear in \(T_1\) , can be neglected.
Then, the \(T_2^2>T_{1}^2\)  condition can be realized.
On the other hand, for the large \(I_{2}\) , the \(T_{1}\)  and \(T_2\)  are roughly given by
\(T_{1}^2&\sim \frac{3m_h^2v^2}{2I_{2}(m_{H}^2 + \lambda _2 v_2^2 )} ,\quad T_2^2\sim \frac{\lambda _2v_2^2v^2 + \frac{3I_{2} \lambda _2^2v^2v_2^2}{\pi ^2}\left( \ln \frac{T_2^2\alpha _B}{Q^2} - \frac{3}{2}\right)}{ I_{2}^2m_W^2},\) 

where the second term in the \(T_2\)  is one-loop effect with the \(I_{2}\)  dependence.
One-loop effect in the \(T_{1}\)  could be neglected because such a term becomes small by using the condition of the height of potential in Eq. (REF ): \(v_2^2<m_hv/\sqrt{2}\) .
By using these terms \(T_{1}\)  and \(T_2\) , the ratio of these terms is
\(\frac{T_2^2}{T_{1}^2} \sim \frac{2\lambda _2(m_{H}^2 + \lambda _2 v_2^2 )v_2^2\left(1 + \frac{3I_{2} \lambda _2}{\pi ^2}\left( \ln \frac{T_2^2\alpha _B}{Q^2} - \frac{3}{2}\right)\right)}{ 3I_{2}m_W^2 m_h^2}.\) 

Due to this ratio, the denominator has the isospin effect, however, the \(I_{2}\)  dependence in this ratio disappears when the last term in the numerator is large.
Then the value of \(\lambda _2\)  is an important source to realize the two-step phase transition.
For the heavy \(m_H\)  case, the \(T_2^2/T_1^2\)  ratio is roughly given by
\(\frac{T_2^2}{T_{1}^2} \sim \frac{ \lambda _2 v_2^2 (1+2I_{2}) + \frac{3m_h^2}{2\pi ^2}\left( \ln \frac{T_2^2\alpha _B}{Q^2} - \frac{3}{2}\right) }{m_h^2}.\) 

In this case, the \(I_2\)  is important to realize the two-step phase transition.
The numerical results of the isospin dependence are given by Fig. REF .
<FIGURE>The horizontal axis is the left (right) figure is \(\lambda _2\)  (\(m_H\) ), and the vertical axes are \(I_{2}\) .
In the black squares, the EWPT is multi-step and has not been completed in the current universe.
According to these figures, the multi-step EWPT could be generated by the large \(\lambda _2\) .
Furthermore, the small or large \(I_{2}\)  value is the source of the multi-step EWPT.
These results could follow the above discussion.
From that, we can find that the EWPT can assure the \(Z_2\)  symmetry after the spontaneous EWSB in the model with \(I_{2}<5\) .

-0.4The Collider Constraints
<FIGURE>Although the multiplet carries the \(Z_2\)  symmetry, which forbids the mixing between the neutral scalar component and the SM Higgs boson, there are still constraints on the model parameters. At the LHC, the Drell-Yan process produces a pair of the charged components of the multiple, and the charged scalar subsequently decays into the SM weak gauge bosons and the neutral scalar. Thus the charged scalars have been tightly constrained except that the charged scalars and neutral scalars have degenerate masses [17]}. So in the following, we focus on the multiplet in which all the components have degenerate masses In fact, even the masses are degenerate for each component, there are still mass splittings among components due to radiative correction from the electroweak loop [18]}. However, we neglect such small corrections.. Furthermore, we also assume that the neutral scalar only contributes to part of the relic abundance, so the constraints from the relic density could be avoided.
One unavoidable constraint from the loop correction is the constraints from the \(h\rightarrow \gamma \gamma \)  measurements at the LHC.
The deviation from the SM prediction in the \(h\rightarrow \gamma \gamma \)  is given by
\(\delta R_{h\rightarrow \gamma \gamma } = \frac{\Gamma _{h\rightarrow \gamma \gamma }^{\rm SM+BSM}-\Gamma _{h\rightarrow \gamma \gamma }^{\rm SM}}{\Gamma _{h\rightarrow \gamma \gamma }^{\rm SM}},\) 

where the \(\Gamma _{h\rightarrow \gamma \gamma }^{\rm SM+BSM}\)  is the decay rate of \(h\rightarrow \gamma \gamma \)  in this model, which is given by [19]}
\(\Gamma _{h\rightarrow \gamma \gamma }^{\rm SM+BSM}= \frac{v^2\alpha _e^2m_h^3}{256\pi ^3}\left|\frac{4}{3}A_{1/2}\left(\frac{4m_t^2}{m_h^2}\right)+A_{1}\left(\frac{4m_W^2}{m_h^2}\right)+\sum _{n=-I_{2}}^{I_{2}}\frac{\lambda _{12}v^2}{2m_{\Phi _2}^2}\left(n+Y_{2}\right)A_{1}\left(\frac{4m_{\Phi _2}}{m_h^2}\right)\right|^2\) 

where \(\alpha _e=e^2/4\pi \) the loop function \(A_{1/2}\) , \(A_1\)  and \(A_0\)  is defined as [20]}
\(A_1(\tau ) &= 2 + 3 \tau + 3 \tau (2 -\tau ) f(\tau ),\\A_{\frac{1}{2}}(\tau ) &= -2 \tau \left\lbrace  1 +(1 -\tau ) f(\tau )\right\rbrace ,\\A_0(\tau ) &=\tau \left(1 -\tau f(\tau )\right)\) 

with
\(f(\tau ) = \left\lbrace \begin{array}{ll}\arcsin ^2\sqrt{1/\tau } & \quad \tau \ge 1 \\-\frac{1}{4}\left( \ln \frac{1+\sqrt{1-\tau }}{1-\sqrt{1-\tau }} - i \pi \right) & \quad \tau <1 \\\end{array}\right. .\) 

The current bound for the effective coupling of the Higgs boson to the photon \(\kappa _\gamma \)  is \(1.02^{+0.16}_{-0.14}\)  at 95\(\%\)  confidence level [21]}.
Fig. REF  represents the condition of two-step phase transition and the \(h\gamma \gamma \)  observation in the plane of (\(m_{\Phi _2}\) , \(\lambda _{12}\) ).
\(m_{\Phi _2}\)  is the degenerate mass of additional charged-scalar bosons.
The left and right panels are the results in the model with (\(I_{2}\) , \(Y_{2}\) )=(1/2,1/2) and (1,0), respectively.
The blue region can satisfy the conditions of a two-step phase transition.
The black region is allowed region with \(|\delta R_{h\rightarrow \gamma \gamma }|<0.16\) .
Due to this figure, the two-step phase transition can be realized in the allowed region for \(h\gamma \gamma \)  observation.
The numerical results in Figs. REF  and REF  can be avoided the current bound of \(h\rightarrow \gamma \gamma \)  observation.
 First-order phase transition and GW Signatures
In this part, we will show the testability of the parameter region of the left panel of Fig. REF .
Fig. REF  represents the results of the S/N ratio, which is given by in Eq. (REF ), of the green and red squares in the left panel of Fig. REF .
The left (right) two figures represent the S/N ratio for LISA [10]} (BBO [11]}).
The experimental period for the S/R ratio is chosen as 5 years.
The red dashed line in these figures represents \(S/N=10\) .
The parameter regions above the red dashed line can realize \(S/N>10\) .
Therefore, such a parameter region can be tested by each experiment.
<FIGURE>The top (bottom) two figures have \(v_b\)  = 0.2 (0.5), respectively.
The low speed of bubble wall velocity is required to generate the baryon asymmetry of the universe [24]}.
The upper right parts of the white regions in these figures correspond to the black squares in Fig. REF , which represent that the EWPT has not been completed in the current universe.
The parameter region close to the upper right white region can be tested at the LISA.
On the other hand, the BBO can test the almost parameter region realizing, especially the multi-step EWPT in Fig. REF -(b).
Therefore, we can use the GW observation experiments to explore the possibility of the multi-step EWPT to restore the \(Z_2\)  symmetry after the spontaneous EWSB.
The first-order EWPT in Fig. REF -(a) can be realized by the negative \(\mu _2^2\)  case. However, the first-order phase transition is not strong enough to generate the detectable GW spectrum because the \(\lambda _{12}v^2\)  is typically smaller for negative \(\mu _2^2\)  than for positive \(\mu _2^2\)  to determine the \(m_H\) .

 Possibility of two first-order EWPTs
For the red star point in Fig. REF  with \(v_2=200~{\rm GeV}\) ,
the two peaks of GW spectra could be produced by the two first-order phase transitions in (green \(\rightarrow \)  magenta) and (magenta \(\rightarrow \)  red) paths, depicted in Fig. REF .
Parameters at the red star point are \((v_2, m_H,\lambda _{2}, I_{2}, Y_{2})\)  = (200 GeV, 180 GeV, 0.08, 1/2, 1/2).
The GW spectra of this parameter are shown by Fig. REF -(1).
The red (black) lines are the GW spectrum from the first-step (second-step) EWPT in Fig. REF -(b).
The phase transition parameters are \((T_t, \alpha , \tilde{\beta }, v_b)=\)  (149 GeV, \(3.79\times 10^{-4}, 8.37\times 10^5\) , 0.5) [first step] and (142 GeV, \(4.70\times 10^{-3}, 4.91\times 10^3\) , 0.5) [second step]. The blue and yellow lines correspond to the sensitivity of the LISA and BBO experiments, respectively.
The S/N ratio for these spectra are given as \(S/N\)  (LISA) \(\sim 10^{-16}\) , \(S/N\)  (BBO) \(\sim 10^{-7}\)  [first step] and \(S/N\)  (LISA) \(\sim 10^{-7}\) , \(S/N\)  (BBO) \(\sim 0.565\)  [second step].
Therefore, it could not be easy to detect these GW spectra in Fig. REF -(1).
On the other hand, we take the second benchmark point with larger \(v_2\)  than the star point: \(v_2=223\)  GeV.
<FIGURE>The GW for the second benchmark point is described in Fig. REF -(2), and the phase transition parameters are \((T_t, \alpha , \tilde{\beta }, v_b)=\)  (153 GeV, \(5.03\times 10^{-4}, 9.58\times 10^4\) , 0.5) [first step] and (109 GeV, \(3.06\times 10^{-2}\) , 72.2, 0.5) [second step].
The S/N ratio for each path in Fig. REF -(2) are \(S/N\)  (LISA) \(\sim 10^{-16}\) , \(S/N\)  (BBO) \(\sim 10^{-9}\)  [first step] and \(S/N\)  (LISA) \(\sim \)  58.0, \(S/N\)  (BBO) \(\sim \)  238 [second step].
Therefore, we could detect one peak of GW spectrum at least. However, these experiments are difficult to observe both GW spectra from the two first-order EWPTs in Fig. REF -(b).
The reasons why the GW spectrum from the first step is so small are (i) first-step EWPT does not have non-thermal contributing to cubic term \(h_2^3\) , and (ii) \(\lambda _2\)  could not be small enough to generate the strongly first-order EWPT.
The detectable GW spectrum could require a large value of \(v_C/T_C\) ; in other words, it could require a small value of \(\lambda _2\) .
However, the condition of \(T_2>T_{1}\)  in Eq. (REF ) cannot be satisfied by such a small value of \(\lambda _2\) .
In the case of large \(v_2\)  value, multi-step phase transitions are realized even for small values of \(\lambda _2\) , and two first-order EWPTs may be realized.
But the condition of the height of the potential between the magenta and red points in Fig. REF  becomes strict by the large \(v_2\) : \( \lambda _2 v_2^4 < m_h^2v^2/2\) .
Therefore, the parameter region for producing the detectable two peaks of the GW spectrum is very narrow.
Furthermore, we study the supercooling phase transition (\(\alpha \gg 1\) ).
The benchmark point of the supercooling phase transition is \((v_2, m_H,\lambda _{2}, I_{2}, Y_{2})\)  = (225 GeV, 168.64 GeV, 0.08, 1/2, 1/2).
Then the percolation temperature of the second path is \(T_p\sim 4\)  GeV, and \(\alpha \sim 3880\) .
However, even in the supercooling case, the first step of the EWPT is not strong.
It is difficult to detect both peaks of the GW spectra coming from both steps of the EWPT in Fig. REF -(b); however, at least one peak of the spectra could be detected.
Suppose the multi-step EWPT is unable to restore the \(Z_2\)  symmetry after the spontaneous EWSB, e.g., Fig. REF -(e) phase transition. In that case, we could detect both peaks of GW spectra produced by two first-order EWPTs [25]}, because the condition of the height of the potential is not necessary to realize the multi-step EWPT.

Summary
We have investigated that the two-step phase transition can be naturally realized in the scalar potential with the discrete \(Z_2\)  symmetry due to degenerate extreme along field directions. To lay out the general conditions on such a two-step phase transition, we consider the most general electroweak \(N\) -plet scalar extension. For the first time, we provide the general form of the Deybe masses and obtain the two-step phase transition conditions for general isospin and hypercharge scalars. We found that there is a strong correlation between the coupling \(\lambda _{12}\)  and the heavy scalar mass, and the coupling \(\lambda _{12}\)  needs to be in a moderate region because of the global minimum and \(T_2 > T_1\)  conditions. As the heavy scalar mass increase, the coupling \(\lambda _{12}\)  is required to take a larger value, and similarly for larger isospin multiplet.
Although both steps of the EWPT could be first-order, it is difficult to detect two peaks of GW spectra produced from these first-order EWPTs.
The first-order EWPT along the first step can be realized by thermal loop effects, and then typically, the GW spectrum frequency is \({\cal O} (1)\)  Hz. The detectable GW spectrum from this step can be produced if there is a small \(\lambda _2\)  or large \(v_2\)  value. Since a barrier in the second step of the EWPT occurs by the tree-level effects, the detectable GW could be generated through the second step, of which the frequency is less than about \(10^{-1}\)  Hz. Therefore, we may distinguish different patterns of the EWPT by the GW observation.
Even in the case that the multi-step EWPT is able to restore the discrete symmetry, at least one peak of the spectrum could be detected by future GW observation experiments. According to the numerical results, the second step produces the GW spectrum that could be detected in the BBO experiment. Furthermore, it is possible to explore such parameter region of the moderately large \(\lambda _{12}\)  by the Higgs diphoton data in future hadron and lepton colliders. From these experiments, we would explore the possibility of multi-step phase transition comprehensively.
This work is supported in part by the National Key Research and Development Program of China Grant No. 2020YFC2201501. Q.H.C. is supported in part by the National Science Foundation of China under Grant Nos. 11675002, 11635001, 11725520, and 12235001. J.H.Y. is supported in part by the National Science Foundation of China under Grants No. 12022514, No. 11875003, and No. 12047503, and CAS Project for Young Scientists in Basic Research YSBR-006, and the Key Research Program of the CAS Grant No. XDPB15.

Field-dependent masses in several multiplets
The field-dependent mass (matrix) can be calculated by taking the second derivatives of the potential as
\(M_{ij}^2(h_1, h_2) = \frac{\partial ^2 V_0}{\partial \phi _i \partial \phi _j} \Bigg |_{h_1, h_2},\) 

where \(\phi \)  can be scalars, fermions, and gauge bosons. We list the mass matrix squared under various representations \((I_2, Y_2)\)  as following:
-0.5Notation

NS, CS, and NCPC are abbreviations of neutral scalar, charged scalar, and neutral CP conserving scalar. By default, \(n_{\text{NS}} = n_{\text{NCPC}} = 1,~ n_{\text{CS}} = 2\) .

\(\chi \)  is components of \(\Phi _1\)  except \(h_1\) . By default, \(n_\chi = 3\) .

\(c^+ = 2\operatorname{Re}(c) = c + c^*,~ c^- = 2\operatorname{Im}(c) = -i(c - c^*)\)  for parameter c.

-0.6Real singlet
\((0, 0)\) :
\(M_{\rm {NS}}^2 &= \left(\begin{array}{cc}3\lambda _{1}h_1^2 + a_0h_2^2 - \mu _1^2 & 2a_0h_1h_2 \\2a_0h_1h_2 & 12b_{0}h_2^2 + a_0h_1^2 - \mu _2^2\end{array}\right), \nonumber \\M_{\chi }^2 &= \lambda _{1}h_1^2 + a_0h_2^2 - \mu _1^2.\) 
-0.7Complex singlet
\((0, 0)\) :
\(M_{\rm {NS}}^2 &= \left(\begin{array}{ccc}(M_{\rm {NS}}^{11})^2 & (a_0+a_1^{\prime +})h_1h_2 & -a_0^{\prime -}h_1h_2 \\(a_0+a_1^{\prime +})h_1h_2 & (M_{\rm {NS}}^{22})^2 & -\frac{1}{2}a_0^{\prime -}h_1^2 - 3\left(b_0^{\prime -} + \frac{1}{2}b_0^{\prime \prime -}\right)h_2^2 \\-a_0^{\prime -}h_1h_2 & -\frac{1}{2}a_0^{\prime -}h_1^2 - 3\left(b_0^{\prime -} + \frac{1}{2}b_0^{\prime \prime -}\right)h_2^2 & (M_{\rm {NS}}^{aa}))^2\end{array}\right), \nonumber \\M_{\chi }^2 &= \lambda _{1}h_1^2 + \frac{1}{2}a_0h_2^2 - \mu _1^2,\) 

where
\((M_{\rm {NS}}^{11})^2 & = 3\lambda _{1}h_1^2 + \frac{1}{2}(a_0+a_1^{\prime +})h_2^2 - \mu _1^2, \nonumber \\(M_{\rm {NS}}^{22})^2 & = 3(b_{0}+b_0^{\prime +}+b_0^{\prime \prime +})h_2^2 + \frac{1}{2}(a_0+a_1^{\prime +})h_1^2 - \mu _2^2, \nonumber \\(M_{\rm {NS}}^{aa})^2 & = (b_0 - 3b_0^{\prime +})h_2^2 + \frac{1}{2}(a_0 - a_0^+)h_1^2 - \mu _2^2.\) 
-0.8Doublet
\((1/2, 1/2)\) :
\(M_{\rm {NS}}^2 &= \left(\begin{array}{cccc}(M_{\rm {NS}}^{11})^2 & (a_0 + a_1^{\prime +})h_1h_2 & -\frac{1}{2}a_1^{\prime -}h_2^2 & a_1^{\prime -}h_1h_2 \\(a_0 + a_1^{\prime +})h_1h_2 & (M_{\rm {NS}}^{22})^2 & -a_1^{\prime -}h_1h_2 & \frac{1}{2}a_1^{\prime -}h_1^2 \\-\frac{1}{2}a_1^-h_2^2 & -a_1^-h_1h_2 & (M_{\rm {NS}}^{\phi _0\phi _0})^2 & a_1^{\prime +}h_1h_2 \\a_1^{\prime -}h_1h_2 & \frac{1}{2}a_1^{\prime -}h_1^2 & a_1^{\prime +}h_1h_2 & (M_{\rm {NS}}^{aa})^2\end{array}\right), \nonumber \\M_{\rm {CS}}^2 &= \left(\begin{array}{cc}\lambda _{1}h_1^2 + \frac{1}{2}(a_0 + a_1)h_2^2 - \mu _1^2 & -\frac{1}{2}(a_1 - 2a_1^{\prime })h_1h_2 \nonumber \\-\frac{1}{2}(a_1 - 2a_1^{\prime *})h_1h_2 & b_{0}h_2^2 + \frac{1}{2}(a_0 + a_1)h_1^2 - \mu _2^2\end{array}\right),\) 

where
\((M_{\rm {NS}}^{11})^2 &= 3\lambda _{1}h_1^2 + \frac{1}{2}(a_0 + a_1^{\prime +})h_2^2 - \mu _1^2, \nonumber \\(M_{\rm {NS}}^{22})^2 &= 3b_{0}h_2^2 + \frac{1}{2}(a_0 + a_1^{\prime +})h_1^2 - \mu _2^2, \nonumber \\(M_{\rm {NS}}^{\phi _0\phi _0})^2 &= \lambda _{1}h_1^2 + \frac{1}{2}(a_0 - a_1^{\prime +})h_2^2 - \mu _1^2, \nonumber \\(M_{\rm {NS}}^{aa})^2 &= b_{0}h_2^2 + \frac{1}{2}(a_0 - a_1^{\prime +})h_1^2 - \mu _2^2.\) 
-0.9Real triplet
\((1, 0)\) :
\(M_{\rm {NCPC}}^2 &= \left(\begin{array}{cc}3\lambda _{1}h_1^2 - a_0h_2^2 - \mu _1^2 & -2a_0h_1h_2 \\-2a_0h_1h_2 & 12b_{0}h_2^2 - a_0h_1^2 + \mu _2^2\end{array}\right), \nonumber \\M_{\chi }^2 &= \lambda _{1}h_1^2 - a_0h_2^2 - \mu _1^2, \nonumber \\M_{a}^2 &= 4b_{0}h_2^2 - a_0h_1^2 + \mu _2^2,~ n_{a} = 2.\) 
-0.10Complex triplet
\((1, 0)\) :
\(M_{\rm {NS}}^2 &= \left(\begin{array}{ccc}(M_{\rm {NS}}^{11})^2 & \frac{1}{2}(2a_0 + a_1 - 2a_0^{\prime +})h_1h_2 & a_0^{\prime -}h_1h_2 \\\frac{1}{2}(2a_0 + a_1 - 2a_0^{\prime +})h_1h_2 & (M_{\rm {NS}}^{22})^2 & \frac{1}{2}a_0^{\prime -}h_1^2 - 3(b_0^{\prime -} - b_0^{\prime \prime -})h_2^2 \\a_0^{\prime -}h_1h_2 & \frac{1}{2}a_0^{\prime -}h_1^2 - 3(b_0^{\prime -} - b_0^{\prime \prime -})h_2^2 & (M_{\rm {NS}}^{aa})^2\end{array}\right), \nonumber \\M_{\rm {CS}}^2 &= \left(\begin{array}{ccc}(M_{\rm {CS}}^{G^\pm })^2 & -\frac{1}{2\sqrt{2}}a_1h_1h_2 & -\frac{1}{2\sqrt{2}}a_1h_1h_2 \\-\frac{1}{2\sqrt{2}}a_1h_1h_2 & (M_{\rm {CS}}^{s^+})^2 & a_0^{\prime }h_1^2 - (b_{1} + 2b_{0}^{\prime } - b_{0}^{\prime \prime })h_2^2\\-\frac{1}{2\sqrt{2}}a_1h_1h_2 & a_0^{\prime }h_1^2 - (b_{1} + 2b_{0}^{\prime } - b_{0}^{\prime \prime })h_2^2 & (M_{\rm {CS}}^{s^-})^2\end{array}\right), \nonumber \\M_{G^0}^2 &= \lambda _1h_1^2 + \frac{1}{4}(2a_0 + a_1 - 2a_1^{\prime +})h_2^2 - \mu _1^2,~ n_{G^0} = 1,\) 

where
\((M_{\rm {NS}}^{11})^2 &= 3\lambda _{1}h_1^2 + \frac{1}{4}(2a_0 + a_1 - 2a_0^{\prime +})h_2^2 - \mu _1^2, \nonumber \\(M_{\rm {NS}}^{22})^2 &= 3(b_{0} + b_{1} + b_0^{\prime +} - b_0^{\prime \prime +})h_2^2 + \frac{1}{4}(2a_0 + a_1 - 2a_0^{\prime +})h_1^2 - \mu _2^2, \nonumber \\(M_{\rm {NS}}^{aa})^2 &= (b_{0} + b_{1} - 3b_0^{\prime +})h_2^2 + \frac{1}{4}(2a_0 + a_1 + 2a_0^{\prime +})h_1^2 - \mu _2^2, \nonumber \\(M_{\rm {NS}}^{G^\pm })^2 &= M_{G^0}, \nonumber \\(M_{\rm {NS}}^{s^+})^2 &= \frac{1}{2}(2b_{0} - b_0^{\prime \prime +})h_2^2 + \frac{1}{2}(a_0 + a_1)h_1^2 - \mu _2^2, \nonumber \\(M_{\rm {NS}}^{s^-})^2 &= \frac{1}{2}(2b_{0} - b_0^{\prime \prime +})h_2^2 + \frac{1}{2}a_0h_1^2 - \mu _2^2.\) 
-0.11Complex triplet
\((1, 1)\) :
\(M_{\rm {NCPC}}^2 &= \left(\begin{array}{cc}3\lambda _{1}h_1^2 + \frac{1}{2}a_0h_2^2 - \mu _1^2 & a_0h_1h_2 \\a_0h_1h_2 & 3b_{0}h_2^2 + \frac{1}{2}a_0h_1^2 - \mu _2^2\end{array}\right), \nonumber \\M_{G^0}^2 &= \lambda _{1}h_1^2 + \frac{1}{2}a_0h_2^2 - \mu _1^2,~ n_{G^0} = 1, \nonumber \\M_{a}^2 &= b_{0}h_2^2 + \frac{1}{2}a_0h_1^2 - \mu _2^2,~ n_{a} = 1, \nonumber \\M_{\rm {CS}}^2 &= \left(\begin{array}{cc}\lambda _{1}h_1^2 + \frac{1}{2}(a_0 + a_1)h_2^2 - \mu _1^2 & -\frac{1}{2\sqrt{2}}a_0h_1h_2 \nonumber \\-\frac{1}{2\sqrt{2}}a_0h_1h_2 & b_{0}h_2^2 + \frac{1}{4}(2a_0 + a_1)h_1^2 - \mu _2^2\end{array}\right), \nonumber \\M_{{s^{\pm \pm }}}^2 &= (b_{0} + 2b_{1})h_2^2 + \frac{1}{2}(a_0 + a_1)h_1^2 - \mu _2^2,~ n_{s^{\pm \pm }} = 2.\) 
Stationary conditions and neutral scalar masses 
We modify the tree-level conditions in Fig. REF  by using the effective potential with the loop correction at zero temperature in Eq. (REF ), with the non-diagonal terms neglected.
We take \(Q=v=246\)  GeV and the tadpole conditions are defined as follows:
\(\left.\frac{\partial V_{\rm eff}}{\partial \left\langle \Phi _1 \right\rangle } \right|_{ \left\langle \Phi _1 \right\rangle =v, \left\langle \Phi _2 \right\rangle =0} &=-\mu _1^2+\lambda _1 v^2 +\frac{1}{16\pi ^2} \Bigg [ \frac{6\lambda _1+\lambda _{12}}{4} f_+(m_{\Phi _{1,r}}^2,m_{\Phi _{2,r}}^2) + \frac{1}{4} (6\lambda _1-\lambda _{12}) f_-(m_{\Phi _{1,r}}^2,m_{\Phi _{2,r}}^2) \nonumber \\& \quad +(2(2I_{2}-1)-1) \frac{\lambda _{12} m_{N_{2,r}}^2}{2}\left(\ln \frac{m_{N_{2,r}}^2}{Q^2}-1\right) + \frac{6m_W^4}{ v^2}\left(\ln \frac{m_W^2}{Q^2}-\frac{1}{3}\right) \nonumber \\&\quad + \frac{3m_Z^4}{v^2}\left(\ln \frac{m_Z^2}{Q^2}-\frac{1}{3}\right) - \frac{12m_t^4}{v^2}\left(\ln \frac{m_t^2}{Q^2}-1\right) \Bigg ]= 0,\) 
\(\left.\frac{\partial V_{\rm eff}}{\partial \left\langle \Phi _2 \right\rangle } \right|_{ \left\langle \Phi _1 \right\rangle =0, \left\langle \Phi _2 \right\rangle =v_2} &= -\mu _2^2+\lambda _2 v_2^2 +\frac{1}{16\pi ^2} \Bigg [\frac{6\lambda _2+\lambda _{12}}{4} f_+\left((m_{\Phi _{1,b}})^2,(m_{\Phi _{2,b}})^2\right)) \nonumber \\&\quad +\frac{1}{4}(\lambda _{12} - 6\lambda _2) f_-\left((m_{\Phi _{1,b}})^2,(m_{\Phi _{2,b}})^2\right) +3\frac{\lambda _{12}}{2} (m_{N_{1,b}}^s)^2\left(\ln \frac{(m_{N_{1,b}})^2}{Q^2}-1\right) \nonumber \\& \quad +(2(2I_{2}-1)-1)\lambda _{2} (m_{N_{2,b}})^2\left(\ln \frac{(m_{N_{2,b}})^2}{Q^2}-1\right) \nonumber \\&\quad + \frac{3Y_{2}^4m_Z^4v_2^2}{v^4}\left(\ln \frac{Y_{2}^2v_2^2m_Z^2}{Q^2v^2}-\frac{1}{3}\right) + \frac{6I_W^4m_W^4v_2^2}{ v^4}\left(\ln \frac{I_W^2m_W^2v_2^2}{v^2Q^2}-\frac{1}{3}\right) \Bigg ]= 0,\) 

where
\(f_\pm (m_i^2, m_j^2) &= m_i^2\left(\ln \frac{m_i^2}{Q^2}-1\right) \pm m_j^2\left(\ln \frac{m_j^2}{Q^2}-1\right), \\\Delta m_\Phi ^2 &= m_{\Phi _2}^2 - m_{\Phi _1}^2= \sqrt{\left({\cal M}_{\Phi _1\Phi _1}^2-{\cal M}_{\Phi _2\Phi _2}^2\right)^2+4{\cal M}_{\Phi _1\Phi _2}^4},\\\mathcal {M}_{h, H}^{2}&=\left(\begin{array}{ccc}{\cal M}_{\Phi _1\Phi _1}^2 & {\cal M}_{\Phi _1\Phi _2}^2 \\{\cal M}_{\Phi _2\Phi _1}^2 & {\cal M}_{\Phi _2\Phi _2}^2\end{array}\right),\) 

and
\(m_{\Phi _{1,r}}^2 &= -\mu _1^2 + 3\lambda _1 v^2 , &&m_{\Phi _{2,r}}^2 = -\mu _2^2 + \lambda _{12} v^2/2, \nonumber \\m_{N_{1,r}}^2 &= -\mu _1^2 + \lambda _{1} v^2, &&m_{N_{2,r}}^2 = -\mu _2^2 + \lambda _{12} v^2/2, \nonumber \\m_{\Phi _{1,b}}^2 &= -\mu _1^2 + \lambda _{12} v_2^2/2 , &&m_{\Phi _{2,b}}^2 = -\mu _2^2 + 3\lambda _{2} v_2^2 ,\nonumber \\m_{N_{1,b}}^2 &=-\mu _1^2 + \lambda _{12} v_2^2/2 , &&m_{N_{2,b}}^2 = -\mu _2^2 + \lambda _{2} v_2^2/2.\) 

To replace the Higgs boson mass \(m_h\)  and additional neutral CP-even boson mass \(m_{ H}\) , we use the following second derivative of the effective potentials:
\(m_h^2\equiv \left.\frac{\partial ^2 V_{\rm eff}}{\partial \left\langle \Phi _1 \right\rangle ^2} \right|_{ \left\langle \Phi _1 \right\rangle =v, \left\langle \Phi _2 \right\rangle =0} &=2 \lambda _1v^2 +\frac{v^2}{32\pi ^2} \Bigg [A^r_2\ln \frac{m_{\Phi _{1,r}}^2m_{\Phi _{2,r}}^2}{Q^4} -A^r_3\ln \frac{m_{\Phi _{1,r}}^2}{m_{\Phi _{2,r}}^2}\nonumber \\&\quad + 12\lambda _1^2 \ln \frac{m_{N_{1,r}}^2}{Q^2}+ 4\lambda _{12}^2(2(2I_{2}-1)-1) \ln \frac{m_{N_{2,r}}^2}{Q^2} + 12 \frac{m_Z^4}{v^4}\left(\ln \frac{m_Z^2}{Q^2}+\frac{2}{3}\right)\nonumber \\&\quad +24 \frac{m_W^4}{ v^4}\left(\ln \frac{m_W^2}{Q^2}+\frac{2}{3}\right) -48 \frac{m_t^4}{ v^4}\ln \frac{m_t^2}{Q^2} \Bigg ] \\m_{ H}^2 \equiv \left.\frac{\partial ^2 V_{\rm eff}}{\partial \left\langle \Phi _2 \right\rangle ^2} \right|_{ \left\langle \Phi _1 \right\rangle =v, \left\langle \Phi _2 \right\rangle =0} &= -\mu _2^2 + \lambda _{12}v^2 +\frac{1}{64\pi ^2} \Bigg [ B_1f_-(m_{\Phi _{2,r}}^2,m_{\Phi _{1,r}}^2)+ B_2f_+(m_{\Phi _{2,r}}^2,m_{\Phi _{1,r}}^2) \nonumber \\&\quad + 4(2(2I_{2}-1)-1) \lambda _2 m_{N_{2,r}}^2 \left(\ln \frac{m_{N_{2,r}}^2}{Q^2} - 1\right) + 12 \frac{m_Z^4Y_{2}^2}{v^2}\left(\ln \frac{m_Z^2}{Q^2}-\frac{1}{3}\right) \nonumber \\&\quad +24 \frac{m_W^4I_W^2}{ v^2}\left(\ln \frac{m_W^2}{Q^2}-\frac{1}{3}\right) \Bigg ] ,\) 

where
\(A^r_2 &=\frac{1}{4}\left( (6\lambda _1+\lambda _{12})^2 +(6\lambda _1-\lambda _{12})^2\right), \quad A^r_3 = -\frac{36\lambda _1^2-\lambda _{12}^2}{2}\\B_1&= \frac{2}{\Delta m^2}\left[\left(-6\lambda _2+\lambda _{12}\right) \left({\cal M}_{\Phi _1\Phi _1}^2-{\cal M}_{\Phi _2\Phi _2}^2\right)+ 4\lambda _{12}^2 v^2 \right]\\B_2 &= 6\lambda _2+\lambda _{12} .\) 

The model parameters can be replaced by Eqs. (REF ), (REF ) and (REF ):
\((\mu _1^2, \mu _2^2, \lambda _1, \lambda _2, \lambda _{12} )\rightarrow (v, v_2 , m_h , \lambda _2, m_{ \Phi _2}^2 ).\) 
