1.
There is a lower bound in the distribution of \(m_{\mathcal {S}}\) , i.e., \(m_{\mathcal {S}}\gtrsim 500\)  GeV. This is because in order to trigger a phase transition, the local minimum at the broken phase \((v(T),~w(T))\)  should be
larger than the local minimum at the symmetric phase \((0,~w_0(T))\)  when the temperature is higher than critical temperature. To see this, we show in the upper left plot of Fig. REF  the potential difference \(\Delta V_{bs}(T)\equiv V_{\rm eff}(v(T),~w(T),~T) - V_{\rm eff}(0,~w_0(T),~T)\) 
at temperature \(T=300\)  GeV (the distribution of critical temperature can be found in Fig. REF ) as a function of \(m_{\mathcal {S}}\) , with \(\theta =0.4\)  and \(m_{\chi }=300\)  GeV. The blue, green, and red points in Fig. REF  represent the results with \(w=300\)  GeV, 400 GeV, and 500 GeV, respectively. As shown in the drawing, the potential difference slowly increases with \(m_{\mathcal {S}}\)  and remains
negative for \(m_{\mathcal {S}}\lesssim 500\)  GeV. At larger values of \(m_{\mathcal {S}}\) , the potential difference sharply increases to a positive value.
Therefore, the lower bound \(m_{\mathcal {S}}\gtrsim 500\)  GeV is to guarantee a successful phase transition.
2.
There is a lower bound in the distribution of \(\left| \theta \right|\) , i.e., \(\left| \theta \right|\gtrsim 0.2\) .
Such a relatively large lower bound on the absolute value of the mixing angle is straightforwardly derived from the requirement of a sufficiently strong first-order EWPT.
As shown in Eq. (), the mixing angle \(\left| \theta \right|\)  directly controls the Higgs-portal interacting strength \(\left| \lambda _m \right|\)  (whose distribution can be found in the left plot of the second row of Fig. REF ).
For \(\left| \theta \right|\lesssim 0.2\)  (or, equally, \(\left| \lambda _m \right|\lesssim 1\) ),
the Higgs boson \(h\)  and singlet scalar \(s\)  have nearly no mixing and the phase transition takes place mostly along the \(h\)  direction.
A larger \(|\theta |\)  would induce a significant mixing between the \(h\)  and \(s\)  fields and curve the tunneling path to be along a linear combination of the two fields (see Fig. REF ) and finally lead to a strong EWPT (see also Ref. [1]} for similar conclusions).
3.
The lower bound on \(w\)  and the decrease in the distribution of \(\left| \theta \right|\)  at larger values are from the bound of unitarity (last condition of Eq. (REF )).
We plot this bound condition in various planes in Fig. REF , the colored regions are excluded by the constraints.
In the upper right plot of Fig. REF , we show the unitarity bounds in the \(\theta \) -\(m_{\mathcal {S}}\)  plane, fixing \(m_{\chi }=300\)  GeV.
The values of \(w\)  in this plot are 100 GeV (green), 200 GeV (yellow), 300 GeV (blue), 400 GeV (light blue),
600 GeV (pink), and 1000 GeV (sandy-brown), respectively. As shown in this plot, when \(w=100\)  GeV, the heavy scalar mass \(m_{\mathcal {S}}\)  is restricted to lie below about 500 GeV. As one increases \(w\) , more parameter space opens up.
Sufficient sample points for strong EWPT are available for \(w\gtrsim 200\)  GeV.
These indicate that the unitarity bound on large values of \(m_{\mathcal {S}}\)  can be avoided by increasing \(w\) .
In the lower left plot of Fig. REF , we plot the unitarity bounds in the \(w\) -\(m_{\mathcal {S}}\)  plane,
with \(\theta =0.2\)  (blue), 0.4 (pink), and 0.6 (green), respectively. As shown in this plot, the available parameter space
is largely reduced with the increase mixing angle, which explains the decreasing behavior in the \(\left| \theta \right|\)  distribution at large values. When \(m_{\mathcal {S}}\gtrsim 500\)  GeV, nearly all of the parameter space is excluded by the unitarity bound for \(w\lesssim 200\)  GeV. The constraints become independent of \(m_{\mathcal {S}}\)  for \(w\gtrsim 500\)  GeV and \(\theta \gtrsim 0.4\) .
4.
The distribution of pseudo-Goldstone DM mass \(m_{\chi }\)  has a peak around \((200-400)\)  GeV and then decreases at larger values. The decreasing behavior of \(m_{\chi }\)  distribution involves a few contributions. One is the unitarity bound; it is easy to check that the unitarity bound becomes tighter for a larger DM mass \(m_{\chi }\) . Another reason is the condition (REF ), which is to ensure that the stationary point located along the \(h\)  axis cannot be a local minimum. We show this constraints in the lower right plot of Fig. REF ,
the colored regions are excluded by the constraints. In this plot, the value of \(w\)  is fixed at 1000 GeV,
the green, yellow, pink, and blue regions represent the bounds with \(m_{\mathcal {S}}=\)  500 GeV, 700 GeV, 1000 GeV, and 1500 GeV, respectively.
We see that for \(\theta \gtrsim 0.2\) , the parameter space with \(m_{\chi }\gtrsim 1000\)  GeV is excluded. However, we should note that this bound is imposed only when there is a stationary point along \(h\)  axis, which requires \(\mu _h^2>0\) .
From the distribution of \(\mu _h^2\)  shown in the upper left plot of Fig. REF , we find that only about half of the total sample points have a stationary point along \(h\)  axis. We also note that such a requirement is somehow too stringent because of the loop corrections on the potential.
Further constraint may arise from the requirement of \(d\Delta V_{bs}(T)/dT^2>0\)  near the critical temperature [2]}, we find this bound is relatively weak and do not give a detailed account here.
5.
We see that within the parameter ranges considered here, the distributions of \(w\)  and \(m_{\mathcal {S}}\)  do not decrease with increasing values. Such distributions are consistent with a recent global fit performed in Ref. [3]}.

<FIGURE>