Lost in Space: Companions' Fatal Dance around Massive Dying Stars

To model the change of the orbital parameters of the two-body system (central exploding star and a companion, which is assumed to be either a stellar companion or a planet), we solve the equations of motion of the system numerically in two dimensions. This way, unlike via the analytic methods presented in Veras et al. (2011), we can examine a secondary whose mass is nonnegligible compared to that of the primary. The barycenter is not fixed on the primary, and thus we can analyze the SN II events in a binary star system. The equations of motion for the primary and secondary components are

$$\begin{eqnarray}&&{\ddot{{\boldsymbol{r}}}}_{\mathrm{pri}}^{{\prime} }=-{M}_{\sec }\ G\ \displaystyle \frac{{{\boldsymbol{r}}}_{\mathrm{pri}}^{{\prime} }}{{r}_{\mathrm{pri}}^{{\prime} 3}},\,\,\,\,\,{\ddot{{\boldsymbol{r}}}}_{\sec }^{{\prime} }=-{M}_{\mathrm{in}}\ G\ \displaystyle \frac{{{\boldsymbol{r}}}_{\sec }^{{\prime} }}{{r}_{\sec }^{{\prime} 3}},\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant, and ${{\boldsymbol{r}}}_{\mathrm{pri}}^{{\prime} }$ and ${{\boldsymbol{r}}}_{\sec }^{{\prime} }$ are the position vector of the primary and secondary components in the inertial frame, respectively. M_sec is the secondary's mass, which is a constant value, while the mass inside the secondary's orbit, M_in, changes due to the envelope loss of the primary. To calculate M_in, we develop a homologous envelope loss model; see details in the next section.

Working in Python 3.8.8 and opting for SciPy's integrate.solve_ivp() function with an explicit Runge–Kutta method of order 8 (DOP853), we integrate the primary and the secondary's barycentric position and velocity vectors through the explosion and beyond. The method works with an adaptive step size to reach the desired precision. M_in is updated at every time step during numerical integration.

When a massive star's life ends in a type II supernova explosion, the stellar core collapses and forms a neutron star. Meanwhile, the envelope falls onto this dense core and bounces back. Through this process, a shock wave is initiated, which is then responsible for the ejection of the stellar outer shell. The shock wave propagates through the envelope and reaches the stellar surface. The mass of the ejecta could be as high as 95% of the initial stellar mass (Smartt et al. 2009). Here, we assume that the envelope expands in a spherically symmetric, homologous fashion (e.g., Arnett 1980; Branch & Wheeler 2017), which results in the following: (1) the velocity of each ejected layer is a linear function of distance; (2) the velocity of the outermost layer is independent of time; and (3) the density profile of the ejected material is also time-invariant.

The assumed spherical geometry of the expanding ejecta, at least during the first ∼100 days after explosion, is consistent with spectropolarimetric observations of hydrogen-rich (type II) core-collapse supernovae (e.g., Branch & Wheeler 2017; Nagao et al. 2021). The observed continuum polarization of the majority of such supernovae is found to be close to zero shortly after explosion, but tends to increase as the ejecta expands. This suggests that the inner parts of the ejecta may deviate from spherical symmetry, even though interaction with the surrounding circumstellar matter (CSM) may also contribute to the observed polarization. Motivated by the initial spherical symmetry of the ejecta as observed (and admitting that an asymmetric geometric configuration would introduce additional complexity in our calculations), we assume spherical symmetry of the expanding supernova ejecta in our calculations.

In the following, we summarize the homologous envelope expansion model applied, e.g., in Vinkó et al. (2004). Due to the homologous expansion, the velocity of the envelope at a given distance r(t) from the center is $v(r,t)=\left(r(t)/R(t)\right){v}_{\max }$, where R(t) is the time-dependent radius of the outermost layer of the expanding ejecta and v_max is the maximum expansion velocity at R. Throughout the expansion, v_max is assumed to be constant. Initially the envelope radius R = R₀ is equal to the stellar radius, and is assumed to be R₀ = 500 R_⊙ (∼2.23 au) throughout this study. We assume that the inner 10% of the star contains a constant density core, which has a fractional radius r_c = 0.1R.

Introducing the relative (comoving) distance coordinate as x = r(t)/R(t), the envelope velocity at x is $v(x,t)=x\cdot {v}_{\max }$. Using the above equation and defining Δt as the time elapsed since shock breakout, the distance of a given mass layer from the SN center can be expressed as

$$\begin{eqnarray}&&r(t)=r(0)+v(r,t){\rm{\Delta }}t=x({R}_{0}+{v}_{\max }{\rm{\Delta }}t).\end{eqnarray} \tag{ 2 }$$

When this layer reaches the orbit of the secondary body, $r={r}_{\sec }$,

$$\begin{eqnarray}&&x({r}_{\sec })\equiv {x}_{\sec }=\displaystyle \frac{{r}_{\sec }}{{R}_{0}+{v}_{\max }{\rm{\Delta }}t}.\end{eqnarray} \tag{ 3 }$$

To determine the mass inside r_sec, M_in, three cases should be distinguished: ${x}_{\sec }\gt 1$, or either ${x}_{\sec }\gt {x}_{{\rm{c}}}$ or ${x}_{\sec }\lt {x}_{{\rm{c}}}$. In the first case, all of the stellar mass resides inside the secondary's orbit, and therefore M_in = M_ej + M_n, where M_ej is the mass of the stellar ejecta, and M_n the mass of the remnant neutron star. In the latter cases, we take into account the change of envelope density due to homologous expansion. The density profile of the core is considered to be spatially constant (ρ₀), while the density of the envelope follows a power law,

$$\begin{eqnarray}&&\rho ={\rho }_{0}{\left(\displaystyle \frac{x}{{x}_{{\rm{c}}}}\right)}^{-n},\end{eqnarray} \tag{ 4 }$$

where n = 7 is assumed if x > x_c, and n = 0 otherwise. The core density function changes in time according to

$$\begin{eqnarray}&&{\rho }_{0}(t)=\displaystyle \frac{3{M}_{{\rm{c}}}}{4\pi }{\left({x}_{{\rm{c}}}({R}_{0}+{v}_{\max }{\rm{\Delta }}t)\right)}^{-3}.\end{eqnarray} \tag{ 5 }$$

If ${r}_{\sec }\gt {r}_{{\rm{c}}}$ then the mass residing within the secondary's orbit is

$$\begin{eqnarray}&&{M}_{\mathrm{in}}={M}_{{\rm{n}}}+{M}_{{\rm{c}}}+{\int }_{{r}_{{\rm{c}}}}^{{r}_{\sec }}4\pi {r}^{2}\rho (r){dr},\end{eqnarray} \tag{ 6 }$$

where r_c is the radius of the core. If ${r}_{\sec }\gt {r}_{{\rm{c}}}$, the above equation, after some trivial algebra, gives

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{in}} & = & {M}_{{\rm{n}}}+{M}_{{\rm{c}}}+4\pi {\rho }_{0}{R}_{0}^{3}{x}_{{\rm{c}}}^{n}{\displaystyle \int }_{{x}_{{\rm{c}}}}^{{x}_{\sec }}{x}^{2-n}{dx}\\ & = & {M}_{{\rm{n}}}+{M}_{{\rm{c}}}\left[1+\displaystyle \frac{3}{n-3}\left(1-{\left(\displaystyle \frac{{x}_{\sec }}{{x}_{{\rm{c}}}}\right)}^{3-n}\right)\right].\end{array}\end{eqnarray} \tag{ 7 }$$

That said, if ${x}_{\sec }\lt {x}_{{\rm{c}}}$,

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{in}} & = & {M}_{{\rm{n}}}+4\pi {\rho }_{0}{R}_{0}^{3}{\displaystyle \int }_{0}^{{x}_{\sec }}{x}^{2}{dx}\\ & = & {M}_{{\rm{n}}}+{M}_{{\rm{c}}}{\left(\displaystyle \frac{{x}_{\sec }}{{x}_{{\rm{c}}}}\right)}^{3}.\end{array}\end{eqnarray} \tag{ 8 }$$

With the assumed density profile, the mass of the ejecta, M_ej, can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{ej}} & = & 4\pi {R}_{0}^{3}{\rho }_{0}\left({\displaystyle \int }_{0}^{{x}_{{\rm{c}}}}{x}^{2}{dx}+{x}_{{\rm{c}}}^{n}{\displaystyle \int }_{{x}_{{\rm{c}}}}^{1}{x}^{2-n}{dx}\right)\\ & = & 4\pi {R}_{0}^{3}{\rho }_{0}\left(\displaystyle \frac{{x}_{{\rm{c}}}^{3}}{3}+\displaystyle \frac{{x}_{{\rm{c}}}^{n}-{x}_{{\rm{c}}}^{3}}{3-n}\right).\end{array}\end{eqnarray} \tag{ 9 }$$

From this expression, one can derive the core density as

$$\begin{eqnarray}&&{\rho }_{0}=\displaystyle \frac{{M}_{\mathrm{ej}}}{4\pi {R}_{0}^{3}}{\left(\displaystyle \frac{{x}_{{\rm{c}}}^{3}}{3}+\displaystyle \frac{{x}_{{\rm{c}}}^{n}-{x}_{{\rm{c}}}^{3}}{3-n}\right)}^{-1}.\end{eqnarray} \tag{ 10 }$$

Now the core mass can be calculated as

$$\begin{eqnarray}&&{M}_{{\rm{c}}}=\displaystyle \frac{4\pi {R}_{0}^{3}{x}_{{\rm{c}}}^{3}}{3}{\rho }_{0}=\displaystyle \frac{{M}_{\mathrm{ej}}}{1+\tfrac{3}{n-3}\left(1-{x}_{{\rm{c}}}^{n-3}\right)}.\end{eqnarray} \tag{ 11 }$$

The upper panel of Figure 1 shows an example of the evolution of the radial density profile of the stellar core and envelope at three different phases. The lower panel of Figure 1 shows some examples of the time evolution of the total mass residing inside the secondary for three different expansion velocities.

Zoom In Zoom Out Reset image size

The total integration time for each specific model was 5 million days (about 13,700 yr). On average, by the end of the simulations, the measured mass-loss rate is 10⁻¹¹–10⁻¹⁸ M_⊙ day⁻¹, meaning that the mass loss saturates; see the lower panel of Figure 1. This allows the monitored values (semimajor axis, a, eccentricity, e, the velocity of the primary, v_pri, and the secondary, v_sec) to relax by the end of the simulation. Using relative and absolute tolerance values of 10⁻¹² for a fixed-mass two-body problem with an assumption of a high-mass secondary (10 M_⊙), the relative error in total energy is on the order of 10⁻⁹ by the end of integration time.

We monitor the perturbation and evolution of the secondary's eccentricity and semimajor axis throughout simulations using the following equations:

$$\begin{eqnarray}&&e=\sqrt{1+2h{\left(\displaystyle \frac{{{\boldsymbol{r}}}_{\sec }\times {{\boldsymbol{v}}}_{\sec }}{\mu }\right)}^{2}},a=-\displaystyle \frac{\mu }{2h},\end{eqnarray} \tag{ 12 }$$

where ${{\boldsymbol{r}}}_{\sec }$ and ${{\boldsymbol{v}}}_{\sec }$ are the secondary's heliocentric position and velocity vectors, respectively. The heliocentric position and velocity vectors of the primary are derived from the barycentric components. In Equations (12), the energy term reads

$$\begin{eqnarray}&&h=\displaystyle \frac{1}{2}{v}_{\sec }^{2}-\displaystyle \frac{\mu }{{r}_{\sec }},\end{eqnarray} \tag{ 13 }$$

and $\mu =G({M}_{\sec }+{M}_{\mathrm{in}})$. Here, r_sec and v_sec are the secondary's heliocentric distance and speed, respectively.

To map a wide range of scenarios, we investigate different initial conditions. Model set A maps various orbital elements (initial semimajor axis a_init, initial eccentricity e_init, and initial true anomaly ν) and secondary masses (M_sec). In this set, we also test different envelope and remnant masses. Neutron star mass is constrained by the Chandrasekhar limit, and an upper limit of 3 M_⊙ is adopted based on theory and observational evidence (Kalogera & Baym 1996; Clark et al. 2002). As for the velocity of the outermost envelope layer, i.e., the ejection velocity, we adopt the 1000–10,000 km s⁻¹ range based on Hamuy & Pinto (2002). In the investigated semimajor axis range this results in the ejecta reaching the secondary well before it completes its first orbit after the explosion. The parameters used in Model set A can be found in Table 1. This gives us a total of 1,875,000 models to investigate. Only those models are analyzed where the secondary mass is below that of the primary.

Table 1. Parameters Used in Model Set A

ainit	$[{R}_{0}-100\ \mathrm{au}]$ a
einit	0.0, 0.1, 0.2, 0.4, 0.8
ν	0.0, π/4, π/2, 3π/4, π
vmax	[1000–10,000] km s−1 a
Mej	6.5 M⊙, 10 M⊙, 15 M⊙, 22 M⊙
Mn	1.5 M⊙, 2 M⊙, 3 M⊙
Msec	1 M⊕, 10 M⊕, 1 MJ, 10 MJ,
	1 M⊙, 2 M⊙, 3 M⊙, 5 M⊙, 10 M⊙, 20 M⊙

Notes. M_⊕ denotes Earth, M_J Jupiter, and M_⊙ solar mass.

^a Twenty-five values are sampled linearly from the given ranges.

Download table as:  ASCIITypeset image

To map the effect of true anomaly on the evolution of the system, we run Model set B, an additional 720,500 models. The settings of these models are M_ej = 6.5 M_⊙ and M_n = 3 M_⊙. A total of 1441 different true anomaly values are sampled linearly from the [0–2π] range, and 25 a_init values from the range [R₀–100 au]. We assumed that the secondary orbits outside the envelope in each model, as whether planets can survive in stellar envelopes still remains a question (see, e.g., Setiawan et al. 2010, Bear et al. 2011, and more recently Chamandy et al. 2021, Lagos et al. 2021, and Szölgyén et al. 2022). We investigate both planetary (${M}_{\sec }=1{M}_{\oplus }$) and stellar (${M}_{\sec }=1,3,5\,\,\mathrm{and}\,7{M}_{\odot }$) companions, whose eccentricity is 0.4 or 0.8. The velocity of the envelope is set to be 1000 km s⁻¹ or 10,000 km s⁻¹.

A given model is considered to be bound if the secondary's eccentricity value remains under unity, and unbound if it grows beyond unity by the end of the simulation. Models with a planetary and a stellar-mass companion show qualitatively different outcomes: bound planets, (unbound) free-floating planets, bound stars, and free-floating stars.

In this section, we investigate cases in which a bound system is formed. In the center, there is always a neutron star, while the companion can be either a planet or a star, depending on the initial condition. Figure 2 shows the orbital parameters (a and e) of the companion planet or star for different initial conditions, calculated according to Equation (12) at the end of the integration when these orbital parameters are relaxed.

Zoom In Zoom Out Reset image size

With regards to the incidence of a bound state, we find that by increasing the mass of the neutron star, the chance of forming a bound orbit also increases, as shown in panels a1–b3 of Figure 2. This finding is in agreement with Veras et al. (2011) for planetary companions and is also found to be valid for stellar companions according to our simulations. The formation of a bound state requires a planet initially close to the apocenter (see details in Section 3.3.). We also find that the semimajor axis of bound planets increases with increasing the mass of the ejecta, which can be seen in panels a1–a3 of Figure 2.

By analyzing semimajor axis histograms normalized by a_init, a clear trend can be identified in the planetary case. The ratio of the final-to-initial semimajor axis values shows discrete peaks (due to the assumption of discrete ejecta mass values), and the value of the peak centers increases with ejecta mass. Interestingly, there is also a minimum value of about 2 for this ratio. By contrast, for stellar-mass companions, such a trend cannot be identified.

The final eccentricities of bound planets, shown in panels c1–c3, are found to be clumpy and the number of clumps increases with increasing M_n. This is due to the fact that we use discrete M_ej and M_n values, and the number of discrete M_ej/M_n bound models increases with M_n. There is an eccentricity minimum (being the most probable value), for all pairings, that also decreases with M_n and increases with M_ej. Final a and e values are independent of planetary mass.

We emphasize that the initial eccentricity of a planetary companion should be at least 0.4 to form a bound system assuming plausible SN II explosion scenarios. By analyzing the results we reveal that the orbit of a planetary companion is circularized in all cases. An exception to the above statement occurs for the most massive neutron star models (M_n = 3 M_⊙) with the least massive ejecta mass models (M_ej = 6.5 M_⊙, shown in panel c3). In this case, a bimodal peak can be found at e = 0.4 and e = 0.9, corresponding to e_init = 0.8 and e_init = 0.4, respectively, with the latter case meaning that the eccentricity of the planet increases.

Assuming a stellar-mass companion, the clumpy structure in eccentricity is absent, as shown in panels d1–d3 of Figure 2. We find that a bound system can form assuming arbitrary initial eccentricity (even a circular orbit with e_init = 0). As a consequence both circularization and the excitation of eccentricity occur. Stellar systems do not require a specific range of true anomalies to stay bound. Contrary to the planetary-mass companion models, e increases with increasing M_sec. A wider final stable orbit also coincides with a more eccentric one.

For stable binary systems, we observe that the system gains non-zero peculiar velocity, even though the explosion was assumed to be completely symmetric. For a planetary-mass secondary companion, this peculiar velocity is found to be on the order of 10⁻⁶ km s⁻¹, which is observationally insignificant. To find the largest peculiar velocity binary, one must consider the smallest possible companion, highest ejecta mass, and eccentricity, and the secondary should be at pericenter, where the velocity of the components is the largest. For stellar binaries, we find significantly higher velocities, up to 37 km s⁻¹, which can presumably be detected.

Investigation of the effect of ejecta velocity showed that higher v_max always results in a lower e value for the bound system. This finding is independent of the companions' mass, and thus is true for both planetary and stellar-mass companion models.

For unbound systems, the final eccentricity of the secondary grows above unity, corresponding to a parabolic or hyperbolic unbound orbit. In these scenarios, we analyze the distribution of the final velocities of both the primary and secondary components. Figure 3 shows v_pri and v_sec for an unbound state both for planetary and stellar companions. The histograms of e, v_pri, and v_sec are also given in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For planetary companion models, v_pri are grouped by M_sec and the velocity values increase with M_sec; see panel a of Figure 3. However, v_sec is independent of the mass of the planet. Concerning the true anomaly, the sensitive parameter however is only the secondary's velocity. One can see in panel b of Figure 3 that the closer v is to π (i.e., apocenter), the lower the secondary velocity, while the primary velocity can span a wide range of values.

For stellar-mass companion models (panels c and d of Figure 3) v_pri values are three orders of magnitude higher than that of planetary-mass companion models, while v_sec values display the same range of values for both models. Notably, low primary and high secondary velocity systems are absent, as are systems with high primary and low secondary velocities. Both velocities are sensitive to stellar mass and true anomaly, as shown in panels c and d of Figure 3. Based on the clear separation of colors representing the secondary's mass in panel c, a higher M_sec results in a higher v_pri and slightly lower v_sec. Similarly to planetary-mass companion models, low v_sec develops if the secondary is close to apocenter.

An interesting feature of the secondary eccentricities shown on panels a1 and a2 of Figure 4 is that there exists a most likely value for planetary systems around e = 4. This trend is absent for stellar companions, where the larger the eccentricity the lower its occurrence.

Two populations of primary velocity (with narrow and wide distributions) can be distinguished; see panel b1 of Figure 4. The width of the distribution of v_pri increases with M_sec. Very low primary velocities (∼10⁻⁶−10⁻⁴ km s⁻¹) correspond to ${M}_{\sec }=1\,{M}_{\oplus }$ and ${M}_{\sec }=10\,{M}_{\oplus }$, while an order of magnitude larger primary velocities (up to 0.1 km s⁻¹) correspond to ${M}_{\sec }=1\,{M}_{{\rm{J}}}$ and ${M}_{\sec }=10\,{M}_{{\rm{J}}}$. The primary velocities for stellar companion models, shown in panel b2 of Figure 4, also exhibit distinct populations based on M_sec, where a higher M_sec also yields a higher v_pri value. The most likely value for primary velocity in the stellar case is ∼4 km s⁻¹.

The secondary's velocity histograms of planetary and stellar companion models (panels c1 and c2 in Figure 4) are remarkably similar: the most likely secondary velocity is ∼18 km s⁻¹, while the highest measured value is ∼270 km s⁻¹. A notable feature of both primary and secondary velocities is that few models can result in velocities higher than v_pri = 30 km s⁻¹ and ${v}_{\sec }\,$= 100 km s⁻¹ for a stellar companion. For a planetary-mass companion, v_pri is smaller by three orders of magnitude.

By analyzing the primary and secondary velocities as a function of v_max, we found no clear trend. This means that neither the primary nor the secondary velocities are sensitive to the speed of the exploding envelope in the investigated parameter space.

In accordance with Veras et al. (2011), we also found that in order to have a bound system, an eccentric orbit is required for planetary-mass companions. For stellar companion models, however, a companion on a circular orbit initially can also remain bound to the neutron star. For an eccentric orbit, only a certain range of initial true anomaly (true anomaly at the instance of explosion) values, ν_b, results in bound systems.

We thoroughly investigated the effect of true anomaly on planetary and binary systems with an initial eccentricity of either 0.4 or 0.8 and an ejecta speed of either 1000 or 10,000 km s⁻¹. For the planet, we assumed 1 M_⊕, though we find that the results can be generalized to the whole examined planetary mass range. Results are shown in Figure 5 for three a_init. Four panels are shown for each model: a, e , v_pri, and v_sec as a function of ν.

Zoom In Zoom Out Reset image size

The higher the semimajor axis, the fewer systems stay bound in all model sets, as visible on panels a1, b1, c1, and d1 in Figure 5. It is clearly seen that ν_b shifts toward higher initial true anomaly values if a larger a_init is used for low v_max models. However, this trend diminishes for higher ejecta velocities; therefore, ν_b becomes symmetric around the apocenter. As one can see in Figure 5, for ${v}_{\max }=1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$, a_init affects e, (panels a2 and b2); however, e becomes independent of a_init for ${v}_{\max }=10,000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (see panels c2 and d2). The width of the bound region is found to be proportional to e_init independent of v_max (meaning that more bound systems can form with higher e_init). If ${v}_{\max }=1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$, ν_b ∈ [165°–235°], while for ${v}_{\max }=10,000\,\mathrm{km}\,{{\rm{s}}}^{-1}$, ν_b ∈ [155°–200°]. In general, a lower a_init results in lower eccentricities, clearly visible for unbound orbits. However, with increasing v_max, the eccentricity becomes independent of a_init.

The primary and secondary velocities measured in unbound models show a strong dependence on a_init and ν. In general, the component velocities are found to decrease with a_init. Interestingly, v_sec curves are highly symmetric around the apocenter for high v_max values, for which symmetry starts to break for a lower ejecta velocity. By contrast, v_pri curves are never symmetric. Note that, with e_init = 0.8, the velocities of components are almost doubled compared to e_init = 0.8, while ${v}_{\max }=1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ results in less than 50% variation in the secondary, and less than 10% variation in primary velocity. The above findings correspond quite well to the analytically derived true anomaly ranges in Veras et al. (2011). The asymmetry seen in Figure 5 can be explained by the difference in the definition of ν; see details in Section 4.

Figure 6 is similar to Figure 5, the only difference being that we modeled a binary system in which the mass of the secondary is set to be 5 M_⊙. One can see that ν_b spans a wider range of values, [110°–340°] (panel a3), but can be as low as [140°–225°] (panel d3). Moreover, it is clearly seen in the velocity panels (a3, a4, b3, b4, c3, c4, d3, and d4) that almost 50% of models remain bound. Concerning the symmetry and asymmetry of velocities, we found that the secondary velocities are highly asymmetric assuming ${v}_{\max }=1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$, contrary to the planetary case. However, for ${v}_{\max }=10,000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ the symmetry is somewhat restored.

Zoom In Zoom Out Reset image size

Another important behavior of binary systems is that the number of bound systems is about two times higher for e_init = 0.4 than for e_init = 0.8, contrary to the planetary-mass companion models. This simply means that a binary system can remain bound if e_init is low. Finally, we found that all other conclusions drawn from the planetary companion models are also valid here.

Figure 7 shows the evolution of the secondary's eccentricity as a function of time for all possible outcomes of models for both planetary (panel a) and stellar companions (panel b). In agreement with Veras et al. (2011), we found that e_init ≥ 0.4 is required for the formation of a bound planetary system. We also confirm that at early phases the eccentricity decreases, then attains a minimum value (which never reaches zero) and starts to increase for systems with e_init ≥ 0.8.

Zoom In Zoom Out Reset image size

With regards to the stable orbits' dependence on true anomaly, ν_b, Veras et al. (2011) found that ν_b is symmetric around the apocenter. However, according to our findings ν_b is asymmetric and is pushed toward higher ν values. This contradiction can be resolved by the fact that the true anomaly is measured at the beginning of the explosion in our models. If ν were measured at the instance when the outer envelope layer reaches the planet, ν_b would be symmetric. This means that asymmetry is simply caused by orbital movement of the secondary during the early phase of envelope expansion.

Our numerical approach made it possible to examine cases where the secondary mass is nonnegligible relative to that of the primary. With regard to these stellar-mass companion models (panel b of Figure 7), we got novel results. Although bound systems show similar behavior to planetary models, decreasing eccentricity is not observed for free-floating stars. Another major difference compared to planetary systems is that there is no minimum requirement of e_init to form a bound orbit.

In Veras et al. (2011) a simple linear mass-loss model was assumed, in which case $\dot{{M}_{\mathrm{in}}}$ is constant. In our homologous expansion model described in Section 2, we apply a more elaborate version of an explosion model, in which the magnitude of $\dot{{M}_{\mathrm{in}}}$ is not time-invariant (see Figure 1). In this case $\dot{{M}_{\mathrm{in}}}$ slowly increases then gains a constant value, which slowly saturates to zero at later stages. To compare the two explosion models, we conduct a set of representative planetary and stellar-mass companion models. Three linear models are compared with different mass-loss rates in each case. Figure 8 shows two planetary-mass companion models assuming ${v}_{\max }=5000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (panel a) and ${v}_{\max }=1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (panel b). In both runs, we found that the linear models can give qualitatively and quantitatively different results compared to the homologous model. If the constant mass-loss rate is slightly smaller than in the homologous model, then the semimajor axis of the bound state can be smaller. However, if the mass-loss rate underestimates that of the homologous expansion, the system can become unbound. The linear model departs from the homologous model if it assumes the constant mass-loss rate value of the homologous model. The observed differences hold regardless of the ejecta speed. We found an easy method to get a somewhat similar solution for constant and homologous mass-loss models: appropriate $\dot{{M}_{\mathrm{in}}}$ can be derived by fitting a straight line connecting the primary mass and ∼0.6 M_ej points.

Zoom In Zoom Out Reset image size

For stellar-mass companion models, the homologous and linear models give qualitatively the same results. Variation in $\dot{{M}_{\mathrm{in}}}$ of the linear model, however, can give slightly different orbital parameters in bound states, or different velocity vectors in unbound states.

Although we use an elaborate mass-loss approximation, the homologous expansion model intrinsically assumes spherical symmetry. This simplification cannot explain most of the observations (e.g., planet-hosting PSR B1257+12's peculiar velocity of 326 km s⁻¹ (see Yan et al. 2013), or the velocities presented in Hobbs et al. 2005). A pulsar with a high peculiar velocity (e.g., 1100 km s⁻¹ in the case of PSR B1508+55; Gvaramadze et al. 2008) can be indirect proof of a strongly asymmetric explosion. Asymmetric envelope ejection complicates the perturbation of the orbital elements of the secondary, as well as the change in orbital elements of the primary (see, e.g., Parriott & Alcock 1998; Namouni 2005; Namouni & Zhou 2006). In this case, our prediction for v_pri and v_sec are presumably inadequate.

With regard to the envelope size, R₀, we assumed it to be independent of primary mass, which predicts certain M_ej−M_n pairings, which might affect the results.

By the nature of a one-dimensional approximation, the secondary-envelope interactions and the subsequent orbital perturbations are also neglected in this study. For example, local perturbation in the envelope mass distribution by the secondary can break the assumed spherical symmetry. As a result, both the secondary and primary orbital elements are subject to change. Moreover, in a three-dimensional model, an additional orbital element, i.e., inclination, should be taken into account, which further complicates the picture, though the numerical N-body simulations for this limited number of bodies is not computationally demanding.

In a more sophisticated model, the secondary-envelope interactions (see, e.g., Debes & Sigurdsson 2002) could result in mass loss of the secondary, which effect was completely neglected in our orbital solutions. To investigate this phenomenon, a mass-loss and accretion model for the secondary is required. Presumably, this effect is more important for stellar binary configurations.

Handling the abovementioned limitations may affect the orbital parameters of bound systems and velocity components of unbound systems, or can even give qualitatively different results. Future investigations will also address subsequent supernova explosions in binary systems and hierarchical triple system configurations (circumbinary—continuing the work of Fagginger Auer & Portegies Zwart 2022, circumprimary, circumsecondary planets, and star-giant planet-moon systems). These hierarchical systems can be investigated with our homologous expansion model. Adding another planet (Debes & Sigurdsson 2002) and adding a Kuiper belt object (Bonsor et al. 2011) have been investigated in detail for white dwarfs and a white dwarf with a main-sequence companion. We emphasize that both studies use extremely low mass-loss rates (10⁻⁸−10⁻⁶ M_⊙ yr⁻¹, 2/3 × 10⁻⁵ M_⊙ yr⁻¹), which resemble asymptotic giant branch mass loss.

Our results imply that the formation of close-in (a_init < 2 au) pulsar planets requires additional physics, as bound systems have, in general, a much wider separation. Note that PSR B0329+54 b (Starovoit & Rodin 2017) can be explained by our simple model if the system was initially in a high eccentricity state (e_init ≥ 0.4) and the planet was very close to the stellar surface. Since the final-to-initial semimajor axis ratio is found (in the investigated parameter space) to be always at least 2, pulsar planets should have orbited the progenitor inside half of their current semimajor axis. We emphasize that a bound system also requires fortunate explosion timing, i.e., the planet should be close to its apocenter.

With regards to the first exoplanetary system discovered around a pulsar, PSR B1257+12 (Wolszczan & Frail 1992; Wolszczan 1994, 2012), we think that the probability of capturing a single planet is low, and multiple captures are highly improbable. The missing physics to explain the existence of such systems could be a second-generation planet formation hypothesis, in which the pulsar planets form after the explosion from the fallback material. Although disks formed this way would be highly irradiated by the PSR (Martin et al. 2016), they are found to be similar in mass to protoplanetary disks around low-mass young stars (Lin et al. 1991; Menou et al. 2001; Williams & Cieza 2011). Planet formation scenarios in PSR disks have been studied and summarized in Phinney & Hansen (1993), Thorsett et al. (1993), and Podsiadlowski (1993).

An interesting consequence of our results is that rogue planets most probably originate from a low eccentricity orbit. These planets have an average velocity of 1–275 km s⁻¹. In these models, the newborn pulsar gains a peculiar velocity regardless of the symmetric nature of the explosion. Therefore, low peculiar velocity (up to 0.1 km s⁻¹ for Jupiter-sized, and 10⁻³ km s⁻¹ for Earth-sized planets) pulsars most probably originate from pulsar-planet systems.

Bound stellar systems can later become neutron star–white dwarf pairs or binary neutron stars depending on secondary mass. Simulations in this study reveal that the eccentricity and semimajor axes of these systems show a homogeneous distribution on a wide range. Postexplosion binary systems will gain a peculiar speed of up to 37 km s⁻¹ in the investigated parameter space. This can be explained by the fact that the exploding shell inherits the velocity of the primary at the moment of explosion, while the binary gains equal but opposite momentum due to the conservation of momentum. This can be a novel channel for the formation of young expanding shells missing the central pulsar, besides the obvious explanations of the nondetection of the pulsar signal or an asymmetric explosion. Another possible outcome can be the production of peculiar speed (a few tens of km s⁻¹) binary stellar corpses after a second explosion. The formation of binary neutron star systems is going to be the subject of a future study, where a subsequent explosion of the secondary will be modeled. Based on our results we predict that the binary neutron star velocity can be two times as high as observed in the single-explosion models, which is still below 100 km s⁻¹.

Disruption of binary systems can produce two groups of stellar remnants. In the first group, there is a pulsar and a white dwarf, for which case the progenitor of the white dwarf is under the mass limit of an SN II explosion: it is ≤8 M_⊙. In this case, our results show that the most likely primary pulsar velocities are 1–25 km s⁻¹, while white dwarf velocities are 1–128 km s⁻¹.

In the second group, two neutron stars are born, and the primary neutron star velocity is most likely in the range of 5–40 km s⁻¹, while the velocity maximum of the secondary is decreased to 100 km s⁻¹. The increased velocity of the primary neutron star is due to the fact that the secondary should have a mass above 8 M_⊙.

We found an upper limit of neutron star velocity (275 km s⁻¹) based on our secondary velocity measurements. This prediction implicitly assumes that the secondary's subsequent SN II explosion is fully symmetric, so does not change the velocity vector of the neutron star. Contradicting our results, PSR B1508+55 has a 1100 km s⁻¹ peculiar velocity (Gvaramadze et al. 2008); however, the assumption of a strongly asymmetric explosion might resolve this inconsistency.