Accretion-induced Collapse of Dark Matter-admixed Rotating White Dwarfs: Dynamics and Gravitational-wave Signals

It is widely believed that dark matter (DM) constitutes the major mass-energy component of galaxy clusters (Clowe et al. 2006) and large-scale structures of the universe (Davis et al. 1985). Besides terrestrial experiments, physicists are tackling the DM problem through astrophysical observations. It has been shown that in a region with a high concentration of DM particles, DM could be captured by normal matter (NM; Sulistiyowati & Ibrahim 2014; Arun et al. 2019). Therefore, it is natural to expect stellar objects that are composed of NM and DM. There have been extensive theoretical studies on the possible effects of DM admixture on stellar evolution (Lopes & Lopes 2019; Clea et al. 2020; Raen et al. 2021). Unusual stellar objects consistent with these models might hint at the existence of DM-admixed stars. Furthermore, there are studies utilizing DM-admixed star models to understand the properties of DM. For instance, Leung et al. (2022) proposed a method for inferring the DM particle mass by measuring the tidal deformability of neutron stars. Using the DM-admixed neutron star model, Bramante et al. (2013) and Bell et al. (2013) gave constraints on the bosonic DM particle mass and annihilation cross section. These examples show that DM-admixed stellar objects could be a promising channel to probe astrophysical DM.

The majority of studies on white dwarfs (WDs) assume they are not rotating, but observational evidence shows the opposite (Spruit 1998; Kawaler 2004). It has been suggested that WDs gain angular momentum through accretion from a companion star (Langer et al. 2004; Yoon & Langer 2004) or mergers between two or more WDs (Gvaramadze et al. 2019; Pshirkov et al. 2020). Therefore, rotation is an important ingredient of the full picture of WD structure and evolution (Yoon & Langer 2004, 2005). In addition, rotating WDs have been proposed to be progenitors of superluminous thermonuclear supernovae because rotating WDs can support more mass than their traditional Chandrasekhar limit (Pfannes et al. 2010; Wang et al. 2014; Fink et al. 2018). Recently, the effects of the strong magnetic field on the equilibrium structures of WDs have been studied (Franzon & Schramm 2015; Bera & Bhattacharya 2016; Chatterjee et al. 2017), for which WD rotation plays a critical role.

It is widely believed that a WD undergoes thermonuclear explosion when its mass approaches the canonical Chandrasekhar limit. However, if the WD contains an oxygen-neon (O-Ne) core, accretion-induced collapse (AIC) is possible as its mass increases toward the Chandrasekhar limit through stable accretion from a companion object (Nomoto & Kondo 1991; Wang 2018; Ruiter et al. 2019), though a binary WD merger seems to be another possible scenario (Liu & Wang 2020). The collapse is triggered by electron capture in the degenerate matter, weakening the electron degenerate pressure (Brooks et al. 2017). On the other hand, pycnonuclear burning is also possible in such an extremely dense core. Hence, the ultimate fate of an O-Ne WD would depend on the competition between nuclear runaway and electron capture (Wang & Liu 2020). However, it was later found that the central temperature of O-Ne WDs is insufficient for explosive O-Ne burning (Wu & Wang 2018). Even if deflagration occurs, it fails to unbind the WD, which directly leads to a collapse for a wide range of parameters (Zha et al. 2019b; Leung & Nomoto 2019; Leung et al. 2020).

Besides the iron core collapse of massive stars, the AIC of WDs has been proposed as another channel for forming neutron stars. However, AIC is much less luminous than typical core-collapse supernovae. The small amount of nickel synthesized indicates that AICs are usually faint transients (Darbha et al. 2010). On the other hand, AIC emits radio signatures (Moriya 2016) and has been hypothesized as a source candidate for fast radio bursts (Margalit et al. 2019) and millisecond pulsars (Wang et al. 2022). Electromagnetic-wave detection of AIC would be a challenging but possible task. One possible way to search for AIC is by neutrino detection because a neutrino burst should accompany AIC after the WD dynamical collapse (Dar et al. 1992). The burst luminosity could be as large as 10⁵⁵ erg s⁻¹ (Dessart et al. 2006; Zha 2019). On the other hand, the collapse dynamics of the compact iron core are expected to produce strong gravitational-wave (GW) signals (Ott et al. 2005; Ott 2009). There have been some efforts to predict the GW signature from an AIC. Dessart et al. (2006) simulated 2D AIC with neutrino transport and estimated the GW emission from AIC via the Newtonian quadrupole formula. They concluded that LIGO-class detectors could detect Milky Way AIC events. Abdikamalov et al. (2010) found that the GW signals from an AIC show a generic "Type III" shape, though detailed neutrino physics have been omitted.

Although DM-admixed neutron stars have been studied and applied to explain anomalous compact objects, there is still no in-depth research on their formation channel. Even though Leung et al. (2019) and Zha et al. (2019a) numerically investigated DM-admixed AIC (DMAIC), they assumed the DM component to be spherically symmetric and nonmoving. As pointed out by Leung et al. (2019), the stationary DM approximation may break down if the dynamical timescales for DM and NM become comparable, and the dynamical modeling of the DM becomes important. They also pointed out that there is a moment during the collapse in which the NM has a mass density comparable to that of the DM. Also, Chan (2021) showed that fermionic DM with sub-gigalectronvolt particle mass would produce a massive and extended component comparable in size to that of the NM. In such a scenario, modeling the DM dynamics would be necessary. In this study, we extend the multidimensional simulations by Zha et al. (2019a) to include also the dynamical evolution of the DM component. Our study aims to investigate if DMAIC could make a DM-admixed neutron star, with the DM motion taken into account, and to predict the corresponding GW signature to facilitate the search for DM through observing AIC in the future. The structure of this paper is as follows: we present the method used to obtain the progenitor and the tools used for the simulation in Section 2. We then present the results of the collapse dynamics and GW signals in Section 3. Finally, we conclude our study in Section 4.

We compute DM-admixed rotating WDs (DMRWDs) as DMAIC progenitors by solving the Newtonian hydrostatic equations, including the centripetal force

$$\begin{eqnarray}&&\begin{array}{l}\vec{{\rm{\nabla }}}{P}_{1}=-{\rho }_{1}\vec{{\rm{\nabla }}}{\rm{\Phi }},\\ \vec{{\rm{\nabla }}}{P}_{2}=-{\rho }_{2}\vec{{\rm{\nabla }}}{\rm{\Phi }}+[{\rho }_{2}{\omega }_{2}{\left({\rm{s}}\right)}^{2}{\rm{s}}]\hat{{\rm{s}}},\\ {{\rm{\nabla }}}^{2}{\rm{\Phi }}=4\pi G({\rho }_{1}+{\rho }_{2}).\end{array}\end{eqnarray} \tag{ 1 }$$

Here, the subscript i = 1(2) denotes the DM (NM) quantities, and ρ, P, ω, and Φ are the density, pressure, angular speed, and gravitational potential of the fluid element. s is the perpendicular distance from the rotation axis, and $\hat{{\rm{s}}}$ is the unit vector orthogonal to and pointing away from that axis. The angular velocity is assumed to be a function of s only. We consider the Newtonian framework because the rotation speed and compactness of WDs are low.

We follow Eriguchi & Mueller (1985), Hachisu (1986), and Aksenov & Blinnikov (1994) to integrate the equations of equilibrium:

$$\begin{eqnarray}&&\begin{array}{l}{H}_{i}+{\rm{\Phi }}+{\delta }_{i2}{h}_{i}^{2}{\psi }_{i}={C}_{i},\\ \displaystyle \int \displaystyle \frac{{{dP}}_{i}}{{\rho }_{i}}={H}_{i},\\ \displaystyle \int \omega {\left({\rm{s}}\right)}_{i}^{2}{\rm{s}}{ds}=-{h}_{i}^{2}{\psi }_{i},\end{array}\end{eqnarray} \tag{ 2 }$$

where C_i is an integration constant, H is the enthalpy, ψ is the rotational potential, and h² is a constant to be determined (Hachisu 1986). We solve the equilibrium equations for the DM and NM using a two-fluid, self-consistent field method (Chan 2022).

We have considered rotation profiles for the NM from Hachisu (1986) and Yoshida (2019) including (1) the rigid rotation

$$\begin{eqnarray}&&\omega {\left({\rm{s}}\right)}_{2}^{2}={{\rm{\Omega }}}_{2}^{2},\end{eqnarray} \tag{ 3 }$$

and (2) the "Kepler" profile

$$\begin{eqnarray}&&\omega ({\rm{s}}{)}_{2}^{2}\propto 1/{\left({d}^{3/2}+{{\rm{s}}}^{3/2}\right)}^{2},\end{eqnarray} \tag{ 4 }$$

which resembles a rapidly rotating core surrounded by an envelope rotating at its Keplerian limit. Here, d is the rotating core radius.

We integrate the angular velocity to obtain the effective potential of the rigid rotation

$$\begin{eqnarray}&&{\psi }_{2}={{\rm{s}}}^{2}/2.\end{eqnarray} \tag{ 5 }$$

The effective potential for the Kepler rule is

$$\begin{eqnarray}\begin{array}{rcl}{\psi }_{2} & = & -\displaystyle \frac{1}{9}\left[-\displaystyle \frac{6\sqrt{{\rm{s}}}}{{{\rm{s}}}^{\tfrac{3}{2}}+{d}^{\tfrac{3}{2}}}+\displaystyle \frac{1}{d}\mathrm{ln}\left(\displaystyle \frac{{\left(\sqrt{d}+\sqrt{{\rm{s}}}\right)}^{2}}{d+{\rm{s}}-\sqrt{{sd}}}\right)\right.\\ & & \left.-\,\displaystyle \frac{2\sqrt{3}}{d}{\tan }^{-1}\left(\displaystyle \frac{1-2\sqrt{s/d}}{\sqrt{3}}\right)\right].\end{array}\end{eqnarray} \tag{ 6 }$$

To simulate DMAIC, we solve the two-dimensional Euler equations assuming axial symmetry:

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}{\rho }_{i}+{\rm{\nabla }}\cdot (\alpha {\rho }_{i}{{\boldsymbol{v}}}_{i})=0,\\ {\partial }_{t}({\rho }_{i}{{\boldsymbol{v}}}_{i})+{\rm{\nabla }}\cdot [\alpha {\rho }_{i}({{\boldsymbol{v}}}_{i}\otimes {{\boldsymbol{v}}}_{i})]+{\rm{\nabla }}(\alpha {P}_{i})\\ =-\,\alpha ({\rho }_{i}-{P}_{i}){\rm{\nabla }}{\rm{\Phi }}.\end{array}\end{eqnarray} \tag{ 7 }$$

Here, $\alpha =\exp (-\phi /{c}^{2})$ is the lapse function with c being the speed of light. It is used to mimic general relativistic time-dilation effects and has been applied to study the first-order quantum chromodynamics phase transition in core-collapse supernovae (Zha et al. 2020). We also solve the advection equation for the NM total internal energy density ${\tau }_{2}={\rho }_{2}{\epsilon }_{2}+{\rho }_{2}{v}_{2}^{2}/2$ and electron fraction Y_e :

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}{\tau }_{2}+{\rm{\nabla }}\cdot [\alpha ({\tau }_{2}+{P}_{2}){{\boldsymbol{v}}}_{2}]=-\alpha {\rho }_{2}{{\boldsymbol{v}}}_{2}\cdot {\rm{\nabla }}{\rm{\Phi }},\\ {\partial }_{t}({\rho }_{2}{Y}_{e})+{\rm{\nabla }}\cdot (\alpha {\rho }_{2}{Y}_{e}{{\boldsymbol{v}}}_{2})=0.\end{array}\end{eqnarray} \tag{ 8 }$$

The gravitational potential Φ is solved by a multipole solver, for which we adopt the one by Couch et al. (2013), which can reduce error by computing the potential at the cell center while using the mass density at that point. To mimic general relativistic strong-field effects, we use the modified case A potential (Muller et al. 2008) as an additional correction to the Newtonian potential:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}\to {\rm{\Phi }}-\langle {\rm{\Phi }}\rangle +{{\rm{\Phi }}}_{\mathrm{TOV},1}+{{\rm{\Phi }}}_{\mathrm{TOV},2},\\ \langle {\rm{\Phi }}\rangle =-4\pi {\displaystyle \int }_{0}^{\infty }{dr}^{\prime} r{{\prime} }^{2}\displaystyle \frac{\langle {\rho }_{1}+{\rho }_{2}\rangle }{| r-r^{\prime} | }.\end{array}\end{eqnarray} \tag{ 9 }$$

Here, r is the radial distance and 〈ρ₁ + ρ₂〉 represents the angular average of the total density. Φ_TOV,i for i = 1, 2 are the relativistic corrections:

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{\mathrm{TOV},i}=-4\pi {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{dr}^{\prime} }{r{{\prime} }^{2}}\displaystyle \frac{1}{{{\rm{\Gamma }}}_{i}^{2}}\left(\displaystyle \frac{{m}_{\mathrm{TOV},{\rm{i}}}}{4\pi }+r{{\prime} }^{3}{P}_{i}\right)\\ \left(1+{\epsilon }_{i}+\displaystyle \frac{{P}_{i}}{{\rho }_{i}}\right),\\ {m}_{\mathrm{TOV},{\rm{i}}}=4\pi {\displaystyle \int }_{0}^{r}{dr}^{\prime} r{{\prime} }^{2}{{\rm{\Gamma }}}_{i}{\rho }_{i}(1+{\epsilon }_{i}),\\ {{\rm{\Gamma }}}_{i}=\sqrt{1+{v}_{r,i}^{2}-\displaystyle \frac{2{m}_{\mathrm{TOV},{\rm{i}}}}{r}},\end{array}\end{eqnarray} \tag{ 10 }$$

where v_r,i is the radial velocity. We adopt a finite-volume approach to solve the hydrodynamic equation in spherical coordinates (Mignone 2014). We use the piecewise parabolic method (Colella & Woodward 1984) to reconstruct primitive variables at the cell interface and the HLLE Riemann solver (Toro 2009) to compute fluxes across cell boundaries. The reconstruction and flux evaluation are done on a dimension-by-dimension basis. We discretize the temporal evolution using the method of lines where the strong stability-preserving five-step, fourth-order Runge–Kutta method is implemented (Gottlieb et al. 2011). In addition to the (modified) Euler equation, we also append the internal energy equation for the NM:

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}({\rho }_{2}{\epsilon }_{2})+{\rm{\nabla }}\cdot \{\alpha [({\rho }_{2}{\epsilon }_{2})+{P}_{2}]{{\boldsymbol{v}}}_{2}\}\\ =\,{{\boldsymbol{v}}}_{2}\cdot {\rm{\nabla }}(\alpha {P}_{2})-\alpha {P}_{2}({{\boldsymbol{v}}}_{2}\cdot {\rm{\nabla }}{\rm{\Phi }}).\end{array}\end{eqnarray} \tag{ 11 }$$

It not only allows one to interpolate the internal energy density to the cell interface so that computational cost could be reduced but also reduces the error of due to advection (Zingale et al. 2020). We adopt a computational grid similar to that in Skinner et al. (2016), in which an analytic function describes the positions of the radial cell interfaces as

$$\begin{eqnarray}&&{r}_{i}={A}_{t}\sinh ({x}_{t}i/{A}_{t}).\end{eqnarray} \tag{ 12 }$$

Here, i is the cell index. We set x_t = 0.5 (in code unit) ⁴  and A_t = 150 so that a central resolution of around 0.74 km is provided, while a total of 500 computational grids are used to contain the progenitor. We use 20 grids to resolve the polar direction, which we find sufficient for ensuring convergence in GW signals for both the NM and DM.

After mapping the density profiles of the NM and DM components computed from Section 2.1, we assign an initial temperature profile to the NM (Dessart et al. 2006):

$$\begin{eqnarray}&&T(\rho )={T}_{c}{\left(\rho /{\rho }_{c}\right)}^{0.35}.\end{eqnarray} \tag{ 13 }$$

Here, T_c = 10¹⁰ K and ρ_c = 5 × 10¹⁰ g cm⁻³ are the central temperature and density, respectively. The core electron capture process initiates the AIC. We implement the parameterized electron capture scheme described in Liebendorfer (2005) to simulate such a process. In their work, Y_e depends on ρ₂ as

$$\begin{eqnarray}\begin{array}{rcl}x({\rho }_{2}) & = & \max \left[-1,\min \left(1,\displaystyle \frac{2-\mathrm{log}{\rho }_{2}-\mathrm{log}{\rho }_{\alpha }-\mathrm{log}{\rho }_{\beta }}{\mathrm{log}{\rho }_{\alpha }-\mathrm{log}{\rho }_{\beta }}\right)\right],\\ & & {Y}_{e}(x)=\displaystyle \frac{1}{2}({Y}_{b}+{Y}_{a})+\displaystyle \frac{x}{2}({Y}_{b}-{Y}_{a})\\ & & +\,{Y}_{c}[1-| x| +4| x| (| x| -1/2)(| x| -1)].\end{array}\end{eqnarray} \tag{ 14 }$$

Here, $\mathrm{log}{\rho }_{\alpha }$, $\mathrm{log}{\rho }_{\beta }$, Y_a , Y_b , and Y_c are fitting parameters and are obtained by Leung et al. (2019) and Zha et al. (2019a). We first assign an initial equilibrium Y_e profile to the NM using Equation (14). We then start the electron capture process by updating Y_e at each time step using the same equation. We force Y_e to strictly decrease with time. We terminate the electron capture process once the core bounce condition (Liebendorfer 2005) is achieved, which is when the core NM entropy is larger than 3 k_B , where k_B is the Boltzmann constant.

We use the quadrupole formula in the weak-field approximation to compute the GW strain (Finn & Evans 1990; Moenchmeyer et al. 1991):

$$\begin{eqnarray}&&{h}_{+}=\displaystyle \frac{3}{2}\displaystyle \frac{G}{{{Dc}}^{4}}{\sin }^{2}\theta \displaystyle \frac{{d}^{2}}{{{dt}}^{2}}{I}_{{zz}}.\end{eqnarray} \tag{ 15 }$$

Here, D = 10 kpc is the assumed distance and θ is the orientation angle of the collapsing DMRWD, and I_zz is the moment of the inertia tensor:

$$\begin{eqnarray}&&{I}_{{zz}}=\displaystyle \frac{1}{3}{\int }_{\mathrm{All}\ \mathrm{Space}}({\rho }_{1}+{\rho }_{2}){r}^{2}{P}_{2}(\cos \theta )d\tau .\end{eqnarray} \tag{ 16 }$$

To simulate AICs, we first use the ideal degenerate Fermi gas equation of state (EOS) for equilibrium structure construction. Following the subsequent collapse dynamics, we use the nuclear matter EOS given by Shen et al. (2011), widely used in simulating core-collapse supernovae and neutron star dynamics. We adopt the ideal degenerate Fermi gas EOS for the DM component (Narain et al. 2006).

We have computed a series of DMRWD models as DMAIC progenitors. The stellar parameters of these progenitors have been listed in Table 1 for reference. The progenitors have DM mass fractions _DM = M_DM/(M_DM + M_NM), where M_NM (M_DM) is the NM (DM) mass, from 0.01–0.2 and include rigidly rotating and differentially rotating DMRWDs with different d as described in Equation (4). In particular, d is chosen so that ${\rho }_{2}(r=d,\theta =\tfrac{\pi }{2})={\alpha }_{d}{\rho }_{2c}$. We choose α_d = 0.1 and 0.01. Another free parameter to be specified for these progenitors is the central angular velocity Ω_c . We adjust this value for rigidly rotating DMRWDs so that the corresponding pure NM progenitor almost rotates at the Keplerian limit and that a total mass of ≈1.8 M_⊙ is achieved for a pure NM, differentially rotating WD. We fix the DM particle mass to be 0.1 GeV for all of these progenitors. As shown in Chan (2021), the fluid component formed by DM particles with such a mass will be more diffusive and comparable in size to that of the NM.

Table 1. Stellar Parameters for Different DMAIC Progenitors

Model	MNM	MDM	αd	log10 ρ1c	Ωc	ReNM	ReDM	DM	tb	log10 ρ1b	log10 ρ2b	MPN
	(M⊙)	(M⊙)		(g cm−3)	(s−1)	(km)	(km)		(ms)	(g cm−3)	(g cm−3)	(M⊙)
Rigid-NM	1.477	0.000	⋯	⋯	10.8	1105	⋯	0	53.151	⋯	14.318	1.217
Rigid-0.01	1.447	0.015	⋯	8.816	10.8	1098	400	0.01	53.424	11.078	14.317	1.194
Rigid-0.03	1.416	0.044	⋯	8.980	10.8	1027	568	0.03	53.717	11.093	14.315	1.170
Rigid-0.05	1.397	0.074	⋯	9.046	10.8	980	695	0.05	53.902	11.097	14.316	1.157
Rigid-0.07	1.384	0.104	⋯	9.086	10.8	948	807	0.07	54.035	11.100	14.315	1.146
Rigid-0.09	1.374	0.136	⋯	9.114	10.8	923	904	0.09	54.139	11.102	14.315	1.139
Rigid-0.1	1.307	0.152	⋯	9.125	10.8	910	954	0.1	54.183	11.103	14.315	1.136
Rigid-0.2	1.343	0.336	⋯	9.193	10.8	857	1379	0.2	54.496	11.107	14.313	1.119
Kepler-NM-d001	1.770	0.000	0.01	⋯	32.5	1826	⋯	0	35.124	⋯	14.355	1.597
Kepler-0.01-d001	1.725	0.017	0.01	8.853	32.5	1766	420	0.01	35.352	11.073	14.352	1.553
Kepler-0.03-d001	1.662	0.051	0.01	9.016	32.5	1494	587	0.03	35.647	11.086	14.351	1.498
Kepler-0.05-d001	1.623	0.085	0.01	9.082	32.5	1343	710	0.05	35.835	11.091	14.352	1.463
Kepler-0.07-d001	1.596	0.120	0.01	9.122	32.5	1256	812	0.07	35.969	11.094	14.350	1.439
Kepler-0.09-d001	1.577	0.156	0.01	9.150	32.5	1190	904	0.09	36.072	11.096	14.350	1.419
Kepler-0.1-d001	1.569	0.174	0.01	9.162	32.5	1159	948	0.1	36.116	11.097	14.350	1.411
Kepler-0.2-d001	1.518	0.379	0.01	9.232	32.5	1020	1352	0.2	36.424	11.101	14.349	1.362
Kepler-NM-d001	1.771	0.000	0.1	⋯	45.2	1106	⋯	0	32.311	⋯	14.354	1.598
Kepler-0.01-d01	1.727	0.017	0.1	8.860	45.2	1098	417	0.01	32.555	11.070	14.353	1.555
Kepler-0.03-d01	1.677	0.052	0.1	9.026	45.2	1062	579	0.03	32.788	11.084	14.354	1.511
Kepler-0.05-d01	1.647	0.087	0.1	9.094	45.2	1034	700	0.05	32.924	11.088	14.351	1.483
Kepler-0.07-d01	1.625	0.122	0.1	9.135	45.2	1007	801	0.07	32.024	11.092	14.351	1.460
Kepler-0.09-d01	1.609	0.159	0.1	9.164	45.2	987	892	0.09	33.101	11.095	14.352	1.446
Kepler-0.1-d01	1.602	0.178	0.1	9.176	45.2	980	941	0.1	33.133	11.096	14.352	1.439
Kepler-0.2-d01	1.556	0.389	0.1	9.247	45.2	916	1343	0.2	33.364	11.101	14.355	1.396

Note. They include rigid (labeled Rigid) and differentially (labeled Kepler) rotating DMRWDs. All progenitors have an NM central density of 5 × 10¹⁰ g cm⁻³. The DM particle mass is 0.1 GeV. In this table, R_eNM (R_eDM) is the equatorial radius of the progenitor for the NM (DM) component. ρ_1c is the DM central density, _DM is the DM fraction, and t_b is the bounce time. ρ_2b (ρ_1b) is the maximum NM (DM) density at the core bounce. M_PNS is the proto-neutron star mass, defined as summing all the NM mass with ρ₂ > 10¹¹ g cm⁻³ at the end of the simulation.

Download table as:  ASCIITypeset image

3.1.1. Collapse Dynamics

We first focus on the collapse dynamics of DMAIC. From Table 1, we observe that the admixture of DM delays the time of core bounce and reduces the proto-neutron star mass, which is similar to the results by Leung et al. (2019) and Zha et al. (2019a). We show the maximum NM density evolution for the rigidly rotating DMAIC models in the right panel of Figure 1. Despite having different initial and proto-neutron star masses (see Table 1), the maximum NM density evolution is almost identical for all DMAIC models. The final maximum NM densities are also insensitive to the DM mass fraction _DM. Furthermore, we find that AIC is successful for all DMAIC progenitors. This differs from the results presented by Leung et al. (2019) and Zha et al. (2019a) because we assume different DM particle masses. Their work assumed a heavy (1 GeV) DM particle mass, leading to a more compact DM core with a large central density. Thus, it significantly impacts the NM density profile near its center. The NM density decreases sharply due to the strong gravitational force provided by the compact DM core. Electron capture is less efficient in their model, so the NM component's effective adiabatic index remains near $\tfrac{4}{3}$. We assumed a light (0.1 GeV) DM particle mass in our study. The DM component is more diffusive and extended. Hence, it brings a less significant impact to the NM density profile near its core. We show how the NM density profile changes with increasing DM mass fraction _DM in the right panel of Figure 2. We observe that when more DM is admixed, the core of the NM component remains almost unchanged. Since the collapse dynamics of a WD are governed by the dense core, where ρ₂ is large enough to initiate electron capture, it is natural to expect generic collapse dynamics for all rigidly rotating DMRWDs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find that the DM component collapses with the NM component to form a bound DM core. We show the DM density profile evolution in the left panel of Figure 3. The DM density evolves similarly to the NM density, but it remains stable after the NM core bounce. We show the DM density profile evolution for a particular model Rigid-0.01 in the right panel of Figure 3 as an example. The DM radius contracts from ∼290 km at $\bar{t}=-0.004\,{\rm{s}}$, to ∼180 km at $\bar{t}=0.01\,{\rm{s}}$. Although the DM radius increases at $\bar{t}=0.02\,{\rm{s}}$, the DM component gradually contracts to ∼200 km at $\bar{t}=0.03\,{\rm{s}}$ and pulsates around ∼180–200 km. This suggests that a bound DM component has formed with negligible mass loss. We show the DM velocity profile evolution of the same DMRWD model in the left panel of Figure 3. The post-bounce velocity shock breaks through the DM surface around $\bar{t}=0.01\,{\rm{s}}$. However, the shock is too weak to unbind the DM component. The shock gradually weakens and becomes a sound wave that propagates inside the DM component. This also explains the pulsation of the DM component between $\bar{t}=0.03$ and 0.049 s.

Zoom In Zoom Out Reset image size

3.1.2. Formation of DM-admixed Neutron Stars

3.1.3. GW Signatures

The nonluminous nature of the DM makes it difficult to be detected through conventional telescopes. The weak electromagnetic signatures from a typical AIC also hinder indirect DM detection by comparing AIC luminosities. Therefore, we rely on the GW signatures generated by both the NM and DM components.

Equation (14) suggests that the moment of inertia tensor I_zz is separable into individual DM and NM components:

$$\begin{eqnarray}\begin{array}{rcl}{I}_{{zz}} & = & {I}_{{zz},1}+{I}_{{zz},2},\\ {I}_{{zz},i} & = & \displaystyle \frac{1}{3}{\displaystyle \int }_{\mathrm{All}\ \mathrm{Space}}{\rho }_{i}{r}^{2}{P}_{2}(\cos \theta )d\tau .\end{array}\end{eqnarray} \tag{ 17 }$$

Since the DM only interacts with NM through gravity, the Euler equation for the DM component does not contain any nontrivial NM-related terms except the gravitational potential Φ. Hence, the GW signature from the AIC of a DMRWD can be separated into the DM and NM contributions:

$$\begin{eqnarray}\begin{array}{rcl}{h}_{+} & = & {h}_{+,1}+{h}_{+,2},\\ {h}_{+,i} & = & \displaystyle \frac{3}{2}\displaystyle \frac{G}{{{Dc}}^{4}}{\sin }^{2}\theta \displaystyle \frac{{d}^{2}}{{{dt}}^{2}}{I}_{{zz},i}.\end{array}\end{eqnarray} \tag{ 18 }$$

To compute h_+,i , we make use of Equation (16) in Ott et al. (2004) and substitute all the components of v and ρ by the corresponding DM/NM values.

It is also a common practice to study GW strains by time-frequency analysis. To obtain the GW spectrogram, we perform a windowed Fourier transform:

$$\begin{eqnarray}&&{\tilde{h}}^{* }(f,t)={\int }_{-\infty }^{\infty }{h}_{+}(\tau )w(t,\tau )\exp (-2\pi {if}\tau )d\tau .\end{eqnarray} \tag{ 19 }$$

Here, w(t, τ) is the window function, and we choose the Hann window.

We first show the AIC GWs generated by the rigidly rotating DMRWD models shown in Figure 4. The GWs are all generic Type I waveforms (Fryer & New 2011). There are no considerable differences in the GW signature with respect to all DM-admixed models. This contrasts with the results presented by Zha et al. (2019b), where they show enhanced amplitudes during $\bar{t}=0\,{\rm{s}}$. This is because the contributions to the GW strains are mainly from the innermost core (∼10 km). We have shown in the previous section that the effects of admixing 0.1 GeV DM on the NM density profile are mainly at the NM outer envelope. The NM collapse dynamics are also generic for all DM-admixed models. We append the NM density contour plots of models NM-Rigid and DM-Rigid-0.2 in Figure 5 for comparison. We observe that the dense core, which corresponds to the major part of the proto-neutron star of the DM-admixed model, is almost identical to that of the pure NM counterpart. This explains why the GW signatures from rigidly rotating DMRWDs are all generic.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, the situation is different for differentially rotating progenitors. We show the GW strains of the Kepler-rotating and α_d = 0.1 model in Figure 6. We find that the DM admixture indirectly suppresses the post-bounce third and fourth peaks of the GW strains. This could also be observed as the gradual disappearance of the third and fourth spectral peaks in Figure 7. Therefore, the GW strains of DMAIC are qualitatively different from that of the pure NM model. We find that the spectral peaks exist for the pure NM model because the reflected shock waves pass through the NM core and make it pulsate non-radially. The corresponding pulsation amplitudes for the DM-admixed models are smaller, resulting in weaker GW signatures. We find similar results for the Kepler-rotating and α_d = 0.01 models, except that the fourth spectral peak never exists for the pure NM, and hence, the DM-admixed models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The DM component is more diffusive for fermionic DM with a particle mass of 0.1 GeV when compared to those with heavier DM particle mass. As such, the collapse dynamics of the DM component could not produce GW amplitudes comparable to that of the NM component. The effects of DM admixture on the total GW signatures are, therefore, indirect. To quantitatively determine whether such effects could be observable, we compute the mismatch ${\mathfrak{M}}$, which quantifies how dissimilar two waveforms are (Reisswig & Pollney 2011; Richers et al. 2017):

$$\begin{eqnarray}&&{\mathfrak{M}}=1-\max \left(\displaystyle \frac{\langle {h}_{a},{h}_{b}\rangle }{\sqrt{\langle {h}_{a},{h}_{a}\rangle \langle {h}_{b},{h}_{b}\rangle }}\right).\end{eqnarray} \tag{ 20 }$$

The second term here contains the match between two waveforms h_a and h_b :

$$\begin{eqnarray}&&\langle {h}_{a},{h}_{b}\rangle ={\int }_{0}^{\infty }\displaystyle \frac{4{\tilde{{h}_{a}}}^{* }\tilde{{h}_{b}}}{s}{df}.\end{eqnarray} \tag{ 21 }$$

Here, s is the estimated noise amplitude spectral density of the Advanced LIGO (Barsotti et al. 2018). ${\tilde{h}}^{* }$ is the Fourier transform of the GW strain, which is just Equation (19) but with w(t, τ) = 1. The mismatch is maximized over the relative phase, amplitudes, and arrival times. We follow Zha et al. (2019b) to set the integration limit of Equation (21) to be from 100–2000 Hz. The computations are facilitated through the open-source package PyCBC (Nitz et al. 2022). We extract GW waveforms for all the models listed in Table 1 with a time window of −0.01 s $\lt \,\bar{t}\lt 0.05\,{\rm{s}}$ and compute the mismatches with respect to the pure NM model. The results are listed in Table 2. The mismatches for the rigidly rotating DMAIC models are small, which is no surprise because the GW waveforms of the DM-admixed models in such a scenario are very similar to that of the pure NM counterpart. The mismatches for the Kepler-rotating DMAIC models, however, are relatively large. The presence of 1% of DM can be inferred from future GW detection produced by DMAIC if Advanced LIGO can distinguish two waveforms with an accuracy better than 14%.

Table 2. GW Mismatch (in %) with Respect to the Pure NM Model for DMAICs with Different Initial Rotation Profiles

⋯	Rigid	Kepler-d001	Kepler-d01
DM-0.01	0.306	13.570	19.217
DM-0.03	0.870	26.151	42.386
DM-0.05	1.263	29.103	41.173
DM-0.07	1.596	36.149	41.169
DM-0.09	1.844	38.256	44.517
DM-0.1	1.992	40.916	43.785
DM-0.2	2.784	41.845	52.794

Note. See Table 1 for the simulation parameters of these models.

Download table as:  ASCIITypeset image

The properties of a Fermionic DM-admixed compact star were shown to be sharply changing around the DM particle mass of 0.1 GeV (Leung et al. 2022). To better capture the transitional effects from a sub-gigaelectronvolt to gigaelectronvolt mass, we include progenitor models admixed with fermionic DM of particle mass 0.3 GeV. Furthermore, the progenitors are all differentially rotating DMRWDs with α_d = 0.5. For reference, we include the parameters of our appended models in Table 3. We generally find similar collapse dynamics for the DM and NM components as those of the diffusive DM limit. For instance, we find a delay in the NM bounce time and the successful formation of a DM-admixed neutron star. An in-depth discussion of the collapse dynamics of DMAIC under the compact DM limit is therefore omitted.

Table 3. Same as Table 1, but for Differentially Rotating DMAIC Progenitors That Have α_d = 0.5, Ω_c = 45.2 s⁻¹, and a DM Particle Mass of 0.3 GeV

Model	MNM	MDM	log10 ρ1c	ReNM	ReDM	DM	tb	log10 ρ1b	log10 ρ2b	MPNS	${\mathfrak{C}}$
⋯	(M⊙)	(M⊙)	(gcm−3)	(km)	(km)	⋯	(ms)	(gcm−3)	(gcm−3)	(M⊙)	(10−3)
Kepler-NM-d05	1.771	⋯	⋯	916	⋯	0.00	29.087	⋯	14.351	1.604	5.71
Kepler-0.01-d05	1.672	0.017	9.968	929	157	0.01	30.595	12.720	14.345	1.516	5.32
Kepler-0.03-d05	1.525	0.047	10.206	948	187	0.03	32.426	12.826	14.339	1.379	4.75
Kepler-0.05-d05	1.410	0.074	10.313	961	202	0.05	33.778	12.856	14.336	1.267	4.33
Kepler-0.07-d05	1.313	0.099	10.383	967	212	0.07	34.899	12.866	14.332	1.170	4.01
Kepler-0.09-d05	1.229	0.122	10.434	967	220	0.09	35.888	12.871	14.332	1.090	3.75
Kepler-0.1-d05	1.191	0.132	10.454	967	223	0.1	36.348	12.873	14.333	1.052	3.64
Kepler-0.2-d05	0.896	0.224	10.588	935	246	0.2	40.396	12.872	14.331	0.767	2.83

Note. We also append the NM compactness ${\mathfrak{C}}=2{{GM}}_{\mathrm{NM}}/{R}_{\mathrm{eNM}}{c}^{2}$.

Download table as:  ASCIITypeset image

3.2.1. GWs from the DM Component

In this section, we focus on the GW signature produced by a DMAIC event. We show the GW strains of a particular model Kepler-0.05-d05 in Figure 8. The DM GW strain has amplitudes comparable to that of the NM and produces secondary oscillations on top of the NM GW strain. This is in contrast to the diffusive DM limit. The DM component, which couples to a highly differentially rotating NM configuration set by α_d = 0.5, becomes spheroidal in shape. The vigorous nonspherical collapse dynamics thus generate a considerable magnitude of GWs. We analyze the GW strains by plotting their spectrograms in Figure 9. We find that the admixture of DM greatly suppresses the NM GW strains around the NM core bounce, which is no surprise because the NM mass and compactness are reduced substantially when DM is admixed (see Table 3). We also find that the DM GW strains show up as a continuous low-frequency (<1000 Hz) signal in the spectrogram before $\bar{t}=0.1\,{\rm{s}}$. This is important because convection-related GW signals are emitted after $\bar{t}=0.1\mbox{--}0.2\,{\rm{s}}$ (Zha 2019). Therefore, any low-frequency signals observed before $\bar{t}=0.1\,{\rm{s}}$ could be direct evidence of a compact DM admixture. The DM GW waveforms (see Figure 10) are consistent with the Type III collapsing polytrope waveforms presented in Fryer & New (2011).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.2. Detection Prospects

To study the detectability of the DM GW signals, we compute the dimensionless characteristic GW strain (Flanagan & Hughes 1998):

$$\begin{eqnarray}&&{h}_{\mathrm{char}}=\sqrt{\displaystyle \frac{2}{{\pi }^{2}}\displaystyle \frac{G}{{c}^{3}}\displaystyle \frac{1}{{D}^{2}}\displaystyle \frac{{{dE}}_{\mathrm{GW}}}{{df}}}.\end{eqnarray} \tag{ 22 }$$

Here, $\tfrac{{{dE}}_{\mathrm{GW}}}{{df}}$ is the GW spectral energy (Murphy et al. 2009):

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{GW}}}{{df}}=\displaystyle \frac{3}{5}\displaystyle \frac{G}{{c}^{5}}{\left(2\pi f\right)}^{2}| {\tilde{h}}_{+}{| }^{2}.\end{eqnarray} \tag{ 23 }$$

We compare h_char f^−1/2 with the Advanced LIGO noise spectral density $\sqrt{s(f)}$ in Figure 11. In the same figure, we mark vertical lines corresponding to the peak frequencies of the DM GW waveforms (see Figure 12). We choose the sampling window as $\bar{t}\gt 0\,{\rm{s}}$. The DM characteristic GW strains corresponding to the frequency peaks are above the Advanced LIGO sensitivity curve, which is true for all of our considered models for all our considered models, assuming D = 10 kpc. Hence, the GW signature of a collapsing, compact DM in a Milky Way DMAIC event should be detectable by Advanced LIGO. Our results represent the first-ever numerical calculation of the GW waveforms of a collapsing DM core in a compact star. Finally, we show the detectability of DM GWs from rigidly rotating progenitors in Appendix B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size