\section{Introduction}

The LIGO/Virgo events GW170817~\cite{TheLIGOScientific:2017qsa,Abbott:2018wiz} and GW190425~\cite{Abbott:2020uma} are in general considered as binaries of neutron stars (BNSs), as~well as the possibility of being binaries of hybrid stars. However, in~\cite{Abbott:2020khf}, a~new binary coalescence event GW190814 is reported, while one is a black hole with 23 $M_{\odot}$, and~the other is a compact object with only 2.6 $M_{\odot}$ (with a range between 2.50 and 2.67$M_{\odot}$). Then, the~question is what is the identity of the secondary compact object? There are several possibilities: It could be the lightest black hole, or~the heaviest neutron star, or~a dark star. In~addition, there are also possibilities for this compact object to be a hybrid star made of neutron and dark matter~\cite{Zhang:2020pfh}.  Below, we will extend a few more discussions on the above possibilities, and~then focus on elaborating more on the scenarios of dark and hybrid stars in the rest of the~paper.

{ {The scenario of black hole---
}The black holes observed by LIGO/Virgo are all known to have masses more than 5 $M_{\odot}$. However, the~remnant of GW170817, which is identified to be the BNS by parameter estimation (PE), is estimated to have a mass of about 2.7 $M_{\odot}$, and~is believed to be a black hole~\cite{Abbott:2018wiz}. Thus, the~black hole of mass less than 5 $M_{\odot}$ can be formed via the coalescence of BNS.
Since no information about tidal Love number is given in the released PE of GW190814, and~no electromagnetic follow-up is reported, we cannot exclude the possibility that this 2.6 $M_{\odot}$ object is the lightest black hole, which could likely be either formed astronomically via the coalescence of BNS, or~formed in the early universe as a primordial black hole~\cite{Yuan:2020iwf}. Alternatively, the~light black hole may come from core-collapse supernova of a massive normal star~\cite{Belczynski:2011bn}.

{The scenario of neutron star---}Before this GW190814 event, masses of neutron stars are generally considered to be less than 2.3 $M_{\odot}$ \cite{Margalit:2017dij, Rezzolla:2017aly}. In~fact, within~our galaxy, the~most massive known pulsar has $2.14^{+0.10}_{-0:09} M_{\odot}$ at $68.3\%$ credible interval measured by the Shapiro time-delay of its white dwarf companion~\cite{Cromartie:2019kug}. Most known equations of state (EoSs) for the dense nuclear matter used in the astronomical search of neutron stars cannot sustain such high mass as 2.3 $M_{\odot}$, although~this can also be seen as a result of tuning the parameters of the EoSs to not violate the upper mass limit of the astronomical observations.
On the other hand, rapid spinning of neutron star~\cite{Breu:2016ufb,Most:2020bba} can increase the maximal mass up to $\sim$20\%, which will lead to $\sim$2.7 $M_{\odot}$. This case cannot be totally ruled out since there is no constraint on the spin of the secondary object in GW190814. However, such high spin NSBH system may hardly exist because of the expected subsequent collapse into black hole due to the loss of spin via gravitational wave~radiation.

{The scenario of dark star and hybrid star---}Currently, dark stars have not yet been confirmed in observations so that it remains to pin down the accessible model space of dark matter and the associated EoSs by the future observations in the coming era of gravitational wave astronomy. In~light of the varieties of dark matter models, it is easy for the resultant dark stars to cover a wide mass range, say 1 to 5 $M_{\odot}$, or~even one order higher with proper parameters. This makes the dark stars or the hybrid stars of dark and nuclear matters to be highly possible candidates to explain the companion star of GW190814. In~the remaining of this paper, we will explore this possibility by studying the mass--radius relations for the dark and hybrid stars based on various dark matter models inspired by the particle physics. Specifically, we consider the massive bosonic field $\phi$ with the following self-interactions, $\phi^4$, $\phi^6$, their linear combinations and $\phi^{10}$.

One way to characterize a compact star is through its Tidal Love Number (TLN)\footnote{The TLN denoted by $\Lambda$ is defined by
$
Q_{ab}=- M^5  \Lambda \; {\cal E}_{ab},
$
 where $M$ is the mass of the star, $Q_{ab}$ is the induced quadrupole moment, and~${\cal E}_{ab}$ is the external gravitational tidal field strength.} \cite{Hinderer:2007mb,Postnikov:2010yn}, which~shows the star’s tendency to deform under an external quadrupolar tidal field.
It is known that the TLN of black holes in Einstein gravity is vanishing, and~the overall TLN effect for a compact binary coalescence event is a weighted average of the individual TLNs. Thus, for~GW190814, the~overall tidal effect is insignificant in the resultant gravitational waveform due to the high mass ratio between black hole and the companion star. For~example, even if the 2.6 $M_{\odot}$ object has a large TLN such as 30,000, when combining with zero TLN of the companion black hole of 23 $M_{\odot}$, the~overall TLN $\tilde \Lambda$\footnote{$\tilde{\Lambda}={16\over 13}{(M_1+12 M_2) M_1^4 \Lambda_1+(M_2+12 M_1) M_2^4 \Lambda_2 \over (M_1+M_2)^5}$ where $M_i$ and $\Lambda_i$ with $i=1$ or $2$ is the mass and TLN for the $i$-th component compact object.} contributed to the gravitational waveform is about 43, which can hardly be observed. Since there is no available information on TLN from this event anyway, we need more future events with much smaller masses where the TLN judgement can be applied. For~the above reason, we will not present the mass-TLN relation in this~work.

The rest of the paper is organized as follows. In~the next section, we will sketch how we extract the EoSs of the bosonic dark matter models in the isotropic limit. In~Section~\ref{sec 3}, we discuss the stability of the dark and hybrid stars based on the famous Bardeen--Thorne--Meltzer (BTM) criteria, especially with the emphasis on the saddle instability when the mixed phase rule does not apply. The~key result is exposed in Section~\ref{sec 4} by plotting the mass--radius relation for various EoSs extracted from the respective dark matter models. Based on the mass--radius relation, we discuss the relevance to and interpretation of GW190814. Finally, we conclude our paper in Section~\ref{sec 5}.
\section{EoS for Bosonic SIDM in the Isotropic~Limit} \label{sec 2}

Most of the dark matter models are motivated by particle physics, which have either weak or no interaction with the standard model particles. The~former is called the WIMP, namely weak interacting massive particles, and~the latter will also include the self-interactions to explain the core profile of dark halos well so that it is usually called SIDM, self-interacting dark matter.  In~this paper, we will mostly focus on SIDM. The~simplest model of SIDM is the massive $\phi^4$ bosonic field theory considered in~\cite{Colpi:1986ye}. Naively, one should solve the combined scalar-tensor field equations for the possible compact dark star configurations. However, if~we assume the scalar profile inside the star varies very slowly, we can neglect the spatial profile and obtain an isotropic dark star configurations. In~the low energy regime, these kinds of isotropic configurations are more favored energetically than the nonisotropic ones. Thus, for~simplicity, we will only consider such kind of dark and hybrid stars. One additional advantage for solving the isotropic dark stars is we can first extract the EoS in the isotropic limit, then solve the Tolman--Oppemheimer--Volkoff (TOV) equations for the mass--radius relation. This is numerically far easier than solving the scalar-tensor field equations by the shooting method\footnote{A way for solving boundary value problems by changing them into initial value problems. One 'shoot' out trajectories in different directions until a trajectory with the desired boundary value is found, which usually costs a long machine time.}.

In the following, we sketch how to extract the EoS for the generic bosonic SIDM models by generalizing the procedure given in~\cite{Colpi:1986ye}. The~Lagrangian of the bosonic SIDM model considered in this work is
\begin{equation}\label{effective}
{\cal L}=-{1\over 2} g^{\mu\nu} \partial_{\mu}\phi^* \partial_{\nu} \phi -V(|\phi|)
\end{equation}
from which we can obtain the field equation
\begin{eqnarray}\label{scalar-eq}
0=\partial_\mu (\sqrt{-g}\partial^\mu \phi )-\sqrt{-g} V'(\phi),
\end{eqnarray}
where $V':=\partial_{\phi}V$.

We will consider the following metric ansatz with spherical symmetry for the dark or hybrid star
\begin{eqnarray}\label{metrica}
ds^2=-B(r) d t^2+A(r) d r^2 + r^2 d \Omega.
\end{eqnarray}

This metric is sourced by the spherically symmetric scalar field configuration
\begin{eqnarray}\label{fullscalar}
\phi(r,t)=\Phi(r)e^{-i \omega t},
\end{eqnarray}
which should solve \eq{scalar-eq} in the space-time \eq{metrica}, i.e.,~
\begin{eqnarray}
0&=&\partial_r (r^2 \sqrt{AB}\frac{1}{A}\partial_r\Phi )+\omega^2
(r^2 \sqrt{AB}\frac{1}{B}\Phi )-r^2 \sqrt{AB} V'(\Phi)\nonumber\\
&=&\partial_r (r^2 \sqrt{\frac{B}{A}}\partial_r\Phi )+r^2 \sqrt{AB}\left[
\frac{\omega^2}{B}\Phi - V'(\Phi)\right].\label{scalar}
\end{eqnarray}

On the other hand, the~stress tensor associated with Lagrangian \eq{effective} is
\begin{eqnarray}
 T_{\nu}^{\mu}={1\over 2} g^{\mu\sigma} \left(\nabla_{\sigma}\phi^* \nabla_{\nu} \phi
 +\nabla_{\sigma}\phi \nabla_{\nu} \phi^* \right)- \delta_{\nu}^{\mu}({1\over 2} g^{\rho\lambda} \nabla_{\rho}\phi^* \nabla_{\lambda} \phi+V(|\phi|))
\end{eqnarray}
which satisfies the conservation law
\begin{eqnarray}
 \nabla_{\mu} T_{\nu}^{\mu}=0,
\end{eqnarray}
and also sources the Einstein equation
\begin{equation}\label{Einstein}
G_{\mu\nu}=8\pi G_N T_{\mu\nu}
\end{equation}
with $G_{\mu\nu}$ the Einstein tensor and $G_N$ the Newton constant. The~total configurations for a boson star specified by $A(r)$, $B(r)$, and~$\Phi(r)$ should be determined by solving \eq{scalar-eq} and \eq{Einstein} together. In~general, the~stress tensor for the stationary configurations of \eq{fullscalar} satisfying \eq{scalar} in the space-time \eq{metrica} takes the following form in the co-moving frame
\begin{eqnarray}\label{Tmunu-1}
T_{\mu}^{\nu}=\textrm{diag}(-\rho,p,p_{\perp},p_{\perp}).
\end{eqnarray}

However, in~the isotropic limit, it will further reduce to the form of a perfect fluid, i.e.,~$p_{\perp}=p$.

Now, we consider the following concrete example of bosonic dark matter model with the self-interactions as
\begin{equation}\label{Vn}
V(\phi)=\frac12 m^2|\phi|^2+{1\over n} {\lambda_n \over \Phi_0^{n-4}} |\phi|^n\;.
\end{equation}

This can be thought as a model possessing of a UV $Z_n$ symmetry, which is however mildly broken at low energy by the mass term, that is, the~$Z_n$ symmetry is approximately good at low energy if $\Phi_0 \gg m$. For~this bosonic SIDM model, its stress tensor is given by
\begin{eqnarray}
  T_{\nu}^{\mu}= \left[
 \begin{matrix}
  - {1\over 2}({\omega^2 \over B}+m^2)\Phi^2 - {1\over n} {\lambda_n \over \Phi_0^{n-4}} \Phi^n - \frac{(\partial_{r}\Phi)^2}{2A }  \qquad 0 \qquad 0 \qquad~0\\
   0 \qquad  {1\over 2}({\omega^2 \over B}-m^2)\Phi^2 - {1\over n} {\lambda_n \over \Phi_0^{n-4}} \Phi^n + \frac{(\partial_{r}\Phi)^2}{2A } \qquad 0 \qquad~0\\
   0 \qquad 0\qquad  {1\over 2}({\omega^2 \over B}-m^2)\Phi^2 - {1\over n} {\lambda_n \over \Phi_0^{n-4}} \Phi^n - \frac{(\partial_{r}\Phi)^2}{2A } \qquad~0\\
   0 \qquad 0\qquad 0 \qquad {1\over 2}({\omega^2 \over B}-m^2)\Phi^2 - {1\over n} {\lambda_n \over \Phi_0^{n-4}} \Phi^n - \frac{(\partial_{r}\Phi)^2}{2A }\;
  \end{matrix}
  \right] \label{Tmunu-2}
\end{eqnarray}
which takes the form of \eq{Tmunu-1} as~expected.

To perform the isotropic limit, we first introduce the following dimensionless variables
\begin{equation}
x=rm, \qquad \Omega=\frac{\omega}{m}, \qquad \sigma=\sqrt{4\pi }\frac{\Phi}{M_{Pl}}\;.
\end{equation}

Here, $M_{Pl}$ is the Planckian mass scale which is related to the Newton constant by $G_N=1/M^2_{PL}$. Then, the~scalar field equation \eq{scalar} is reduced into
\begin{eqnarray}\label{key}
0
=\partial_x (x^2 \sqrt{\frac{B}{A}}\partial_x\sigma) +x^2 \sqrt{AB}\left[
\frac{\Omega^2}{B}\sigma  - \frac{\sqrt{4\pi }}{ m^2M_{Pl}}V'\left(\frac{M_{Pl}\sigma}{\sqrt{4\pi }}\right)\right].
\end{eqnarray}

By further introducing a new dimensionless parameter $\Lambda_n$ defined by
\begin{equation}\label{Lambda_n}
\Lambda_n=\Big(\lambda_n\frac{\Phi_0^2}{m^2}\Big)^{2\over n-2}  \Big(\frac{M_{Pl}}{ \sqrt{4\pi }\Phi_0}\Big)^2
\end{equation}
and the new scaled quantities
\begin{equation}
\sigma_*=\sqrt{\Lambda_n}\sigma, \qquad x_*=x/\sqrt{\Lambda_n},
\end{equation}
then the energy density and pressures of of \eq{Tmunu-2} in the form of \eq{Tmunu-1} can be expressed as follows:
\begin{eqnarray}\label{rho-1}
 \rho&=& {m^2 M_{Pl}^2  \over 4\pi\Lambda_n}\Big[ {1\over 2}({\Omega^2 \over B}+1)\sigma_*^2 + {1\over n}\sigma_*^n +\frac1{ \Lambda_n} \frac{(\partial_{x_*}\sigma_*)^2}{2A } \Big], \\
 p&=& {m^2 M_{Pl}^2 \over 4\pi\Lambda_n}\Big[ {1\over 2}({\Omega^2 \over B}-1)\sigma_*^2 - {1\over n}\sigma_*^n +\frac1{ \Lambda_n} \frac{(\partial_{x_*}\sigma_*)^2}{2A }\Big],
 \\
 p_{\perp}&=& {m^2 M_{Pl}^2 \over 4\pi\Lambda_n}\Big[ {1\over 2}({\Omega^2 \over B}-1)\sigma_*^2 - {1\over n}\sigma_*^n -\frac1{ \Lambda_n} \frac{(\partial_{x_*}\sigma_*)^2}{2A }\Big], \label{pperp-1}
\end{eqnarray}
and \eqref{key} is also further reduced to
\begin{eqnarray}\label{reduced key}
0
=\frac1{\sqrt{\Lambda_n}}\partial_{x_*} (x_*^2 \sqrt{\frac{B}{A}}\partial_{x_*}\sigma_*)+\sqrt{\Lambda_n} x_*^2 \sqrt{AB}\left[
(\frac{\Omega^2}{B}-1)\sigma_* - \sigma_*^{n-1}\right]\;.
\end{eqnarray}

From \eq{rho-1} to \eq{pperp-1}, it is then easy to see that the isotropic limit can be achieved if $\Lambda_n$ is large enough~\footnote{Since the UV $Z_n$ symmetry is approximately good  if $\Phi_0\gg m$, then the isotropic limit $\Lambda_n=\Big(\lambda_n\frac{\Phi_0^2}{m^2}\Big)^{2\over n-2}  \Big(\frac{M_{Pl}}{ \sqrt{4\pi }\Phi_0}\Big)^2\gg 1$ is easier to achieve. It is noticed that a large derivative of $ \sigma_*$ may balance the large $\Lambda_n$ term, but~here we only consider the simplest case.} so that the spatial kinetic term can be neglected in comparison to the other terms~\footnote{Naively, the~comparison should also take into account the spatial profile of $A(x)$, $B(x)$ and $\sigma_*(x)$. Here, we assume that $\Lambda_n$ is large enough as argued in the previous footnote so that the variation of the spatial profiles of $A$, $B$ and $\sigma_*$ will not affect the result. This of course can be justified after solving the TOV configurations.}. Thus, the~resultant approximate energy density and isotropic pressure then become
\begin{eqnarray}
 \rho&=& {m^2 M_{Pl}^2  \over 4\pi\Lambda_n}\Big[ {1\over 2}({\Omega^2 \over B}+1)\sigma_*^2 +  {1\over n}\sigma_*^n \Big], \\
 p&=& p_{\perp}={m^2 M_{Pl}^2 \over 4\pi\Lambda_n}\Big[ {1\over 2}({\Omega^2 \over B}-1)\sigma_*^2 -  {1\over n}\sigma_*^n \Big].
\end{eqnarray}

Moreover, in~this isotropic limit, both $\rho$ and $p$ depend only on $\sigma_*$ and not on its spatial derivative, and~the first term of \eq{reduced key} is suppressed with respect to the second term\footnote{The same assumption is adopted as in the previous footnote to justify the suppression.} so that it can be solved for $\sigma_*$ to yield
\begin{equation}\label{phin-s}
\sigma_*^{n-2}=\frac{\Omega^2}{B}-1 \;.
\end{equation}

Using \eq{phin-s}, the~energy density and pressure can be expressed as a pure function of $\sigma_*$,
\begin{eqnarray}
 \rho&=& \rho_{n,0} \Big[ ({1\over 2}+{1\over n})\sigma_*^n + \sigma_*^2 \Big], \label{rho-eos} \\
 p&=& \rho_{n,0} ({1\over 2}-{1\over n})\sigma_*^n ,\label{p-eos}
\end{eqnarray}
where the overall energy density scale is given by
\begin{equation}
\rho_{n,0}= {m^2 M_{Pl}^2  \over 4\pi\Lambda_n}= {4\pi m^2 \Phi_0^2 }({m \over{\sqrt \lambda_n \Phi_0}  })^{\frac4{n-2}}\;.
\end{equation}

Note that \eq{rho-eos} and \eq{p-eos} are already forming a parametric EoS for the dark matter model in the isotropic limit. However, we can further eliminate the parametric function $\sigma_*(r)$ to yield a compact form of EoS as following (in the units of $c=1$ and $\hbar=1$):
\begin{equation}\label{phin}
{\rho\over \rho_{n,0}}={n+2\over n-2}
\Big({p\over \rho_{n,0}}\Big) + \Big({2n \over n-2}  \; {p\over \rho_{n,0}}\Big)^{2\over n}.
\end{equation}

One can then adopt this EoS to solve the TOV configurations for the dark or hybrid stars. Note~that, for~$n=4$, the~above EoS reduces to the known result given in~\cite{Colpi:1986ye}, namely,
\begin{equation}\label{phi4colpi}
{\rho\over \rho_{4,0}}=3\Big({p\over \rho_{4,0}}\Big) +2 \sqrt{p\over \rho_{4,0}}\,.
\end{equation}

For the other $n$'s, the~resultant EoSs are new and not considered before in the literature. It is interesting to see that, as~$n$ goes higher, the~EoS becomes stiffer, e.g.,~as $n\longrightarrow \infty$, the~EoS goes to $p=\rho$, i.e.,~the causality~limit.

Even though we have arrived at the EoS \eq{phi4colpi}, it is written in the unit of $\rho_0$. For~the convenience of later use when solving TOV configurations, we will choose the following astrophysical units associated with the solar mass $M_{\odot}$:
\begin{equation}\label{astrophy}
r_{\odot}=G_N M_{\odot}/c^2,\quad \rho_{\odot}=M_{\odot}/r_{\odot}^3, \quad p_{\odot}= c^2 \rho_{\odot}.
\end{equation}

Then, in~terms of these units, the~EoS \eq{phi4colpi} will turn into the following form~\cite{Zhang:2020pfh}:
\begin{equation}\label{DEoS}
{\rho\over \rho_{\odot}}=3 \; \Big({p\over p_{\odot}}\Big)+{\cal B}_4\; \sqrt{p\over p_{\odot}}\,
\end{equation}
with ${\cal B}_4:= {0.08 \over \sqrt{\lambda_4}}({m\over \textrm{GeV}})^2$ a free  parameter.  Similarly, the~EoS \eq{phin} can be turn into
\begin{equation}\label{EoS_phin}
{\rho\over \rho_{\odot}}={n+2\over n-2}\; \Big({p\over p_{\odot}}\Big)+{\cal B}_n\; \Big({p\over p_{\odot}}\Big)^{2\over n}
\end{equation}
where ${\cal B}_n:=({2n \over n-2}  )^{2\over n} ({\rho_{n,0} \over p_{\odot}})^{1-{2\over n}}$ is again a free parameter related to $m$ and $\lambda_n$.

 \begin{figure}[htbp]
\centering
{\includegraphics[width=7.5cm]{psigma.pdf} }%
\quad
{\includegraphics[height=0.30\textwidth]{prho.pdf} }%
\caption{An example of EoS given by \eq{EoS46-1} and \eq{EoS46-2} with ${\cal B}=0.0006$ and ${\beta}=0.1$. Left: it shows the region of $p<0$ which should be excluded when solving TOV configurations of dark stars although it may yield gravastar configurations. Right: the plot of EoS for the phsyical region. The~unit of $p$ is $p_{\odot}$, and~the unit of $\rho$ is $\rho_{\odot}$, and~$\sigma$ is~dimensionless.  }
\label{prhosigma}
\end{figure}

The above method of extracting the EoS in the isotropic limit from a given  scalar field theory can be applied to the model with more general potentials other than \eq{Vn}.  For~example, for~the model with the following potential,
\begin{equation}
V(\phi)=\frac12 m^2|\phi|^2-\frac14 \lambda_4 |\phi|^4+\frac16 {\lambda_6 \over \Phi_0^2} |\phi|^6\,.
\end{equation}

In this case, we have two scaling parameters $\Lambda_4$ and $\Lambda_6$ as defined in \eq{Lambda_n}. However, it is easier to parameterize the EoS in the isotropic limit by the following two parameters: $\Lambda:=\sqrt{\Lambda_6}$ and $\beta:={\Lambda_4 \over 4 \Lambda}$ so that the extracted EoS takes the following form in terms of the astrophysical units of \eq{astrophy}:
\begin{eqnarray}\label{EoS46-1}
{\rho\over \rho_{\odot}}&=&{\cal B} \left(\frac23 \sigma^6-3  \beta \sigma^4+\sigma ^2\right) \,
\\ \label{EoS46-2}
{p\over p_{\odot}}&=&{\cal B} \left(\frac13 \sigma^6-  \beta \sigma^4\right)
\end{eqnarray}
where ${\cal B}:=\frac{m^2 M^2_{pl}}{4 \pi \Lambda \rho_{\odot}}$.

Some comments about the above EoS are in order: (i) one can see ${\cal B}>0$ but $\beta$ can be either positive or negative. When $\beta=0$, it reduces to the case of pure $\phi^6$ model with ${\cal B}_6=(3{\cal B}^2)^{1/3}$; (ii)~We~should require both $\rho \ge 0$ for positive energy condition and $p\ge 0$ for not considering the gravastars, thus the corresponding physical range of $\sigma$ should be chosen carefully. In~particular, when $\beta>0$, one should require $\sigma \ge \sqrt{3\beta}$ to keep $p\ge 0$. An~example for this case is shown in Figure~\ref{prhosigma}. On~the other hand, both~$\rho$ and $p$ are positive for all ranges of $\sigma$ when $\beta<0$.

Below, we will study the TOV configurations of dark and hybrid stars made of nuclear and dark matters based on the above EoSs of bosonic SIDM model in the isotropic limit. In~particular, we will focus more on the hybrid stars with saddle instability which causes more marginal parameter space to constrain the dark matter models by the observed LIGO/Virgo gravitational wave events. In~this way, we demonstrate how to constrain the particle physics models of dark matter by the gravitational wave astronomical~observations.

\section{BTM Criteria and Saddle Instability for Hybrid~Stars} \label{sec 3}
Neutron stars and some other compact stars are relativistic objects that their structure should be analyzed using general relativity.
TOV equations~\cite{TOV-1,TOV-2,Ott13} can be applied for these kinds of calculations, which are
derived from the Einstein equations and the conservation of the energy--stress tensor,
assuming zero space velocity, spherical symmetry, and~an ideal fluid~model.

In this paper, we will mainly consider the hybrid stars by solving the following TOV configurations with multiple component fluids inside the star~\cite{ Mukhopadhyay:2016dsg, Rezaei:2018cuk}:
\begin{equation}
{dp_I \over dr}=-(\rho_I + p_I) {d \phi\over dr},\quad  {d m_I\over dr} =4\pi r^2 \rho_I, \quad {d \phi \over dr} ={m +4\pi r^3 p \over r(r-2m)},
\end{equation}
where the index $I$ labels the fluid components, and~the total pressure and energy density are given by $p:=\sum_I p_I$  and $\rho:=\sum_I \rho_I$, respectively; and $m(r)=\sum_I m_I(r)$ is the mass profile inside the star, and~the Newton potential $\phi:={1\over 2} \ln(-g_{tt})$ with $g_{tt}$ the $tt$-component of the metric. In~this paper, we~will consider the cases with only two-component hybrid stars made of nuclear matter and dark~matter.

The resulting configurations of hybrid stars depend on if there are strong interactions between component fluids. If~there are strong interactions, it is expected to form a domain wall between phases of different fluids. Following our previous consideration in~\cite{Zhang:2020pfh}, we refer the hybrid stars with neutron core and dark matter crust as scenario I, and~the ones with dark matter core and neutron crust as scenario II. On~the other hand, if~there is almost no interaction among component fluids, then all components prevail inside the star and mix together right from the core. We refer to this case as scenario III~\cite{Zhang:2020pfh}. To~solve the TOV equations, we need to provide the core pressure as the initial condition, and~then, by~changing the initial pressure, we will obtain different TOV configurations to plot the mass--radius relation. For~scenarios I and II, one only needs to provide one initial condition for the core pressure since only one kind of fluid dominating the core, but~needs two for scenario~III.

Not every TOV configuration is stable, and~one needs to judge its stability by either solving the Sturm--Liouville eigenmodes of radial oscillation (another method of judging the equilibrium configuration is used in~\cite{Henriques:1990xg, Valdez-Alvarado:2020vqa, DiGiovanni:2020frc}) or by some equivalent empirical criteria. One set of such criteria for single component fluid is the so-called BTM (Bardeen--Thorne--Meltzer) criteria~\cite{BTM}, which~states that, in~the direction of increasing core pressure along the mass--radius curve, whenever an extremum is passed, one stable mode becomes unstable if the curve bends counterclockwise. Contrarily, one~unstable mode becomes stable if it rotates clockwise. These criteria were originally formulated by starting from the stable planet configuration with low enough core pressure. It, however, will~cause some trouble if one does not start from such kind of planet configuration when solving the TOV equations. Thus, in~practicality, it is useful to formulate the reverse BTM criterion by traveling the mass--radius curve in the direction of decreasing core pressure, and~to require the consistency with the BTM ones. This then leads to following { (Reverse)
 BTM stability criteria:} { Traveling the mass--radius curve along either directions of increasing or decreasing core pressure, whenever passing an extremum, one~stable mode becomes unstable if the cure bends counterclockwise; otherwise, one~unstable mode becomes stable.}

By applying the (reverse) BTM criteria, one can make sure that, for~the unstable regions on the mass--radius curve, not even one starts from a known stable region. However, one can only pin down the stable regions if one starts from a stable region, or~knows how many unstable modes there are. For~example, for~the curve $O'A'D'B'C'$ on the left panel of Figure~\ref{BTM}, one can apply the (reverse) BTM criteria criteria to make sure $O'A'D'$ is unstable regardless if $D'B'C'$ is stable or not, even if seems impossible to apply the original BTM criterion as there is no asymptotic stable region. On~the other hand, for~the curve $OABC$ in Figure~\ref{BTM}, one cannot be sure if $ABC$ is stable even if one can make sure that $OA$ is unstable by applying the (reverse) BTM~criteria.

\begin{figure}[htbp]
\centering
{\includegraphics[height=0.44\textwidth]{BTM.pdf}
\quad
\includegraphics[height=0.26\textwidth]{saddle.pdf}
}
\caption{The various curves and intersections of mass--radius relations used to demonstrate the (reverse) BTM criteria in the main text. Solid lines indicate stable parts, while the dotted lines are for unstable parts. Left panel: Various mass--radius curves. The~thick arrows indicate the directions of increasing core pressures for each branch. Right panel: Summary of the situations of intersections. $X_1$~is a stable point where both the neutron branch (black) and the dark branch (brown) are stable.  $X_2$ is unstable since both branches are unstable.  $X_3$ and $X_4$ show the saddle instability that one of the two branches is unstable. In~some cases, the~two branches may cross twice like points $Y$ and $Z$; then, at~ least one of them is a fake~intersection. }
\label{BTM}
\end{figure}

Importantly, the~above criteria is only checked rigorously for the stars with single component fluid. One should be careful when applying the (reverse) BTM criteria to judge the stability of hybrid stars. For~the hybrid stars of scenarios I and II, there is only one dominant fluid component in each region, one will then expect that the (reverse) BTM criteria should still work without the need of much modification; see~\cite{Kain:2020zjs} for the recent~discussion.

On the other hand, for~the hybrid stars of scenario III, one needs to set the core pressures for both fluid components, say one for nuclear matter and one for dark matter.  Thus, one needs to plot two kinds of mass--radius curves, namely, one by fixing the neutron core pressure but tuning the dark matter core pressure, and~the other by fixing the dark matter core pressure but tuning the neutron core pressure. Let us call the former the dark branch and the latter the neutron branch. In~the left panel of Figure~\ref{BTM}, we show some representatives of both branches, e.g.,~the curve $o'a'b'c'$ is the neutron branch, and~the three curves in the same figure are the dark branches\footnote{Neutron branches could also look like the two midlle curves $O'A'D'B'C'$ and $OABC$, but~with oppisite pressure directions, like the solid curves in {Figures} \ref{phi4}--\ref{phisix}.
}. For~any intersection point of dark and neutron branches, e.g.,~point $X_1$ to $X_4$ on the right panel of Figure~\ref{BTM}, there could be saddle instability\footnote{For a single-component star, the~stability is solely determined by altering the core pressure. When the core pressure has some perturbation, if~the radial oscillation is stable, then the star is stable. However, for~a two-component star, the~core pressure is a sum of two partial pressures of each component. Then, a~necessary condition for such a star being stable is that it is stable when changing either of the partial pressure while keeping the other fixed. It might happen that a star is stable when changing the neutron (or dark matter) pressure but unstable when changing the other partial pressure. This will still lead to an unstable star and we call it the saddle instability. } if one of the branches is unstable at this intersection point.  The~difficulty is how to judge the stability for each branch. Without~rigorous study of the Sturm--Liouville eigenmodes of radial oscillation, it is hard to answer this question. In~this work, we assume that the (reverse) BTM criteria still work for each dark or neutron branch, and~apply the criteria to find out the saddle unstable regions. Furthermore, we shall also assume that the cusp points such as $B$ and $B'$ ($B'$ is not even an extremum) on the left panel of Figure~\ref{BTM} will not induce the change of stability of any radial oscillation eigenmode. Otherwise, there could be more complications~\cite{Alford:2017vca}.

Even with the above assumptions holding when traveling along the mass--radius curve, to~judge the stability with (reverse) BTM criteria, it is better to start from some known stable region. One such region is the stable part of the pure neutron star curve, and~the other is the stable part of the pure dark star curve. We then need to further assume that the small doping with the other component will not change the stability; it then implies that the nearby regions of the stable part of single component fluid stars are also stable. For~example, if~the curve $o'a'b'c'$ is a nearby curve around the pure neutron star, and~its stable part (the solid line $a'b'$) is inherited from the stable part of the pure neutron star, then one can infer the unstable part (the dashed line $o'a'$ and $b'c'$) by applying the (reverse) BTM criteria. Furthermore, the~nearby regions such as $C'B'D'$ and $CBX$ should also be stable. The~reason is as follows: these two regions branch out from an almost pure neutron star curve, and~are at the ends of dark branches, this means that the core pressures of the dark matter component are negligible. Thus, these two parts are indeed almost pure stable neutron stars. We can then infer that: (1) $D'A'O'$ is unstable because $B'D'A'O'$ bends counterclockwise when passing around local extrema $D'$ and $A'$; and (2)$XA$ is stable as there is no local extremum on it. (3) $A'O'$ is unstable because it bends counterclockwise when passing through the local maximum $A$.

Based on the above discussions, we can then apply the reverse BTM criteria to the intersection points of the neutron and dark branches, such as the point A on the right panel of Figure~\ref{BTM} to judge if the intersection point is saddle stable or not. Both curves crossing at $A$ being stable is a necessary condition for $A$ to be stable. In~addition, there are some fake intersection points such as point $X$ for the case when one curve intersects with the other curve more than once. Then, only the end point/branch point is real intersection and others are fake. We have summarized these situations in the right panel of Figure~\ref{BTM}. We will apply the lessons from the above discussions to the case studies in the next section to judge the stability of hybrid stars, especially the saddle instability of scenario~III.

\section{Dark Star and and Hybrid Star~Interpretations} \label{sec 4}

In this section, we consider the dark stars and hybrid stars of all three scenarios based on the EoSs extracted from the bosonic SIDM models discussed in Section~\ref{sec 2} in the isotropic limit. For~the EoS of the nuclear matter, the~standard  phenomenological neutron EoS SLy4~\cite{Douchin:2001sv,SLy4} is employed, which has a maximum mass of 2.05 $M_{\odot}$, and~the radius at 1.4 $M_{\odot}$ is about 12km. This EoS is offered as discrete sets of pressure and energy density, which is convenient to be dealt with numerically.
 The main purpose is to see if some of the dark or hybrid stars can reach the mass of 2.6 $M_{\odot}$ or so to explain the smaller companion compact object of GW190814. It is easy to see that this purpose can be easily achieved for the dark stars and hybrid stars of scenarios I and II by choosing appropriate parameters of SIDM. On~the other hand, it is more difficult for the hybrid stars of scenario III to reach such a mass because of the saddle instability. It is interesting to see that the astronomical observations of gravitational waves can rule out some theoretical scenarios of hybrid stars and the associated dark matter models. Thus, we will focus more on scenario III in the later discussions of this~section.

Now, we present the mass--radius relations case by~case.

 \begin{figure}[htbp]
    \centering
{\includegraphics[height=0.55\textwidth]{rdrn.pdf}}
\quad
{\includegraphics[height=0.55\textwidth]{rdrnL.pdf}}%
\caption{{Mass-Radius} (up) and TLN-Mass (down) relations for the hybrid stars of scenarios I and II. The~red line stands for the pure neutron stars with SLy4 EoS,  while the green and blue lines stand for the pure dark stars of EoS given by \eq{DEoS} with ${\cal B}_4 =0.035$ and $0.045$, respectively. The~other lines represent hybrid stars of neutron and dark matters associated with the above EoSs. For~${\cal B}_4$ = 0.035, they are labelled by {\bf aRN} = $r_W$ (thin brown lines) for hybrid stars of neutron core, and~{\bf aRD} = $r_W$ (thin~ black lines), where $r_w$ is the core radius. Similarly, for~${\cal B}_4 =0.045$, they are by {\bf bRN} = $r_W$ (thick brown lines) and {\bf bRD} = $r_W$ (thick black lines). For~example, {\bf bRN} = 6 labels the hybrid stars of $r_W=6 $ km. The~unstable configurations are denoted by the dashed lines. The~bottom panel is presented to show the orders of the corresponding TLN, and~we find that 2.6 $M_{\odot}$ stars are well below 4000, which are negligible after averaging with a 23 $M_{\odot}$ black hole for GW190814. The~unit of mass is $M_{\odot}$, and~the unit of radius is km, and~TLN $\Lambda$ is~dimensionless.}
\label{rdrn}
\end{figure}

{Bosonic
 $\phi^4$ model. }
We first consider dark and hybrid stars based on the $\phi^4$ EoS, i.e.,~\eq{DEoS}, which is extracted from the $\phi^4$ dark SIDM in the isotropic limit. For~scenarios I and II, i.e.,~assuming interactions exist between neutron and dark matter, it is easy to form a 2.6 $M_{\odot}$ star, as~shown in Figure~\ref{rdrn}, in~which we also show the TLN-mass relations for one's reference. Pure neutron stars of SLy4 EoS are marked in red, and~pure dark stars with ${\cal B}_4=0.035$ and $0.045$ are marked in blue and green, respectively. The~other lines represent hybrid stars of neutron and dark matter with these EoS.  When~${\cal B}_4 =0.035$, they are labelled  by {\bf aRN}  = $r_W$ (brown) for neutron core case and by {\bf aRD} = $r_W$ (black) for dark matter core case. Here, $r_W$ stands for the core radius. Similarly, when ${\cal B}_4 = 0.045$, they~are labelled by {\bf bRN} = $r_W$ (brown) and {\bf bRD} =$r_W$ (black)). For~example, {\bf bRN} = 6 means the radius of neutron core is  $6\textrm{km}$, and~${\cal B}_4=0.045$. The~unstable configurations are denoted by dashed~lines.

We find that this dark star model can cover any mass range by adjusting the free parameter ${\cal B}_4$. In~particular, for~${\cal B}_4 \leq 0.047$, the~maximal mass exceeds $2.6$ $M_{\odot}$. The~maximal mass grows without an upper limit as the value of ${\cal B}_4$ drops, and, for~example, when ${\cal B}_4 =0.035$, we have $3.5$ $M_{\odot}$. Since~for scenarios I and II the 2.6 $M_{\odot}$ hybrid stars can always be achieved as long as the pure dark star has a maximal mass higher than 2.6 $M_{\odot}$, in~the following, we only concentrate on the more nontrivial scenario III, which is more realistic since it is assumed that there is no interaction between SIDM and baryonic~matter.

For scenario III, by~applying the (reverse) BTM stability criteria to the resultant mass--radius relation shown in Figure~\ref{phi4}, no stable stars around 2.6 $M_{\odot}$ can be formed.
The pure SLy4 neutron stars are marked in red as the reference configurations. For~the hybrid stars, ${\cal B}_4=0.045$ case is denoted by the brown lines,  and~$0.035$ case by the green lines. Now, we need the core pressures for both dark and nuclear matters to solve the TOV configurations. By~tuning one of the core pressures and fixing the other, we can obtain the so-called dark branch and neutron branch as discussed in Section~\ref{sec 3}. In~ Figure~\ref{phi4}, the~dark branches are denoted by the dash-dotted lines and the neutron branches by the solid~lines.

To apply the method in Section~\ref{sec 3} to check the saddle (in)stability, we can apply the (reverse) BTM criteria to both neutron and dark branches and then determine the saddle stability for their intersection points. As~discussed in Section~\ref{sec 3}, we also need to assume that the nearby regions of the pure stable star configurations are stable when applying (reverse) BTM criteria.  After~checking this way for most of the intersections in Figure~\ref{phi4}, we can determine the stable regions and unstable regions. For~example, the~point $A$ is unstable because the line $BDA$ turns counterclockwise at point $D$ towards $A$, which~makes the $DA$ part unstable. It turns out that the stable regions are small and confined near the pure neutron star configurations\footnote{It is interesting to notice that the region near pure dark star (like point A) is unstable. This should be understood that, although~the dark matter components are stable, the~neutron parts are unstable, which makes the total configuration unstable.  }, which for ${\cal B}_4=0.045$ is roughly indicated by the regions inside the highlighted blue closed line. { From
 the above analysis, we can conclude that, for~the bosonic $\phi^4$ model, however the parameter ${\cal B}_4$ varies, the~stable region cannot yield the mass of the star more than 2.1 $M_{\odot}$, which is just the maximal mass of pure neutron star associated with SLy4 EoS.} This could be due to the fact that the EoS associated with $\phi^4$ model is not stiff enough when compared with SLY4. This can also be seen from the fact that the radius of a dark star of 2.6 $M_{\odot}$ is at least 25 $\textrm{km}$, far larger than the standard neutron star's radius around 11 $\textrm{km}$.

\begin{figure}[htbp]
\centering
\includegraphics[height=0.5\textwidth]{phi4b.pdf}
\caption{Mass--radius relation for the hybrid stars of scenario III, with~SLy4 EoS for neutron and $\phi^4$ model EoS given by \eq{DEoS} for dark matter. The~pure SLy4 neutron stars are marked in red for one's reference. The~hybrid star configurations are shown in brown for ${\cal B}_4=0.045$, and~in green for ${\cal B}_4=0.035$. Solid lines represent the neutron branches, and~dash-dotted lines the dark branches. The~blue-circled regions roughly indicate the stable hybrid star configurations for ${\cal B}_4=0.045$. Saddle~instabilities can be checked at each intersection point of solid and dash-dotted lines with the same color. For~example, point A is unstable since the solid brown line on it is unstable, though~the dash-dotted line on it is stable. The~unit of mass is $M_{\odot}$, and~the unit of radius is~km.}
\label{phi4}
\end{figure}

\begin{figure}[htbp]
    \centering
\includegraphics[height=0.5\textwidth]{phi6b.pdf}
\caption{Mass--radius relation for the hybrid stars of scenario III, with~SLy4 EoS for neutron and $\phi^6$ model EoS given by \eq{EoS_phin} with $n=6$ for dark matter. The~line styles are the same as in Figure~\ref{phi4} but with the brown lines for ${\cal B}_6=0.011$ and the green lines for ${\cal B}_6=0.008$. The~blue-circled regions roughly indicate the stable hybrid star configurations for ${\cal B}_6=0.011$. Saddle instabilities can be checked at each intersection point of solid and dash-dotted lines with the same color. The~unit of mass is $M_{\odot}$, and~the unit of radius is~km.}
\label{phi6}
\end{figure}

{ Bosonic $\phi^n$ model. }
From the expression of $\phi^n$ EoS, i.e.,~\eq{EoS_phin}, we know it becomes stiffer as $n$ grows. This makes it possible to have more massive stable hybrid star configurations. Indeed, from~Figures \ref{phi6} and  \ref{phi10}, we confirm that the hybrid stars with the same maximal mass have much smaller radius than the $\phi^4$ case. For~example, considering the curves with the maximal mass 2.6 $M_{\odot}$, we find that the radius is 16 $\textrm{km}$ for $\phi^6$ model and 13 $\textrm{km}$ for $\phi^{10}$ one, which are all much shorter than the 25~$\textrm{km}$ of $\phi^4$ model.

In Figure~\ref{phi6}, it shows the mass--radius relations for the hybrid stars of scenario III made of nuclear matter of SLy4 EoS and dark matter of $\phi^6$ model's EoS with ${\cal B}_6=0.011$ (brown lines) and $0.008$ (green lines). The~pure SLy4 neutron stars are marked in red as before. Similarly, the~neutron and dark branches are denoted by solid and dash-dotted lines, respectively; and then we apply the (reverse) BTM criteria to judge the saddle (in-)stabilities as before. We find that the $BDA$ line marked in Figure~\ref{phi6} bends higher than that in Figure~\ref{phi4}, but~still point $A$ is unstable because of the existence of the minimum point $D$, as~Section~\ref{sec 3} tells us. Again, the~stable regions are near the pure star configurations, and~for the case of ${\cal B}_6=0.011$ are circled by the blue closed path.  { Again,
 we find that the maximal mass of the stable hybrid stars of this kind cannot be higher than the maximal mass, i.e.,~2.1 $M_{\odot}$ of pure neutron stars, similar to the $\phi^4$ case. }

{ However, we see the above mass bound is lifted when considering the hybrid stars of scenario III made of nuclear matter SLy4 EoS and the dark matter of $\phi^{10}$ model's EoS, and~the maximal mass can be around 2.6 $M_{\odot}$ to explain GW190814. We have also considered the hybrid star configurations of scenario III for $\phi^8$ model (not shown here) and reach a maximal mass about 2.3 $M_{\odot}$.} The mass--radius relations are shown in Figure~\ref{phi10} for $\phi^{10}$ model's EoS with ${\cal B}_{10}=0.0036$ (brown lines) and $0.0028$ (green lines). The~line styles are the same as in Figures~\ref{phi4} and  \ref{phi6}. The~main reason for the lifting of the mass bound is the disappearance of the local minimum $D$ in Figures~\ref{phi4} and  \ref{phi6}, and~now we see that the intersection point $A$ in Figure~\ref{phi10} is no longer a saddle point so that the stable region is extended beyond the nearby region of the pure neutron star curve (the red curve). The~stable region for ${\cal B}_{10}=0.0036$ is again circled by the blue closed path, however, in~which we see the maximal mass is about  2.6 $M_{\odot}$. This can then be used to explain the heavy companion of the black hole in~GW190814.

\begin{figure}[htbp]
    \centering
\includegraphics[height=0.5\textwidth]{phi10b.pdf}
\caption{Mass--radius relation for the hybrid stars of scenario III, with~SLy4 EoS for neutron and $\phi^{10}$ model EoS given by \eq{EoS_phin} with $n=10$ for dark matter. The~line styles are the same as in Figure~\ref{phi4} but with the brown lines for ${\cal B}_{10}=0.0036$ and the green lines for ${\cal B}_{10}=0.0028$. The~blue-circled regions roughly indicate the stable hybrid star configurations for ${\cal B}_{10}=0.0036$. Saddle instabilities can be checked at each intersection point of solid and dash-dotted lines with the same color. The~unit of mass is $M_{\odot}$, and~the unit of radius is~km.}
\label{phi10}
\end{figure}

{Bosonic $\phi^4+\phi^6$ model. }
The EoS for this model is given in \eq{EoS46-1} and \eq{EoS46-2}. Naively, it seems to be the intermediate model between $\phi^4$ and $\phi^6$ models. However, as we have seen in Section \ref{sec 2}, it contains two tuning parameters $\cal B$ and $\beta $ which lead to some novelty: there are some regions with $p<0$ or even $\rho<0$. We will only consider the region for positive $p$ and $\rho$.

In Figure \ref{phisix}, it shows the mass--radius relations for the hybrid stars of scenario III made of nuclear matter of SLy4 EoS and dark matter of this EoS with ${\cal B}=0.0006$, ${\beta}=0.1$ (brown lines), ${\cal B}=0.0004$, ${\beta}=0.1$ (green lines) and ${\cal B}=0.0006$, ${\beta}=-0.1$ (purple lines). The line styles are the same as in Figure~\ref{phi4}. We then run through the same check of saddle (in-)stabilities as before to pin down the stable and unstable regions by applying the (reverse) BTM criteria. As a result, we find that the maximal mass is mainly affected by the value of ${\cal B}$, with lower ${\cal B}$ leads to higher mass, similar to the behavior of ${\cal B}_n$. On the other hand, $\beta$ affects the compactness, namely, the radius becomes smaller as $\beta$ grows from negative values to the positive ones. However, the change of radius  significantly slows down when $\beta \geq 0.1$.  Though the $\beta=0$ case is not included in Figure \ref{phisix}, this case is equivalent to $\phi^6$ model with ${\cal B}_6=(3 {\cal B}^2)^{1/3}$ as discussed in Section \ref{sec 2}. For example, ${\cal B}=0.0006$ corresponds to ${\cal B}_6=0.010$, which is comparable to the ${\cal B}_6=0.011$ case (brown lines) in Figure \ref{phi6} .

As shown in Figure \ref{phisix}, although the EoS becomes considerably stiffer when $\beta$ increases, we still cannot have stable stars around 2.6 $M_{\odot}$ because there is a local minimum at $D$ so that $A$ is a saddle point. The blue closed path roughly encloses the stable region for ${\cal B}=0.0006$, ${\beta}=0.1$. By varying $\cal B$ and $\beta$, the maximal mass of the stable region is about 2.4 $M_{\odot}$, when ${\cal B}=0.0007$ and ${\beta}=0.1$.

\begin{figure}[htbp]
\centering
\includegraphics[height=0.5\textwidth]{phisixa.pdf}
\caption{Mass--radius relation for the hybrid stars of scenario III, with SLy4 EoS for neutron and $\phi^4+\phi^6$ model EoS given by \eq{EoS46-1} and \eq{EoS46-2} for dark matter. The line styles are the same as in Figure~\ref{phi4} but with the brown lines for ${\cal B}=0.0006$, ${\beta}=0.1$, the green lines for ${\cal B}=0.0004$, ${\beta}=0.1$, and purple lines for ${\cal B}=0.0006$, ${\beta}=-0.1$. The blue-circled regions roughly indicate the stable hybrid star configurations for ${\cal B}=0.0006$, ${\beta}=0.1$. Saddle instabilities can be checked at each intersection point of solid and dash-dotted lines with the same color. The unit of mass is $M_{\odot}$, and the unit of radius is km.}
\label{phisix}
\end{figure}

\section{Conclusions} \label{sec 5}
In this paper, we have extended the usual study of dark and hybrid stars for $\phi^4$ SIDM to more general types of bosonic SIDM models by extracting their EoSs in the isotropic limit so that we can have more access to the complete mass--radius relations for the TOV configurations. Among them, we are especially interested in the $\phi^n$ models which can have a stiff EoS when $n$ is large enough, say~$n=10$. These kinds of models can be motivated by the UV $Z_n$ flavor symmetry, which may be a natural symmetry for some higher theories. It is fascinating to further explore this connection between particle physics and the resultant astrophysical compact objects in the future.

In general, it is easy to tune the parameters in SIDM to yield compact objects with masses comparable to or higher than the ones of neutron stars. Therefore, it is interesting to check the dark star possibilities by the future observations of the compact objects via the gravitational wave observations. Similar conclusions can be reached for the hybrid stars of scenarios I and II made of nuclear and bosonic SIDMs. In particular, many such kinds of dark and hybrid stars can have masses more than 2.6 $M_{\odot}$, and thus can be adopted to explain the recent GW190814 in which one companion compact object with such a mass has been identified by parameter estimation, which is hard to be explained by the usual EoSs for neutron stars.

On the other hand, the hybrid stars of scenario III are subjected to the saddle instability, and it is difficult for such kind of stars to have higher mass such as 2.6 $M_{\odot}$. However, in this paper, we do find such a stable massive configuration when the $n$ of the $\phi^n$ model rises up to $10$ for which the saddle instability around 2.6 $M_{\odot}$ is lifted. Although in practice it is hard to tell if a compact object can be the hybrid stars of scenario III in the near future, we may hope that the obstacle will be overcome in the long run to have precise measurements on the structure of the stars via gravitational wave detection to pin down the scenarios. Despite that, theoretically, it is interesting to see the existence of some mass bound by the saddle instability. In this paper, we simply assume that the BTM criteria still work for the mixing fluids without the mutual interaction, applying it to judge the saddle instability. It will be illuminating to study the saddle instability by directly examining the eigenspectrum of the radial oscillation modes.

Overall, our work demonstrates the intriguing interplay between the particle physics models of dark matter and the astrophysical observations via the gravitational wave detection.  We hope this will encourage more works to explore the dark matter physics via the study of dark and hybrid stars in the new era of gravitational astronomy.

~
\section*{Acknowledgement}
FLL is supported by Taiwan Ministry of Science and Technology (MoST) through Grant No.~109-2112-M-003-007-MY3. KZ (Hong Zhang) is supported by MoST through Grant No.~107-2811-M-003-511. We thank Guo-Zhang Huang and other TGWG members for helpful discussions. We also thank NCTS for partial financial support.


\section{Introduction}

Deep convolutional neural networks(deep CNNs)\cite{b1,b2,b3,b4} have achieved outstanding results in many computer vision tasks. However, deep CNNs comes with a huge computational cost, which limits application on embedded devices(i.e. mobile phone). \\
To expand the application scope of deep CNNs, channel pruning methods were proposed. Traditional channel pruning methods always rely on human-design rules\cite{b6,b7}. Recently, inspired by the Neural Architecture Search(NAS), some AutoML-based pruning work\cite{b8,b9,b12} has been proposed to automatically prune channels without a human-design mode. Considering a network with 10 layers(each layer contains 32 channels), the candidates of each layer and the whole network could be $32$ and $32^{10}$, respectively. Thus, AutoML-based pruning methods can be seen as fine-grained NAS because of more candidates than normal NAS\cite{b28,b27,b23} in each layer. \\
For above mentioned AutoML-based methods, some based on reinforcement learning or evolutionary algorithm\cite{b19,b8,b7} are quite time-consuming due to iterative retrain for each pruned network. To reduce the computation in pruning, some AutoML-based work\cite{b35,b10,b9,b12} shares the weights through training a supernet for all candidate pruned networks which are called subnets. A typical weight-sharing pruning approach contains three steps: training a supernet by iteratively sampling and updating different candidates, searching the best subnet based on the evolutionary or greedy algorithm and training the best subnet from the scratch. However, Chu\cite{b13} considered that the weight-sharing method causes unfull training in the first step since each candidate(subnet) has only a small sampling probability in training. Moreover, unfull training leads to an inaccurate evaluation in the second step, which means some candidates perform well on weight-sharing while bad on training from the scratch. The problem of inaccurate evaluation is particularly obvious in AutoML-based pruning because of more candidates contained in the supernet. Moreover, the width of a supernet also has an impact on the effectiveness of evaluation, which is analyzed in \textbf{Sec. 3.1}\\
To solve the above-mentioned problem, a stage-wise training and searching approach is proposed in this paper. Inspired by \cite{b14}, we consider to divide a deep supernet into several stage-wise supernet(i.e. ResNet50\cite{b3} consists of 4 stages) for reducing the depth of the supernet. Since each stage of the supernet can be trained and searched independently, the candidate number in one stage is an exponential reduction compared with the whole supernet. With a small search space in each stage, the sampling probability of the candidates is raised, which means that each supernet can be fully trained. To alleviate the unfair training result caused by the width of the supernet, both the fullnet and the tinynet, which has the largest and smallest width respectively, are trained in each iteration. Besides, we present a distributed evolutionary algorithm where each stage can be searched independently in terms of an evolutionary algorithm(EA). The constraints(i.e. FLOPs, latency) for each EA are provided by another EA, called EA manager, where the EA manager searches for the best combination of FLOPs in each stage. Due to small and independent stage-wise search space, each EA can be sped up in a parallel way.\\
However, the stage-wise supernet in each stage cannot be trained without internal ground truth. To solve such the problem, an existing pre-trained neural network was used to generate stage-wise feature maps that are viewed as ground truth in each stage. Nevertheless, It is also time-consuming to obtain a pre-trained neural network as the teacher network. Besides, the structural difference between the teacher network and the student networks has a strong impact on distillation results. Hence, a stage-wise inplace distillation method is put forward for the fullnet(the largest width subnet) to supervise the learning of subnets. It is worth noting that the fullnet is jointly trained with other subnets. Thus, there does not exist an extra cost for obtaining the fullnet.\\
Our contribution lies in four folds:
\begin{figure*}[htbp]
\centering
\subfigure[]{
\begin{minipage}{0.4\linewidth}
\includegraphics[width=1.0\linewidth]{weight_sharing1.png}
\end{minipage}}
\subfigure[]{
\begin{minipage}{0.4\linewidth}
\includegraphics[width=1.0\linewidth]{weight_sharing2.png}
\end{minipage}}
\caption{(a) The expectation of top-1 accuracy collected from ResNet56\cite{b3} with different candidate number in one layer. The blue and red dash line denotes ResNet56 trained on CIFAR10\cite{b33} under 100 epoches and 500 epoches, respectively. (b) The expectation of top-1 accuracy collected from ResNet18\cite{b3} , ResNet34\cite{b3}  and ResNet56 with different candidate number in one layer. All models are trained under 100 epoches on CIFAR10\cite{b33}.}
\label{1}
\end{figure*}
\begin{itemize}
   \item We propose a stage-wise training and searching pipeline for both channel pruning and NAS. By splitting a CNN into several stages, the number of stage-wise candidates is exponential reduction in contrast to the net-wise candidates. Hence, each candidate obtains full training, which is the essence of accurate evaluation for searching.
  \item To conveniently provide stage-wise ground truth for each stage, a stage-wise inplace distillation method is presented through the joint training of the fullnet and subnets. Thus, the fullnet can easily supervise the learning of the subnets by offering stage-wise feature maps.
  \item To accelerate the searching process, a distributed evolutionary algorithm is suggested. Each stage can be searched by the EA with constraints given by an EA manager in a parallel way.
  \item Experiments shows that the proposed method can enhance the ranking effectiveness of searching and achieve the state-of-the-art in several datasets.
\end{itemize}
\section{Related Works}
\label{sec:section_name}
\textbf{Neural Architecture Search.} The purpose of neural network structure search is to automatically find the optimal network structure with reinforcement learning(RL)\cite{b38,b21}, evolutionary algorithm(EA)\cite{b22}, gradient \cite{b23,b24,b25} and parameter sharing methods\cite{b26,b27,b28}. RL-based and EA-based methods need to evaluate each sampled network by retraining them on the dataset, which is time-consuming. The gradient-based method can simultaneously train and search the optimal subnet through assigning a learnable weight to each candidate operation. However, the gradient-based approach causes unfair training results because some candidates obtain more learning resources than others\cite{b13}. Moreover, since gradient-based approaches need more memory for training, it cannot be applied to the large-scale dataset. Parameter sharing methods can search on the large-scale dataset by only activating one candidate in each training iteration. Nevertheless, parameter sharing methods cause unfull training results\cite{b13}. Unfull or unfair training results will cause an inaccurate evaluation of searching, which means that the best-searched architecture is not the optimal one after retraining. To solve such a problem, Li\cite{b14} proposed a blockwisely searching method, which can more fully train each sampled subnets. \\
\textbf{Pruning for CNNs.} Pruning some redundant weights is a prevalent method to accelerate the inference of CNNs. According to the different granularity of pruning, it is divided into weight pruning and channel pruning. For weight pruning, the individual weights in the channel are removed based on some rules\cite{b29}, which causes unstructured sparse filters and cannot be accelerated directly on most hardware. Therefore, much recent work focus on channel pruning. Channel pruning methods \cite{b30,b32,b19,b18,b6} can accelerate the inference of CNNs on general-purpose hardware by reducing the number of filters since the remaining filters are structural. Though the above methods achieve remarkable improvement in the practicality of pruning, it still needs human-designed heuristics to guide pruning.\\
\textbf{Auto-ML pruning.} Recently, inspired by NAS work, AutoML pruning methods\cite{b8,b9,b12,b35,b10} have attracted growing interest in automatically pruning for deep CNNs. Different from NAS, the candidate choices are consecutive in the channel pruning task. Compared with pruning methods based on the human-craft rule, AutoML pruning methods aim to search for the best configuration without manual tuning. AMC\cite{b8} adopted a DDPG agent to sample a pruned network. And the performance of the pruned networks is evaluated by training from the scratch, which is time-consuming and cannot be applied to a large-scale dataset. MetaPruning\cite{b9} trained a PruningNet that can predict weights for any pruned networks, while the parameter amount of the PruningNet is several times of the target network, which leads to unfull train. AutoSlim\cite{b12} first trained a slimmable network\cite{b35} where the weights between different widths are shared through the supernet, and then searched the best subnet in terms of the greedy algorithm. However, the width of the convolutional layer in each subnet must be the same in training. This leads to the problem that the best subnet achieves the highest accuracy with weight sharing, but at the same time, the best subnet gets poor performance when trained from the scratch. To keep the consistency of searching and retraining results, the proposed stage-wise pruning method splits a CNN into several stages and trains them separately under the supervision of the fullnet, which will be explained in Sec.3.\\
\textbf{Knowledge Distillation.} Knowledge distillation is used to train the small student model on a transfer feature set with soft labels or intermediate representations provided by the large teacher model. Soft targets lead to the superior performance of knowledge distillation\cite{b45}. However, as the network is designed to become deeper and deeper, it is not enough to just transfer knowledge to a student network from a teacher network by soft targets. To solve such a problem, some previous work\cite{b46,b47,b48,b49,b50} transferred the knowledge for the student network from the internal representation of the teacher network. All existing work assumed that the teacher network has been pretrained. Nevertheless, it is always time-consuming to train from scratch to obtain a pretrained teacher network that is essential for various knowledge distillation methods. For example, it may cost more than 10 GPU days to train a ResNet on ImageNet. Moreover, Liu\cite{b20} found that the architecture of the teacher and the student networks has a huge impact on transferring results. Hence, we proposed a stage-wise inplace distillation method to overcome the gap and to reduce the time consumption.
\begin{figure*}[t]
\centering
\includegraphics[width=1.0\linewidth]{train.pdf}
\caption{Illustration of the stage-wise training. There are three forms of the network, including the fullnet, the subnet and the tinynet. The fullnet infers inputs once to generate and transfer its knowledge to the subnet and the tinynet by minimizing the L2-distance between the their stage-wise output feature maps. It is worthy to note that these three networks are weight sharing.}
\label{2}
\end{figure*}
\section{Stage-wise Pruning}
The problem of inaccurate evaluation caused by weight sharing is introduced in \textbf{Sec. 3.1}. We find that the depth and width of the supernet have an impact on training adequacy. Thus, the stage-wise inplace-distillation is proposed in \textbf{Sec. 3.2} to alleviate the aforementioned drawbacks. To efficiently search the optimal subnet from supernet trained by the stage-wise inplace-distillation, we present a distributed evolutionary algorithm in \textbf{Sec. 3.3.}
\subsection{Challenge of Weight sharing}
\label{sub:challengeofweightsharing}
AutoML pruning methods always need to train a supernet that shares weights for all subnets firstly and almost immediately evaluates the accuracy for each subnet. For many AutoML pruning approaches \cite{b8,b9,b12}, pruning candidates(subnets) directly compare with each other according to evaluation accuracy. The subnets with higher evaluation accuracy are selected and expected to deliver high accuracy after training from the scratch. However, such an intention cannot be necessarily achieved since some subnets which have outstanding performance on shared parameters perform poorly after training from scratch. \\
To visualize the performance drop of weight-sharing, we train a supernet with different candidate numbers. For a trained supernet, we randomly sample a batch of subnets from the supernet and evaluate them on the validation dataset. Statistical accuracy expectation $E(a_{super})$ is utilized to evaluate whether the supernet is adequately trained. $E(a_{super})$ is written as
\begin{equation}
  E(a_{super})= \sum _ { i = 1 } ^ { n }  a_{sub_i}
\end{equation}
where $n$ denotes the number of the randomly sampled subnet and $a_{sub_i}$ represents the accuracy of $i$th subnet. As shown in Figure \ref{1}(a), with the increase of the candidate number, the top-1 accuracy expectation of supernet dramatically degrades under 100 epochs while it falls slightly under 500 epochs. Moreover, we train three different depth supernets and then calculate their $E(a_{super})$s. It is found that the expectation of top-1 accuracy is related to the depth of the CNN, which is shown in Figure \ref{1}(b). Next, the problem of weight sharing in terms of depth and width is discussed.\\
\textbf{The depth of supernet.} CNNs(i.e. ResNet152\cite{b3}) is designed deep to enhance the representative ability, which exponentially increases the subnet number in pruning. The subnet number $N$ that inherit weights from the supernet can be formulated as
\begin{equation}
  ||N|| = g^L
\end{equation}
where $g$ denotes the candidate number for each convolutional layer and $L$ represents the depth of the CNN. For the channel pruning of a deep CNN, the search space $N$ is always large(e.g., $>30^{50}$). Hence, many subnets get unfull training results due to weight sharing, which leads to the ineffectiveness of evaluation. \\
\textbf{The width of layers.} Some Auto-ML pruning work\cite{b12,b9,b35,b11,b10} train a single neural network executable at different widths as the supernet. There is not only cross-layer weight sharing but also within the layer. In one layer, the parameters of different widths(candidates) are shared. For instance, all parameters of 0.25$\times$(width scaled by 0.25 of the original width) are shared with the half parameters of 0.5$\times$. Each sampling on any training step is independent of each other. Thus, for $n$ training steps, the sampling times expectation of $i$th channel can be formulated as
\begin{equation}
  E(c_i) = (1 - \frac{i-1}{m})n
\end{equation}
where $m$ denotes the channel number in the layer. Therefore, the channels with small indexes can be trained more times, which causes unfair training results. Formally, we consider a common supernet that contains $L$ layers, each with $m$ channels. In pruning, a group of ratio sequence $R$ can be obtained under certain constraint(i.e. FLOPs), where each ratio sequence $r = [c_1, ..., c_L] \in R$. Because of independent sampling in each layer, the sampling times expectation of $r$ can be described as
\begin{equation}
  E(r) = E(c_1)E(c_2)...E(c_L)
\end{equation}
According to the inequality of arithmetic and geometric means, we have
\begin{equation}
  E(r) \le (\frac{\sum_{i=1}^{L}E(c_i)}{L})^L
\end{equation}
Equality holds if and only if $c_1=c_2=...=c_L=c$. Hence,
\begin{equation}
  E(r) \le ((1-\frac{cm-1}{m})n)^L
\end{equation}
That is to say, the subnet with uniform sampling obtains the most training resource under certain constraint, which means that the pruning strategy of previous work always prone to get a uniform sampling model.\\
\subsection{Stage-wise Inplace Distillation for Training}
\label{sub:subsection_name}
As mentioned before, too many candidates in training can lead to ineffective evaluation on searching because of unfull training results. To adequately train the supernet, we divide the supernet $S$ into N stages according to the depth. Hence, the search space of supernet $S$ can be represented as
\begin{equation}
  S = [S_1 ,..., S_{i}, S_{i+1}, ..., S_N]
\end{equation}
where $S_i$ denotes the stage-wise supernet of $i$th stage. Then we can train the supernet by training the stages separately. The learning of the stage $i$ can be formulated as
\begin{equation}
  W_i^* = min_{W_i} L_{train}(W_i, S_{}i; X, Y)
\end{equation}
where X and Y denote the input data and the groud truth labels, respectively. Subsequently, the candidates number for $i$ th stage can be written as
\begin{equation}
  ||S_i|| = g^{L_i}
\end{equation}
where $L_i$ denotes the depth of the $i$th stage and it is smaller than $L$. The search space can be extremely reduced when we train stage-wise supernet in each stage independently. \\
However, internal ground truth in Eq.(8) cannot be obtained directly from the dataset. One solution is using block-wise feature maps generated by a pre-trained network to supervise the training of subnets. However, it is time-consuming to obtain a pre-trained network through training from the scratch in practice(e.g. ResNet50 $>$ 10 GPU days). Besides, the architecture of teacher and student networks has a huge impact on transferring results\cite{b20}. \\
To tackle the above problem, the stage-wise inplace distillation is proposed here. The essential idea behind the inplace distillation\cite{b10} is to transfer knowledge inside a single supernet from the fullnet to a subnet in each training iteration. For an individual convolutional layer, the performance of the wider candidate cannot be worse than the slim one. Because the wider one can achieve the performance of the slim one by learning weights from some unuseful channels to zeros. Therefore, the performance of any candidates is bounded as follows
\begin{equation}
 |y^{f}-y^{f}|\le |y^{f}-y^{r}| \le |y^{f}- y^{s}|
\end{equation}
where $y^{r}=\sum_{i=1}^{r} w_ix_i$ is the aggregated feature, $r$, $s$ and $f$ denote the channel number of random sampled, the smallest and the largest candidates, respectively. The rule in Eq.(10) can also be extended to the whole supernet, which means the performance of the subnet with any width is bounded in the tinynet and the fullnet.\\
In stage-wise inplace distillation, we use the stage-wise representation of the fullnet to supervise the training of subnets. The pipeline of the supervision with stage-wise inplace distillation is shown in Figure \ref{2}. The output $\hat{Y}_{i-1}$ of the $(i-1)$th stage from the fullnet is adopted by the input of the $i$th stage of subnets. To supervise the subnets learning from the fullnet, the following MSE loss is considered as the distillation loss in Figure \ref{3}
\begin{equation}
  L_{train}(\hat{Y}_{i-1}, Y_i) = \frac{1}{F}||Y_i-\hat{Y_i}||_2^2
\end{equation}
where $Y_i$ and $\hat{Y_i}$ denote the output of the subnets and fullnet in $i$th stage, respectively, $F$ is the number of the channels in $Y$. \\
As mentioned in \ref{sub:challengeofweightsharing}, the channels with larger indexes suffer from unfull training. To ensure the sufficient training of each channel, an intuitive approach is overtraining. Given a batch of input images and ground truth labels, we first calculate the task loss(e.g. cross-entropy) and the gradients of fullnet through forward and backward propagation, simultaneously. The stage-wise feature maps of fullnet $\hat{Y} = [\hat{Y}_1, ... \hat{Y}_N]$ are saved. Subsequently, under the supervision of $\hat{Y}$, the distillation loss in Eq.(11) and the gradients of the stage-wise subnet are calculated. Furthermore, in a subnet training process, we train the smallest width(tinynet) to improve the performance of the supernet. Based on this, each channel can be trained at least once in one iteration. The ground truth label has been generated in the fullnet training process. Thus, the training of each stage-wise subnet can be sped up in a parallel way. The detailed algorithm is described in Algorithm \ref{alg:train algorithm}. We use \textbf{multiprocessing} module\footnotemark[1] to parallelize our algorithm.
\footnotetext[1]{The multiprocessing module is applied to start the child process and execute our customized tasks in the child process. https://github.com/python/cpython/tree/3.9/Lib/multiprocessing/}
\begin{algorithm}[ht]
\caption{ Framework of supervision with stage-wise inplace distillation.}
\label{alg:train algorithm}
\begin{algorithmic}[1]
\Require
The fullnet $S$, the stage-wise supernets $[S_1,...,S_N]$ and the dataset $(X,Y)$;
\Ensure
The well-trained stage-wise supernets $[S_1,...,S_N]$;
\State \textbf{for} $t=1,...,T$ \textbf{do}
\State \hspace*{1em} Get next mini-batch of data $x$ and label $y$ from $(X,Y)$
\State \hspace*{1em} Execute fullnet $y'=S(x)$, and save stage-wise feature maps $\hat{X}=[x,...,\hat{Y}_{N-1}]$, $\hat{Y}=[\hat{Y}_1,...,\hat{Y}_N]$
\State \hspace*{1em} Caculate cross entropy loss, $loss = CE(y', y)$
\State \hspace*{1em} Clear gradients, $optimizer.zero\_grad()$
\State \hspace*{1em} Accumulate gradients, $loss.backward()$
\State \hspace*{1em} Randomly sample width for convolutional layers and obtain stage-wise subnets, $S_r=[S_{r1},...,S_{rN}]$
\State \hspace*{1em} Uniformly sample smallest width for convolutional layers and obtain stage-wise tinynets, $S_t=[S_{t1},...,S_{tN}]$
\State \hspace*{1em} \textbf{multiprocessing} $s=[S_r, S_t]$, $x_s=[\hat{X},\hat{X}]$, $y_s=[\hat{Y}, \hat{Y}]$ \textbf{do}
\State \hspace*{1em} \hspace*{1em} Execute subnet, $y'=s(x_s)$
\State \hspace*{1em} \hspace*{1em} Caculate distillation loss, $loss = L(y', y_s)$
\State \hspace*{1em} \hspace*{1em} Accumulate gradients, $loss.backward()$
\State \hspace*{1em} \textbf{end multiprocessing}
\State \hspace*{1em} Update weights, $optimizer.step()$
\State \textbf{end for}\\
\Return trained stage-wise suernets $[S_1,...,S_N]$;
\end{algorithmic}
\end{algorithm}
\begin{figure*}[t]
\centering
\includegraphics[width=0.8\linewidth]{search.png}
\caption{Illustration of the distributed evolutionary algorithm. Both EA Manager and EA are modified from an evolutionary algorithm. Given FLOPs constraint for the whole network, EA Manager is responsible for searching for the best combination of stage-wise FLOPs. The feedback of each FLOPs gene in EA Manager is provided by each EA with searching for the smallest distillation loss under the stage-wise FLOPs constraint.}
\label{3}
\end{figure*}{}
\subsection{Distributed Evolutionary for Searching}
\label{sub:distributed_evolutionary_for_searching}
After the stage-wise supernet is trained, the learning ability of a subnet can be evaluated by its loss in each stage. However, each stage-wise supernet still contains about  $30^{10}$ stage-wise subnets. It is infeasible to evaluate all of them. For previous one-shot pruning methods, randomly sampling, EA-based and RL-based methods have been used to sample sub-models from the trained supernet for further evaluation. The most recent work found that EA can search better model than RL but need to spend more time on searching. For the proposed stage-wise inplace distillation, a novel method is suggested to search the best subnet according to the stage-wise performance under certain constraint.\\
Because of inplace distillation training mentioned above, the EA in Figure \ref{3} is applied to search the best stage-wise subnet that has the smallest distillation loss under given FLOPs constraint. In \cite{b9,b21}, the genes of each stage-wise subnet were encoded with a vector of channel numbers in each layer. Different from the above work, we aim to imitate the behavior of the teacher in each state. Thus, Eq.(11) is used to evaluate the learning ability of each gene. Because the supernet is split into several stages, the search space of an individual EA is shrunk about $10^{60}$ times. Then the top k genes with the lowest loss are selected for generating the new genes with mutation and crossover. The mutation is carried out by changing a proportion of elements in the gene randomly. The crossover means that we randomly recombine the genes in two-parent genes to generate offspring. We can easily enforce the constraint by eliminating the unqualified genes. Through further repeating the top k selection process and new genes generation process for several iterations, the gene that meets constraints while achieving the lowest loss can be obtained. How can the optimal stage-wise constraints be assigned for each stage?\\
To automatically find the best assignment plan of stage-wise constraint, a distributed evolutionary algorithm(DEA) is proposed in this section. The workflow of DEA is revealed in Figure \ref{3}. The EA Manager is also a kind of evolutionary algorithm that provides the strategy of FLOPs constraint for other EAs. Different from EA above, the genes in EA Manager are encoded with a vector of FLOPs constraint in each stage. The evaluation of each gene is the sum of distillation losses given by all stage-wise EAs. Subsequently, the top k genes which are kept generate off-springs in way of mutation and crossover. After several repeating, the optimal stage-wise constraints can be obtained from the top 1 genes. The detail is shown in Algorithm \ref{alg:distributed evolutionary algorithm}. And each stage-wise EA is paralleled to accelerate the searching process. Specifically, we first use the teacher network to generate a batch of representation features for each stage. Therefore, each EA can search for the best stage-wise subnet independently. After searching all of the stages, we can assemble the best model by selecting the best stage-wise subnet from each stage.
\begin{algorithm}
\caption{ Framework of distributed evolutionary algorithm.}
\label{alg:distributed evolutionary algorithm}
\begin{algorithmic}[1]
\Require
The constraint $C$, the fullnet $S$, The stage-wise supernets $[S_1,...,S_N]$ and the dataset $(X,Y)$;
\Ensure
The best subnet:$S_{top}$;
\State Execute fullnet and save stage-wise feature maps $\hat{Y}$, $y'=S(X)$, $\hat{Y}=[\hat{Y}_1,...,\hat{Y}_N]$
\State Randomly generate a batch of genes $G$ under constraint $C$, $G=[G_1,...,G_N], s.t. ||G_i||=||C_{i1},...,C_{iN}||=C$
\State \textbf{for} $t=1,...,T$ \textbf{do}
\State \hspace*{1em} \textbf{for} $g=G_1,...,G_N$ \textbf{do}
\State \hspace*{1em} \hspace*{1em} Obtain stage-wise constraint from $g$, $g=[C_1,...,C_N]$
\State \hspace*{1em} \hspace*{1em} \textbf{multiprocessing} $c=C_1,...,C_N$, $x_s=X,...,\hat{Y}_{N-1}$, $y_s=\hat{Y}_1,...,\hat{Y}_N$, $s=S_1,...,S_N$ \textbf{do}
\State \hspace*{1em} \hspace*{1em} \hspace*{1em} Search the best stage-wise subnet $s'$ and caculate distillation loss by $EA$ in Algorithm \ref{alg:evolutionary algorithm}, $(s', L_{p_i})=EA(s,x_s,y_s,c)$
\State \hspace*{1em} \hspace*{1em} \textbf{end multiprocessing}
\State \hspace*{1em} \hspace*{1em}Caculate total loss $L$ for $g$, $L=L_{s_1}+...+L_{s_N}$
\State \hspace*{1em} \textbf{end for}
\State \hspace*{1em} Keep top k genes $G_{topk}$ according to $L$
\State \hspace*{1em} Generate $M$ mutation genes, $G_{mutation}=[G_{m1},...,G_{mM}]$, $s.t. ||G_{mi}||=C$
\State \hspace*{1em} Generate $H$ crossover genes, $G_{crossover}=[G_{c1},...,G_{cH}]$, $s.t. ||G_{ci}||=C$
\State \hspace*{1em} Generate new population $G$, $G=G_{mutation}+G_{crossover}$
\State \textbf{end for}
\State Select $S_{top}=[s_1',...,s_k']$ with smallest $L$\\
\Return $S_{top}$;
\end{algorithmic}
\end{algorithm}

\begin{algorithm}
\caption{ Framework of evolutionary algorithm.}
\label{alg:evolutionary algorithm}
\begin{algorithmic}[1]
\Require
The constraint, $C$, the stage-wise supernet $S$, the stage-wise feature maps $(X, Y)$
\Ensure
The best stage-wise subnet:$S_{top}$ and the stage-wise distillation Loss: $L$;
\State Randomly generate a batch of genes $G$ under constraint $C$, $G=[G_1,...,G_d]$
\State \textbf{for} $t=1,...,T$ \textbf{do}
\State \hspace*{1em} \textbf{for} $g=G_1,...,G_d$ \textbf{do}
\State \hspace*{1em} \hspace*{1em} Construct a stage-wise subnet according to $S$ and $g$, $S_g$
\State \hspace*{1em} \hspace*{1em} Calculate the distillation loss of $S_G$, $L_g = L(S_G(X), Y)$, $L$ from Eq.(11)
\State \hspace*{1em} \textbf{end for}
\State \hspace*{1em} Keep top k genes $G_{topk}$ according to $L_g$
\State \hspace*{1em} Generate $M$ mutation genes under constraint $C$, $G_{mutation}=[G_{m1},...,G_{mM}]$
\State \hspace*{1em} Generate $H$ crossover genes under constraint $C$, $G_{crossover}=[G_{c1},...,G_{cH}]$
\State \hspace*{1em} Generate new population $G$, $G=G_{mutation}+G_{crossover}$
\State \textbf{end for}
\State Select $G_{top}$ with smallest $L_g$\\
\Return $G_{top}$, $L_{G_{top}}$;
\end{algorithmic}
\end{algorithm}
\section{Experiments}
\label{sub:experiments}
In this section, the effectiveness of our proposed stage-wise pruning method is demonstrated. First, We explain the experiment settings on CIFAR-10\cite{b33} and ImageNet 2012 dataset\cite{b15}. Then, we prune ResNet\cite{b3} on CIFAR-10 and visualize the consistency of performance between searching and retraining. Moreover, we apply the stage-wise pruning method to ImageNet 2012 and compare the results with those of other state-of-the-art work. Finally, ablation studies are carried out to find out the influence of using inplace distillation.
\subsection{Setups}
\label{sub:setups}
The stage-wise pruning method consists of three steps:\\
\textbf{Stage-wise training.} According to resolution size of feature maps, we split ResNet\cite{b3} and MobileNet series\cite{b1,b2} into 4 and 5 stages, respectively. The distillation loss of each stage can be calculated by Eq.(11) . To match the channel number of the fullnet, the output of each stage is connected with a $1\times 1$ convolutional layer without BatchNorm and non-linear activation. As MetaPruning\cite{b9}, the width of each convolutional layer is subdivided into 31 ratios from 0.1 to 1.0.\\
On CIFAR-10\cite{b33} dataset, we randomly sample 200 images for each class from training images as validation dataset. The remaining images are used to train the supernet. We use momentum SGD to optimize the weights, with initial learning rate $\eta = 0.025$, momentum $0.9$, and weight decay $3 \times 10^{-4}$. The supernet is trained for 50 epochs with batch size 512 and the learning rate decays $0.1 \times$ per 10 epochs.\\
On ImageNet 2012\cite{b15} dataset, we randomly sample 50 images for each class from training images as validation dataset. The remaining images are used to train the supernet. We use momentum SGD to optimize the weights, with initial learning rate $\eta = 0.1$, momentum $0.9$, and weight decay $3 \times 10^{-4}$. The supernet is trained for 100 epochs with batch size 512 and the learning rate decays $0.1 \times$ when the epoch is 30, 60 or 90.\\
\textbf{Stage-wise searching.} After training the stage-wise supernet as above, the best subnet is searched in each stage. Firstly, we use the fullnet to generate and save the stage-wise feature maps with 2048 batch size. Subsequently, the hyperparameter of each EA and EA Manager is set to 128 population number, 0.1 mutation probability, 10 iterations. We use 4 and 5 multiprocess to speed up the searching for ResNet and MobileNet series, respectively. Each process can use 2 GPUs to infer with 2048 batch size. \\
\textbf{Retraining} After searching the best subnet, we adopt the same training scheme as \cite{b9} on ImageNet 2012 for both ResNet and MobileNet series. The same lines of \cite{b23} is followed for the training scheme of ResNet on CIFAR-10. It is noted that all baseline models are trained under the same scheme mentioned above.
\subsection{Pruning ResNet on CIFAR-10 and Analysis}
\label{sub:pruning_on_cifar_10}
To demonstrate the effectiveness of stage-wise pruning, we prune ResNet-56\cite{b3} under 50\% FLOPs constraint on the small dataset of CIFAR-10. As shown in Table \ref{cifar10_results}, our stage-wise pruning method surpasses the baseline model by about 1.4\%. Moreover, our method outperforms all other pruning methods in terms of top-1 accuracy.\\
\begin{table}
\caption{Pruning results of ResNet-56.}
\label{cifar10_results}
\begin{center}
\begin{tabular}[c]{c|c|c}
\hline
Method&FLOPs(M)&Top-1 Acc(\%)\\
\hline
ResNet-56&125.49&93.27\\
FP\cite{b51}&90.90&93.06\\
RFP\cite{b52}&90.70&93.12\\
HRank\cite{b53}&88.72&93.52\\
EagleEye\cite{b54}&62.23&94.66\\
\textbf{SP(Ours)}&61.36&\textbf{95.03}\\
\hline
\end{tabular}
\end{center}
\end{table}
\begin{figure}
\centering
\includegraphics[width=0.9\linewidth]{rank.png}
\caption{Comparison of ranking effectiveness for Stage-wise Pruning(ours), MetaPruning\cite{b9} and AutoSlim\cite{b12}.}
\label{4}
\end{figure}
To evaluate the consistency of model ranking abilities for our method and other AutoML methods, we visualize the relationship between the proxy performance and actual performance. A PruningNet\cite{b9} and a Universally Slimming Network(USNet)\cite{b10} are trained as supernets under the same training scheme due to fairly compare with MetaPruning \cite{b9} and AutoSlim \cite{b12}. The total distillation loss is viewed as the proxy performance of our method. The other two methods take the top-1 accuracy of each subnet that inherits weights from supernet as proxy performance. Each subnet will be trained from scratch in order to obtain its actual performance. As shown in Figure \ref{4}, the proposed method has a strong correlation between the proxy performance and the actual performance while others barely rank the subnets.
\begin{figure*}[htbp]
\subfigure[]{
\begin{minipage}[t]{0.3\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{resnet50_top.png}
\end{minipage}}
\subfigure[]{
\begin{minipage}[t]{0.3\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{resnet50_mid.png}
\end{minipage}}
\subfigure[]{
\begin{minipage}[t]{0.3\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{resnet50_bottem.png}
\end{minipage}}
\caption{The pruning results of ResNet-50. ResNet-50 is stacked by many blocks which consist of three convolutional layers in the main branch. According to the location, we simply divide the three convolutional layers in each block into top layers, middle layers and bottom layers. (a) The number of channels in top layers. (b) The number of channels in the middle layers. (c) The number of channels in the bottom layers.}
\label{5}
\end{figure*}
\subsection{Pruning MobileNet and ResNet on ImageNet 2012}
\label{sub:pruning_mobilenet_and_resnet}
\begin{figure*}[htbp]
\subfigure[]{
\begin{minipage}[t]{0.3\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{tradeoffres.png}
\end{minipage}}
\subfigure[]{
\begin{minipage}[t]{0.3\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{tradeoffv1.png}
\end{minipage}}
\subfigure[]{
\begin{minipage}[t]{0.3\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{tradeoffv2.png}
\end{minipage}}
\caption{The accuracy-FLOPs tradeoff of Auto-ML pruning methods on (a) ResNet-50, (b) MobileNetV2, (c) MobileNetV1. All models are pruned on ImageNet 2012.}
\label{6}
\end{figure*}
In addition, our method is extended to lighting models on a large-scale dataset, ImageNet 2012. Table \ref{imagenet_results} summarizes our results on MobileNet V1\cite{b1}, MobileNet V2\cite{b2} and ResNet-50\cite{b3}. It is noted that we experiment with both residual and nonresidual networks. We compare our results with uniformly pruned baselines and other recent channel pruning methods. It is shown that our method achieves the best results across different computational budgets. In the case of extreme pruning(i.e. 40M FLOPs), MobileNetV1/MobileNetV2 pruned by SWP outperform baseline model to a considerable degree(9.7\% and 6.4\%). Figure \ref{5} compares the curve of top-1 accuracy and FLOPs for the most recent Auto-ML pruning methods and uniformly pruning method. Our SWP models can achieve better accuracy with lower computation complexity than other methods.\\
\begin{table}
\caption{Results of ImageNet classification. We show the top-1 accuracy of each method under the same or closed FLOPs.}
\label{imagenet_results}
\begin{center}
\begin{tabular}[c]{c|ccc}
\hline
Network&Method&Acc@1&FLOPs\\
\hline
&Baseline$0.75\times$ &68.4\%&325M\\
&AMC\cite{b8} &70.5\%&285M\\
&SN\cite{b35} &69.5\%&325M\\
&MP\cite{b9} &70.4\%&281M\\
&AutoSlim\cite{b12} &69.1\%&325M\\
MobileNet V1&\textbf{SWP(ours)}&\textbf{70.9\%}&285M\\
\cline{2-4}
&Baseline$0.25\times$ &50.6\%&41M\\
&MP\cite{b9}&57.2\%&41M\\
&SN\cite{b35}&53.1\%&41M\\
&\textbf{SWP(ours)}&\textbf{60.3\%}&41M\\
\hline
\hline
&Baseline$0.75\times$&69.8\%&220M\\
&AMC \cite{b8}&70.8\%&220M\\
&MP\cite{b9} &71.2\%&220M\\
&SN\cite{b35}&68.9\%&209M\\
&AutoSlim\cite{b12}&72.4\%&207M\\
MobileNet V2&SWP(ours)&\textbf{73.4}\%&220M\\
\cline{2-4}
&Baseline$0.35\times$ &54.3\%&43M\\
&MP\cite{b9} &58.3\%&43M\\
&\textbf{SWP(ours)}&\textbf{60.7\%}&43M\\
\hline
\hline
&Baseline 1.0$\times$&76.6\%&4.1G\\
&Baseline 0.75$\times$ &74.8\%&2.3G\\
&SN\cite{b35} &74.9\%&2.3G\\
&MP\cite{b9}&75.4\%&2.3G\\
&AutoSlim\cite{b12}&75.6&2.3G\\
&AOFP-C1\cite{b43} &75.63\%&2.58G\\
&C-SGD-50\cite{b42} &74.54\%&1.7G\\
ResNet-50&ThiNet-50\cite{b41} &74.7\%&2.1G\\
&\textbf{SWP(ours)}&\textbf{76.1\%}&2.0G\\
\cline{2-4}
&Baseline 0.5$\times$&72.0\%&1.0\\
&SN\cite{b35}&72.5\%&1.0G\\
&ThiNet-30\cite{b41} &72.1\%&1.2G\\
&MP\cite{b9}&73.4\%&1.0G\\
&AutoSlim\cite{b12}&74.0\%&1.0G\\
&\textbf{SWP(ours)}&\textbf{75.6\%}&1.0G\\
\hline
\end{tabular}
\end{center}
\end{table}
\subsection{Pruning under Latency}
More and more attention is paid to directly optimize the inference time on the target device. Without knowing the implementation details inside the device, SWP learns to prune channels according to the latency estimated from the device. To evaluate the realistic acceleration, we measure the forward time of the search for the best subnet on one 2080Ti GPU under the latency constraint. The results of MobileNet V1/V2 are shown in \textbf{\ref{inference-v1}}/\textbf{\ref{inference-v2}}. Under the same compression ratio and similar inference time, our method can obtain better top-1 accuracy.\\
\begin{table}
\caption{Comparison of the realistic inference time with MobileNet V1.}
\label{inference-v1}
\begin{center}
\begin{tabular}[c]{c|c|c}
\hline
Ratio&Baseline/Pruned time(ms)&Baseline/Pruned Acc@1\\
\hline
$1\times$&0.65/-&70.9\%/-\\
$0.75\times$&0.51/0.51&68.4\%/70.2\%\\
$0.5\times$&0.35/0.35&63.3\%/67.5\%\\
$0.25\times$&0.21/0.20&49.8\%/59.9\%\\
\hline
\end{tabular}
\end{center}
\end{table}
\begin{table}
\caption{Comparison of the realistic inference time with MobileNet V2.}
\label{inference-v2}
\begin{center}
\begin{tabular}[c]{c|c|c}
\hline
Ratio&Baseline/Pruned time(ms)&Baseline/Pruned Acc@1\\
\hline
$1\times$&0.72/-&71.7\%/-\\
$0.75\times$&0.53/0.52&69.8\%/71.2\%\\
$0.5\times$&0.41/0.41&65.4\%/69.3\%\\
$0.35\times$&0.33/0.32&60.3\%/64.8\%\\
\hline
\end{tabular}
\end{center}
\end{table}
\subsection{Combination with NAS}
We extend our effective method to NAS. Many NAS work\cite{b23,b24,b25} only search operations without considering the channel number. Here, we simultaneously take kernel size, layers and channels into account. Our network architectures are designed on MobileNet-V3\cite{b55} that consists of a stack with inverted bottleneck residual blocks(MB-Conv). The search space contains the channel number, kernel size and layer number of each stage. The detailed search space is shown in \textbf{\ref{mobilenetv3}}. \\
\begin{table}
\caption{MobileNetV3-based search space. \checkmark denotes the search space of the stage containing SE module\cite{b57}.}
\label{mobilenetv3}
\begin{center}
\begin{tabular}[c]{c|c|c|c|c|c}
\hline
Stage&Operator&Ratio&Layers&Kernel Sizes&SE module\\
\hline
&Conv&1.0&1&3&\\
1&MBConv&0.1-1.0&1-2&3&\checkmark\\
2&MBConv&0.1-1.0&2-3&3&\checkmark\\
3&MBConv&0.1-1.0&2-3&3,5&\checkmark\\
4&MBConv&0.1-1.0&2-4&3,5&\checkmark\\
5&MBConv&0.1-1.0&2-6&3,5&\checkmark\\
6&MBConv&0.1-1.0&2-6&3,5&\checkmark\\
6&MBConv&0.1-1.0&1-2&3,5&\checkmark\\
&MBConv&1.0&1&1&\\
\hline
\end{tabular}
\end{center}
\end{table}
Several stage-wise supernets are trained on ImageNet dataset\cite{b15} using the same settings as MobileNet-V2. After the training, the optimal subnet is found based on distributed evolutionary algorithm. The retraining setting follows MobileNet-V3: RMSProp optimizer with decay 0.9 and momentum 0.9; batch size with 4096; initial learning rate 0.1 that decays by 0.01 every 3 epochs; dropout of 0.8 that only be applied in the fullnet. It is worth noting that a $3 \times 3$ kernel is central cropped from a $5 \times 5$ kernel and lower-index layers in each stage are always kept. Therefore, the weights also share on the kernel size and layer size dimension.\\
The performance of ours and previous SOTA methods with various search spaces are revealed in \textbf{\ref{NAS_result}}. Our models have 2.0\% and 1.7\% better top-1 accuracy than MobileNetV3-small and MobileNetV3-large with similar FLOPs, respectively. Compared with other mobile-setting NAS methods, our method outperforms them under wide-range FLOPs constraints because of more flexible search space. Specifically, our large-size model obtains a competitive result under fewer FLOPs compared with DNA\cite{b14} that also split a supernet into several stages but use a pretrained teacher network to supervise. It shows that the architecture of the teacher network has a huge impact on the effectiveness of distillation. We will discuss more details in \textbf{Sec 4.7}.
\begin{table}
\caption{Results of NAS methods on Imagenet.}
\label{NAS_result}
\begin{center}
\begin{tabular}[c]{c|cc}
\hline
Method&Acc@1&FLOPs\\
\hline
MobileNetV3-large\cite{b55}&75.2\%&219M\\
OFA\cite{b27}&76.0\%&230M\\
MNasNet\cite{b56} &75.2\%&315M\\
DNA-a\cite{b14} &77.1\%&348M\\
\textbf{SWPNas-large(ours)}&\textbf{76.9\%}&214M\\
\hline
\hline
MobileNetV3-small\cite{b55}&67.4\%&56M\\
Mnas-small\cite{b56} &64.9\%&65M\\
SWPNas-small&69.4\%&56M\\
\hline
\end{tabular}
\end{center}
\end{table}
\subsection{Visualization of Searched Models}
Channel pruning models are visualized and some insights from the results are discussed. We compare our results with default channels and MetaPruning\cite{b9} on ResNet-50. In Figure \ref{6} (a-c), we show the channel number in the top, middle and bottom layers of bottleneck blocks on ResNet-50, respectively. Firstly, it is found that our method is prone to prune more channels from top layers compared with MetaPruning. It is noted that although top layers have a small number of channels, the output feature maps of the top layer can be extracted by the next middle layer where $kernel\ size = 3$. Hence, prune top layers can reduce computational complexity. Secondly, both our method and Metapruning keep more channels for downsampling layers because the feature map size is shrunk. Moreover, our method prunes fewer channels for bottom layers since the feature maps between the subnet and the fullnet should be as close as possible.\\
\subsection{Ablation Study}
\label{sub:ablation_study}
\textbf{The choice of the teacher network.} The influence of the distillation strategy in ResNet-50 is analyzed. In our method, the teacher network and the student network are jointly trained by stage-wise inplace distillation. In Strategy 1, we use a pre-trained network that has the same architecture as fullnet to supervise the training of subnets. In Strategy 2, EfficientNet-B0\cite{b44} of which performance surpasses ResNet-50 with lower parameters are employed as the teacher network. The results are shown in Table \ref{distillation_strategies}. It is found that the performance of the model searched with inplace distillation method is almost the same as the one searched with a pre-trained method. Hence, the fullnet can supervise the training of the subnets well while training itself. It is unnecessary to spend a lot of time to obtain a pre-trained teacher model. Moreover, despite EfficientNet-B0 has outstanding performance compared with ResNet-50, the models searched with EfficientNet-B0 have worse performance. It may be caused by the large gap between the architectures of the teacher and the student networks.\\
To further figure out the reason why using ResNet-50 to search can achieve better performance of models than EfficientNet-B0, we visualize the channel number in bottom layers in each block. As shown in Figure \ref{7}, the model searched with EfficientNet-B0 keeps fewer channels. And EfficientNet-B0 has much fewer channels than ResNet-50. For example, EfficientNet-B0 has only 40 channels while ResNet-50 has 64 channels in the bottom layer of the first stage. Thus, the student does not need many channels to imitate the stage-wise feature maps generated by EfficientNet-B0. However, over-pruning the channels from the bottom layers will result in insufficient information transmission between the stages.\\
\begin{table}[t]
\caption{Comparison of Stage-wise pruning with different distillation strategies.}
\label{distillation_strategies}
\begin{center}
\begin{tabular}[c]{c|ccc}
\hline
Teacher&Student FLOPs&Acc@1\\
\hline
\hline
Ours&2.0G&76.1\%\\
&1.0G&75.6\%\\
\hline
\hline
Strategy 1&2.0G&76.1\%\\
&1.0G&75.6\%\\
\hline
\hline
Strategy 2&2.0G&75.6\%\\
&1.0G&73.1\%\\
\hline
\end{tabular}
\end{center}
\end{table}
\textbf{The stage number.} In previous experiments, the supernet is simply split into 4(ResNet) or 5(MobileNet series) stages according to downsampling. How many split stages are optimal to transfer knowledge from the teacher to the student? We set the stage number as 8, 4, 2, 1 for ResNet-50. It is seen from Table \ref{memory}, the performance of the model searched with 8 stages is almost the same as the one searched with 4 stages and surpasses the one searched with lower than 4 stages. Besides, since we need to save the intermediate feature maps of the teacher network to supervise the training of the student network, the more memory is required if larger stage number is set.\\
\begin{table}[t]
\caption{Comparison of Stage-wise pruning with different stage numbers. We evaluate the models with pruned ResNet-50 under 2.0G FLOPs constraint. The memory is tested under a 224$\times$224 image.}
\label{memory}
\begin{center}
\begin{tabular}[c]{c|ccc}
\hline
Number of stages&Acc@1&Memory\\
\hline
\hline
8&76.18\%&3M\\
4&76.12\%&1.5M\\
2&75.88\%&0.5M\\
1&74.53\%&0.1M\\
\hline
\end{tabular}
\end{center}
\end{table}
\begin{figure}
\centering
\includegraphics[width=0.9\linewidth]{ablation_study.png}
\caption{The number of channels in bottom layers with different teachers.}
\label{7}
\end{figure}
\subsection{Conclusion}
\label{sub:conclusion}
In this work, we have presented stage-wise channel pruning. A stage-wise training process based on inplace distillation and a distributed evolutionary algorithm is proposed in this paper. It is found that the large search space causes low accuracy of the supernet. Hence, we split a supernet into several stage-wise supernets to degrade the complexity of the search space both in training and searching. The experiments show the effectiveness of our proposed method by delivering higher accuracy than the previous work on both CIFAR-10 and ImageNet dataset. The consistency of the proxy and actual subnet performance is greatly improved. Moreover, we extend our method to NAS and also achieve SOTA on the mobile setting. Experiments with various distillation strategies prove that inplace distillation can replace pre-trained distillation, thereby reducing the time to train a teacher network from scratch. We further discuss the impact of stage numbers on search results and found that splitting supernet according to downsampling is the best tradeoff between memory and accuracy.

\section*{Acknowledgement}
\label{sub:acknowledgement}
We gratefully acknowledge the support of National Key R\&D Program of China (2018YFB1308400) and Natural Science Foundation of Zhejiang Province (NO. LY21F030018).

{\small
\bibliographystyle{ieee_fullname}