Impact of Pulsar and Fallback Sources on Multifrequency Kilonova Models

Especially in newly formed NSs like the merged compact object in NS/NS binaries, pulsar emission can vary due to magnetic field restructuring, high-order field configurations, etc. Here we assume a simple dipole magnetic field source with a constant magnetic field following the pulsar luminosity formulae explored in recent studies (Lasky & Glampedakis 2016; Li et al. 2018; Piro et al. 2019). Following Lasky & Glampedakis (2016), the loss of NS remnant angular kinetic energy is balanced by electromagnetic (EM) dipole and GW quadrupole luminosity,

$$\begin{eqnarray}&&-I{\rm{\Omega }}\dot{{\rm{\Omega }}}=\displaystyle \frac{{B}_{p}^{2}{R}^{6}{{\rm{\Omega }}}^{4}}{6{c}^{3}}+\displaystyle \frac{32{{GI}}^{2}{\epsilon }^{2}{{\rm{\Omega }}}^{6}}{5{c}^{2}},\end{eqnarray} \tag{ 2 }$$

where the variables are defined by Lasky & Glampedakis (2016). The pulsar luminosity is given as the EM dipole luminosity multiplied by an efficiency parameter, η, corresponding to the radiation beaming angle (Rowlinson et al. 2014; Lasky & Glampedakis 2016). Assuming GW dominated spin-down, the luminosity is (Lasky & Glampedakis 2016)

$$\begin{eqnarray}&&{L}_{d}(t)=\displaystyle \frac{\eta {B}_{p}^{2}{R}^{6}{\rm{\Omega }}{\left(t\right)}^{4}}{6{c}^{3}}={L}_{0}{\left(1+\displaystyle \frac{t}{{t}_{\mathrm{gw}}}\right)}^{-1},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{L}_{0}=\displaystyle \frac{\eta {{\rm{\Omega }}}_{0}^{4}{B}_{p}^{2}{R}^{6}}{6{c}^{3}},\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{t}_{\mathrm{gw}}=\displaystyle \frac{5{c}^{5}}{128{GI}{\epsilon }^{2}{{\rm{\Omega }}}_{0}^{4}}.\end{eqnarray} \tag{ 4b }$$

To explore the effect of variable remnant lifetime, we introduce a time cutoff, t_cut, to the luminosity,

$$\begin{eqnarray}&&{L}_{d}(t)={L}_{0}{\left(1+\displaystyle \frac{t}{{t}_{\mathrm{gw}}}\right)}^{-1}{\rm{\Theta }}({t}_{\mathrm{cut}}-t),\end{eqnarray} \tag{ 5 }$$

where Θ is the unit step function.

The energy source E_source in Equation (1) for our pulsar models can then be written as

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{source}}=-\displaystyle \frac{{E}_{\mathrm{source}}}{t}+{L}_{d},\end{eqnarray} \tag{ 6 }$$

where L_d is given from Equation (5) and E_d(0) = 0. Note that Equation (6) is a simplification of the ejecta layer equations given by Metzger (2017), neglecting diffusion and changes to the inner velocity from the impulse of the pulsar luminosity. Thus, we have assumed that the pulsar luminosity is not high enough to significantly impact the morphology of the ejecta. The solution to Equation (6) is

$$\begin{eqnarray}&&{E}_{\mathrm{source}}(t)={L}_{0}\displaystyle \frac{{t}_{\mathrm{gw}}}{t}\left[t^{\prime} -{t}_{\mathrm{gw}}\mathrm{ln}\left(\displaystyle \frac{t^{\prime} +{t}_{\mathrm{gw}}}{{t}_{\mathrm{gw}}}\right)\right]\,,\end{eqnarray} \tag{ 7 }$$

where $t^{\prime} =\min (t,{t}_{\mathrm{cut}});$ the derivative with respect to $t^{\prime} /{t}_{\mathrm{gw}}$ of the term in braces is always positive, showing E_d(t) ≥ 0.

In our simulations, E_d(t₀) is added to the innermost spatial cell at a start time of t₀. If t_cut > t₀, the pulsar is still active during the simulation, and Equation (5) is integrated over each time step to add further energy to the innermost cell. In all simulations, we use ${t}_{0}={10}^{4}\,{\rm{s}}$.

Another source of energy can come from fallback after the initial kilonova explosion. Fallback in stellar outbursts follows a simple power law with time ($\dot{m}\propto {t}^{-5/3}$ where $\dot{m}$ is the fallback accretion rate and t is the time; Chevalier 1989). The energy and mass ejected is more complex and high-resolution models have been studied in the case of supernova fallback (Fryer 2009). These models showed that ∼10%–25% of the fallback matter is re-ejected, carrying away roughly 10%–25% of the accretion energy. These properties are directly applicable to the fallback in this scenario and are similar to the simple prescriptions used in the kilonova community, e.g., Li et al. (2018)

$$\begin{eqnarray}&&{L}_{f}(t)=\displaystyle \frac{\eta }{{\eta }_{r}}{L}_{0}\left(\displaystyle \frac{{\dot{m}}_{0}}{{\dot{m}}_{r}}\right){\left(\displaystyle \frac{t}{{t}_{\mathrm{acc}}}\right)}^{-5/3}={\tilde{L}}_{0}{\left(\displaystyle \frac{t}{{t}_{\mathrm{acc}}}\right)}^{-5/3},\end{eqnarray} \tag{ 8 }$$

where η is again an efficiency parameter (roughly 10%–25%), t_acc is initial time (when the luminosity begins to fall off), ${\dot{m}}_{0}$ is the initial accretion rate, and $\dot{{m}_{r}}$ is a reference accretion rate. Reference values η_r and ${\dot{m}}_{r}$ are taken to be 0.1 and 10⁻³ M_⊙ s⁻¹, respectively.

As with our pulsar emission, we do not modify the ejecta mass or momentum in our outflow, focusing instead on the energy injected from this accretion. This luminosity corresponds to a source energy of

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{source}}=-\displaystyle \frac{{E}_{\mathrm{source}}}{t}+{L}_{\mathrm{fallback}},\end{eqnarray} \tag{ 9 }$$

and L_fallback is given by Equation (8). Assuming the luminosity is proportional to t^−5/3 back to t_acc, then E_source can be solved analytically from Equation (9). We assume the luminosity t^−5/3 dependence holds as early as t_acc (which may only be true at t ≫ t_acc Metzger 2017). Assuming the fallback luminosity is constant before t_acc, the solution is

$$\begin{eqnarray}&&{E}_{\mathrm{source}}=\left\{\begin{array}{l}{\tilde{L}}_{0}\tfrac{{t}_{\mathrm{acc}}^{2}}{2t}+{\tilde{L}}_{0}{t}_{\mathrm{acc}}^{5/3}\left(\tfrac{1}{{t}^{2/3}}-\tfrac{{t}_{\mathrm{acc}}^{1/3}}{t}\right)\,,\,\,t\gt {t}_{\mathrm{acc}}\,,\\ \tfrac{{t}_{\mathrm{acc}}{\tilde{L}}_{0}}{2}\,,\,\,t\leqslant {t}_{\mathrm{acc}}.\end{array}\right.\end{eqnarray} \tag{ 10 }$$

Unlike the pulsar models, we do not introduce a cutoff time for the source.

For the fallback models, a considerable amount of energy can be injected into the ejecta on short timescales. Adopting parameters similar to Li et al. (2018), η = 0.1, ${L}_{0}=2\times {10}^{51}$ erg s⁻¹, and ${\dot{m}}_{0}={10}^{-3}$ M_⊙ s⁻¹, from Equation (10), the energy added up to t_acc is 10⁵⁰ erg. The kinetic energy of our smallest-mass ejecta is $\sim {M}_{\mathrm{ej}}{\left({v}_{\max }/2\right)}^{2}/2\approx {10}^{50}$ erg. Assuming all the energy before t_acc goes into boosting the kinetic energy, the resulting average velocity should be increased to about ${\tilde{v}}_{\max }/2=\sqrt{2}{v}_{\max }/2\,\approx 0.42c$. Consequently, we also simulate a suite of fallback kilonova models with ${v}_{\max }/2=0.45$ and with

$$\begin{eqnarray}&&{E}_{\mathrm{source}}=\left\{\begin{array}{c}{\tilde{L}}_{0}{t}_{\mathrm{acc}}^{5/3}\left(\tfrac{1}{{t}^{2/3}}-\tfrac{{t}_{\mathrm{acc}}^{1/3}}{t}\right),\,\,t\gt {t}_{\mathrm{acc}},\\ 0,\,\,t\leqslant {t}_{\mathrm{acc}},\end{array}\right.\end{eqnarray} \tag{ 11 }$$

as an alternative to adding the energy for t < t_acc as a radiative source. We must note that this adjustment to the model does not fully account for the morphological effects of the early source, which would squeeze the ejecta to produce a morphology more similar to those of Metzger (2017) and Li et al. (2018).

We test the effect of varying the pulsar lifetime in 1D spherical models of wind-like outflow, assuming Fe for the wind opacity. The model parameters are summarized in Table 1.

The parameters for the pulsar imply t_gw = 495 s and ${L}_{0}=3.33\times {10}^{44}$ erg s⁻¹, similar to Li et al. (2018).

Figure 2 shows the bolometric luminosity versus time for these cutoff time variations (see model parameters above). The maximum peak luminosity is achieved for ${t}_{\mathrm{cut}}\gtrsim 2\times {10}^{5}\,{\rm{s}};$ increasing the remnant lifetime further only affects the brightness of the tail of the light curve. For ${t}_{\mathrm{cut}}\lesssim 1\times {10}^{5}\,{\rm{s}}$ (∼day), the peak luminosity is more sensitive to the remnant lifetime. This sensitivity can be seen in the change of the 1-day luminosity with respect to the t_cut.

Zoom In Zoom Out Reset image size

For sufficiently low cutoff times, the effect of the pulsar luminosity becomes small relative to the r-process heating. Notably, in Figure 2 the light curve for the lowest cutoff time appears to be monotonically decreasing. This is an effect of the fast expansion speed, low mass, low opacity, and choice of morphology. In particular, for this morphology the time of peak bolometric luminosity follows (Wollaeger et al. 2018)

$$\begin{eqnarray}&&\begin{array}{l}{t}_{\mathrm{peak}}\\ \approx \left(1\ \mathrm{day}\right){\left(\displaystyle \frac{\kappa }{10{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{0.35}{\left(\displaystyle \frac{{M}_{\mathrm{ej}}}{{10}^{-2}{M}_{\odot }}\right)}^{0.318}{\left(\displaystyle \frac{{v}_{\max }}{0.2c}\right)}^{-0.6}.\end{array}\end{eqnarray} \tag{ 13 }$$

Equation (13) gives t_peak ≈ 0.05 day, assuming κ = 0.1 cm² g⁻¹. Alternatively, using the scaling relation for time of peak bolometric luminosity from Grossman et al. (2014), t_peak ≈ 0.09 day. The earlier peak time relative to the models of Grossman et al. (2014) is due to the thermal energy contribution to the light curve (Wollaeger et al. 2018).

With a better understanding of the role of the cutoff time, we can now study the dependence of the pulsar-powered light curves on the magnetic field strength for a range of ejecta masses and two different compositions using just two cutoff times: 2 × 10⁴, 2 × 10⁵ s. We vary the composition by adding a small mass fraction of Nd. Apart from the ejecta mass, magnetic field strength, and composition, the model parameters are the same as in Section 3.1.1. The model parameters are listed in Table 2.

The parameters for the pulsar again imply t_gw = 495 s; the pulsar luminosity is ${L}_{0}\in \{2.88,288,28800\}\times {10}^{43}$ erg s⁻¹. Tables 3 and 4 display the luminosity at day 1 for the SAFe models with 2 × 10⁴ and 2 × 10⁵ s, respectively. Tables 5 and 6 show the same data for the SAFeNd models. For the lowest magnetic field strength, 10¹² G, and lowest remnant cutoff time, increasing the mass produces an upward trend in the luminosity at day 1. For these parameters, the r-process decay energy competes with the pulsar luminosity in setting the kilonova bolometric luminosity at day 1. This is discernible in Figure 2, where increasing the ejecta mass can be seen to mask the pulsar peak around day 1. With a remnant lifetime of 2 × 10⁵ s, increasing the ejecta mass from 0.001 to 0.003 M_⊙ lowers the luminosity at day 1. The diminished luminosity at day 1 indicates delay in the emission from increased optical depth. For the higher magnetic field strengths, the luminosity at day 1 only decreases when more mass is added, resulting from the increased optical depth. The dimming and delaying of the peak luminosity for the models with B = 10¹³ G and ${t}_{\mathrm{cut}}=2\times {10}^{5}\,{\rm{s}}$ can be seen in Figure 3(c).

Zoom In Zoom Out Reset image size

The introduction of the Nd mass fraction, ${\chi }_{\mathrm{Nd}}={10}^{-4}$, does not appear to significantly impact the bolometric luminosity, with respect to the pure Fe models. At sufficiently late time, the bolometric luminosity is set by the r-process decay rate. However, the spectrum at later times tends to be redder; in particular for the low-pulsar luminosity, low-remnant time cutoff spectrum shown in Figure 3(c), the spectral features corresponding to the blue transient are systematically dimmer with ${\chi }_{\mathrm{Nd}}={10}^{-4}$. In Figure 3(d), the brighter, longer duration pulsar luminosity appears to sustain the blue portion of the spectrum further in time.

We have tabulated the day 1 and day 7 bolometric luminosities and UBVRIJHK broadband magnitudes for these models in Tables 10–13 of the Appendix. The intent of tabulating the data at an early and later time is to identify the effect of variations in the model on evolution of the luminosity and magnitudes. As an example, we can see that the sensitivity of the day 7 model data to remnant cutoff time depends on the pulsar magnetic strength. For instance, in Table 12, for B = 10¹² G and M_ej = 0.001 M_⊙, the B-band magnitude is 32.9 at day 7 for models with either ${t}_{\mathrm{cut}}=2\times {10}^{4}$ or 2 × 10⁵ s. The same models with B = 10¹³ G have day 7 B-band magnitudes of 32.9 and 25.3 for ${t}_{\mathrm{cut}}=2\times {10}^{4}$ and 2 × 10⁵ s, respectively.

The morphology of the dynamical ejecta of the 2D models is "model A" used by Wollaeger et al. (2018). The wind superimposed on the model A ejecta has the properties described in Section 3.1.1 for B = 10¹² and 10¹³ G, but without a remnant cutoff time. These field strengths correspond to initial pulsar luminosities of 2.88 × 10⁴³ and 2.88 × 10⁴⁵ erg s⁻¹, respectively. Figure 4 has bolometric luminosity for the range of angular views as shaded regions for each model.

Zoom In Zoom Out Reset image size

These luminosities are "isotropic equivalents," where each is divided by the solid angle per view and multiplied by 4π.

The angular views are divided into 54 polar angular viewing ranges of equal solid angle. The dimmest light curves (or lowest edge of the shaded regions) correspond to views most closely aligned with the merger plane ("edge-on"), which are most obscured by Nd in the dynamical ejecta. The axial view is over an order of magnitude brighter than for the edge-on view for the model with B = 10¹³ G.

Strong viewing-angle dependence is also seen in broadband magnitudes. Figure 5 has plots of UBVRIJHK absolute magnitudes for each model. The viewing-angle dependence is stronger for bluer bands, as expected. Before 2 days, the side view (dashed) is not sensitive to the strength of the pulsar luminosity.

Zoom In Zoom Out Reset image size

For the 1D fallback models, we vary ${v}_{\max }$, $\dot{M}$, and M_ej; otherwise the model parameters of Li et al. (2018) for the fallback luminosity source are adopted, as shown in Equation (8). The ejecta masses, velocity, and composition are the same here as in Section 3.1.2, permitting direct comparison between these fallback models and a subset of the pulsar models. The model parameters are given in Table 7.

Figure 6 displays bolometric luminosity versus time for the different accretion rates. Increasing the ejecta mass acts to dim the light curve, since the fallback source is more obscured at higher optical depth. This indicates that much of the kilonova luminosity in these models is derived from radiated accretion energy.

Zoom In Zoom Out Reset image size

In the Appendix, UBVRIJHK band data at day 1 and day 7 are provided for these fallback models in Tables 14 and 15. The accretion rate impacts the rate of decline in the optical bands. For instance, for M_ej = 0.001 M_⊙ and ${v}_{\max }/2=0.3c$, accretion rates of $\dot{M}=0.001$, $\dot{M}=0.003$, and $\dot{M}=0.01$ M_⊙ s⁻¹ drop by ∼4.2, 3.5, and 2.8 mag in the V-band, respectively. Likewise, the R-band drops by ∼3.9, 3.2, and 2.7 mag, respectively. The J- through K-bands are relatively insensitive to the fallback source at day 1, but their brightness is more affected by day 7, where the variation with accretion rate is > 1 mag. In other words, at early time, the JKH magnitudes appear to be mainly set by M_ej, while at later time, there is a systematic increase in magnitude with increasing accretion rate. For the high-velocity models, day 1 is later in the evolution of the light curves, so these trends are not reflected in the broadband data.

For the fallback models in 2D, we again use the "model A" morphology, as in Section 3.1.3. The wind superimposed on the model A ejecta has the properties described in Section 3.2.1 for $\dot{M}=0.001$ and 0.003 M_⊙ s⁻¹, and with a wind mass of 0.001 M_⊙. Figure 7 has bolometric luminosity for the range of angular views as shaded regions for each model. As in Figure 4, the values are shown as isotropic equivalents. The angular views are again divided into 54 polar angular viewing ranges of equal solid angle, with edge-on views being dimmest. Relative to the 2D pulsar models, the variation with respect to viewing angle does not change as substantially when increasing the source luminosity. This is partly explained by the fallback luminosity scaling linearly with accretion rate, whereas the pulsar source luminosity scales as the square of the magnetic field strength. As in the 2D pulsar models, increasing the remnant source luminosity affects the rise and decline times for near-on-axis angular views of the ejecta.

Zoom In Zoom Out Reset image size

Strong viewing-angle dependence is again seen in broadband magnitudes. Figure 8 has plots of UBVRIJHK absolute magnitudes for each model. The viewing-angle dependence is stronger for bluer bands, as expected. Before 2 days, the side view (dashed) is not sensitive to the strength of the pulsar luminosity. Unlike the 2D pulsar models, the U-band does not significantly shift in time and peak brightness when going to a higher remnant source, which is due to the more modest increase in luminosity relative to the pulsar models.

Zoom In Zoom Out Reset image size

We can compare our two different additional power sources to observations of AT 2017gfo using the results of Li et al. (2018) and Piro et al. (2019) as a starting point. Although there is a qualitative agreement with those results, there are some differences in our conclusions. For the polar ejecta mass of 0.001 M_⊙ given by Li et al. (2018), our ejecta model requires a slightly lower pulsar luminosity and somewhat higher wind outflow speed than previously published: 2 × 10⁴⁴ instead of 3.4 × 10⁴⁴ erg s⁻¹, and 0.45c instead of 0.35c. Additionally, we find that this model fits the later time bolometric luminosity better with a remnant cutoff time of about ∼10 observer days. From the modeling perspective, some differences here include the use of multifrequency opacity tables and a detailed high-latitude r-process heating rate (corresponding to the "Wind 1" composition of Wollaeger et al. (2018)) derived from WinNet (Winteler 2012; Winteler et al. 2012). For the fallback model, similar to the polar model of Li et al. (2018), we find an accretion rate of 0.003–0.0045 M_⊙ s⁻¹ gives bolometric luminosity in the range of the observation, but with a velocity of ∼0.4c instead of 0.25c. Parameters for the pulsar and fallback models are given in Tables 8 and 9, respectively.

Figure 10 has bolometric and broadband luminosity versus time for our attempt to match AT 2017gfo with a pulsar and a fallback model. Despite the decent agreement in bolometric luminosity in Figure 10(a), it is apparent in 10(b) that this model produces blue emission that does not decay quickly enough with respect to the observation. Figure 10(c) has the light curve for the fallback model. The fallback model bolometric luminosity does not fit the data as well at intermediate times, implying the decline in the source luminosity is too rapid to fully account for this emission, as found by Li et al. (2018). However, the trends in the B- and V- bands are much more similar to those of AT 2017gfo, relative to the pulsar model. The I-band of the fallback model is not bright enough at later times, which may be part of the disk.

Zoom In Zoom Out Reset image size