-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 [1]}. 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 [2]}. 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 [3]}
\(\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 [4]}
\(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 [5]}.
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 [6]} (BBO [7]}).
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 [8]}.
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 [9]}, 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 ).\) 
