FIRST J1419+3940 as the First Observed Radio Flare from a Neutron Star Merger

We first probed whether J1419+3940 is consistent with a neutron star merger origin. We adopted the collected radio observations of Law et al. (2018, see their Table 2) and, in addition, the observation of Marcote et al. (2019) with the European VLBI Network. These observations have been carried out at a variety of radio frequencies. For an easier interpretation of the results, we converted the observed fluxes to their expected values at either 0.3, 1.4, or 3 GHz. All observations had a frequency close to one of these three values. Each result was converted to the closest of these three frequencies by assuming that the flux scales with frequency as ${\nu }^{-(p-1)/2}$, where p is the power-law index of the distribution of the accelerated electrons' Lorentz factors, which we take to be p = 2.5 (Piran et al. 2013). The measured fluxes and upper limits as functions of time are shown in Figure 2.

Zoom In Zoom Out Reset image size

We modeled the radio emission from the interaction of neutron star merger ejecta and the circum-merger medium following the prescription of Piran et al. (2013). Based on the numerical outflow simulations of Fernández et al. (2017), we adopted a broken power-law velocity probability density for this outflow, with peak velocity v₀ ≈ 0.2, and indices α₁ = 5 and α₂ = −10 below and above the peak, respectively.

We considered a uniform circum-merger medium with baryon number density n. We parameterized the expansion of the ejecta into this medium with the radial position R of the front of the ejecta from the merger site. At a given R, the front of the outflow accumulated $4/3\pi {R}^{3}n+{M}_{\mathrm{ej}}(\beta )$ mass, where M_ej(β) is the part of the ejecta mass with initial velocity ≥β. Let E(β) be the total kinetic energy of the part of the ejecta mass with initial velocity ≥β. Then the velocity β = β(R) as a function of R can be determined from energy conservation (Piran et al. 2013):

$$\begin{eqnarray}&&M(R){\left(\beta c\right)}^{2}\approx E(\geqslant \beta ).\end{eqnarray} \tag{ 1 }$$

We solved this equation numerically to obtain β(R), and used this to compute $R{\left(t)={\int }_{0}^{t}\beta (R\right)}^{-1}{dR}$, where time t is measured from the start of the outflow.

To determine the radio flux produced by the outflow, we computed two relevant characteristic frequencies following Piran et al. (2013): the typical electron synchrotron frequency

$$\begin{eqnarray}&&{\nu }_{{\rm{m}}}(t)\approx \text{1 GHz}\cdot {n}_{0}^{\tfrac{1}{2}}{\epsilon }_{{\rm{B}},-1}^{\tfrac{1}{2}}{\epsilon }_{{\rm{e}},-1}^{2}\beta {\left(t\right)}^{5}\end{eqnarray} \tag{ 2 }$$

and the self-absorption frequency

$$\begin{eqnarray}&&{\nu }_{{\rm{a}}}(t)\approx \text{1 GHz}\cdot {R}_{17}{\left(t\right)}^{\tfrac{2}{p+4}}{n}_{0}^{\tfrac{6+p}{2(p+4)}}{\epsilon }_{{\rm{B}},\ -1}^{\tfrac{2+p}{2(p+4)}}{\epsilon }_{{\rm{e}},-1}^{\tfrac{2(p-1)}{p+4}}\beta {\left(t\right)}^{\tfrac{5p-2}{p+4}}.\end{eqnarray} \tag{ 3 }$$

Here, _B and _e are the fractions of the total internal energy of the shocked gas carried by the magnetic fields and electrons. In these equations and below we adopted the notation.

The relation of these characteristic frequencies to the observational frequency ν_obs affects the expected flux (see Table 2 in Piran et al. 2013). To compare with available observational data we computed the expected radio flux at 0.3, 1.4, and 3 GHz (see Figure 2). For frequencies 1.4 GHz and 3 GHz, and for 0.3 GHz for t > 5 yr, we found that ${\nu }_{\mathrm{obs}}\gt {\nu }_{{\rm{a}}}\gt {\nu }_{{\rm{m}}}$. In this case we computed the expected flux using

$$\begin{eqnarray}&&{F}_{\mathrm{obs}}(t)={F}_{{\rm{m}}}(t)\cdot {\left(\displaystyle \frac{{\nu }_{\mathrm{obs}}}{{\nu }_{m}(t)}\right)}^{-\tfrac{p-1}{2}},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{F}_{m}(t)=500\,\mu {\rm{Jy}}\,R{\left(t\right)}_{17}^{3}{n}_{-1}^{3/2}\beta {\left(t\right)}^{1}{d}_{27}^{-2}.\end{eqnarray} \tag{ 5 }$$

For ${\nu }_{\mathrm{obs}}=0.3\,\mathrm{GHz}$ and t < 5 yr we found ${\nu }_{{\rm{a}}}\gt {\nu }_{\mathrm{obs}}\gt {\nu }_{{\rm{m}}}$. In this case we approximated the radio flux to be

$$\begin{eqnarray}&&{F}_{\mathrm{obs}}(t)={F}_{{\rm{m}}}(t)\cdot {\left(\displaystyle \frac{{\nu }_{a}(t)}{{\nu }_{m}(t)}\right)}^{-\tfrac{p-1}{2}}\cdot {\left(\displaystyle \frac{{\nu }_{\mathrm{obs}}}{{\nu }_{a}(t)}\right)}^{-\tfrac{5}{2}}.\end{eqnarray} \tag{ 6 }$$

This approximation is valid if the radio emission is before its peak, which we confirmed to be the case here.

We carried out a model fitting in which we identified the parameters of the radio emission model described above that best matched the observed data. We found the best match by minimizing the least-squares error of the model prediction versus the observations using gradient descent. Our best-fit parameters are ejecta mass M₀ = 0.005 M_⊙, characteristic velocity v₀ ≈ 0.3c, circum-merger density n = 5 cm⁻³, and a merger time in 1993. Our model also fit the fraction (Piran et al. 2013) of the kinetic energy in the shocked gas carried by electrons (_e) and magnetic fields (_B), for which we found ${\epsilon }_{{\rm{B}}}\approx {\epsilon }_{{\rm{e}}}\approx 0.2$.

The obtained radio light curves in comparison to observations are shown in Figure 2. We used three different radio bands, 0.3, 1.4, and 3 GHz, where observations were available. Observations with close but different frequencies were scaled to these values. We see that the fit closely follows the data for all three radio frequencies and all times.

We examined what information the host galaxy of FIRST J1419+3940 carries of its origin. FIRST J1419+3940 was detected in a dwarf galaxy, SDSS J141918.81+394035.8 (Abolfathi et al. 2018). The galaxy has an estimated stellar mass of ∼2 × 10⁷ M_⊙, and a star formation rate of ∼0.1 M_⊙ yr⁻¹. Neutron star mergers are typically expected to occur in more massive galaxies as the merger rate is primarily correlated to stellar mass due to the often long delay between star formation and binary merger (Berger 2014; Artale et al. 2020). However, detectability through a radio flare also means that the merger needed to reside in a relatively dense circum-merger medium, introducing selection effects that are difficult to account for. In addition, the original identification of the source is biased as Ofek (2017) searched for low-mass, star-forming galaxies, although this strategy by itself was likely not a major contributor to selection effects. Therefore, here we simply assessed whether similar galaxies can be sources of neutron star mergers.

We computed the expected merger rate in the IllustrisTNG cosmological hydrodynamical simulation (Springel et al. 2018; Nelson et al. 2019). IllustrisTNG uses a comprehensive galaxy formation model to make self-consistent predictions for galactic stellar mass buildup within galaxies. The model itself is tuned to match a number of observational constraints including the galaxy stellar mass function and cosmic star formation rate density (Vogelsberger et al. 2013; Torrey et al. 2014). One result from the IllustrisTNG simulation is star formation rate histories for all simulated systems. Since the overall stellar mass buildup in these simulations matches stellar mass functions across a broad redshift range, the star formation histories are likely reflective of the real universe.

We adopted the star formation rate histories from the IllustrisTNG simulation and convolved them with an assumed delay time distribution to calculate neutron star merger rates. For our merger delay time distribution, we used a power-law form with a cutoff time before which no neutron star mergers occur. We adopted the power law from Naiman et al. (2018) with an exponential slope of γ = −1.12 taken to match the observed SNIa delay time distribution (Maoz et al. 2012). We used a cutoff time of 100 Myr as a fiducial value (Safarzadeh & Berger 2019).

We evaluated the neutron star merger rate across the entire simulated volume, and used this to determine the fraction of neutron star merger events that occur within galaxies of interest for this Letter.

Considering galaxies with stellar mass <10⁸ M_⊙, we found that ∼1% of neutron star mergers are expected to occur in such galaxies. Therefore, while most neutron star mergers will occur in massive galaxies, the observed host galaxy could have been the host of a neutron star merger.

We note that other caveats remain, however. First, as the fiducial neutron star merger delay time of 100 Myr is longer than the typical duration of star formation epochs in dwarf galaxies, the host galaxy's unusually high star formation rate is likely unrelated to the binary's formation. Second, the transient's spatial offset from the galactic center is smaller than typically expected from neutron star mergers. Both of these properties are more typical for the case of long GRBs. Nonetheless, selection effects may be important here too. High star formation may correspond to high interstellar medium densities, necessary for a bright, observable radio flare. Similarly, a larger natal kick pushes the binary neutron star into a sparser environment, diminishing its observability.

Using the European VLBI Network, Marcote et al. (2019) measured the size of the radio source on 2018 September 18 to be about 1.6 ± 0.3 ± 0.2 pc, where the first error bar is the statistical uncertainty, and the second is the systematic uncertainty (derived here from their Table 1). Our best-fit model gives a source size of 1.2 pc on the same date. Therefore, we find a neutron star merger explanation to be consistent with the observed radio source size in 2018.

Both neutron star merger radio flare and long-GRB afterglow models are broadly consistent with the observed data (see Figure 2); therefore, the models' goodness of fit is insufficient to rule out one of them. However, we identified two features of the long-GRB afterglow light curve that are in tension with the data.

• 1.  
  Double-peak structure of early light curve. Law et al. (2018) found that the best-fit afterglow light curve assumes an off-axis GRB, i.e., that the GRB jet initially pointed away from Earth. Such off-axis afterglow light curves have a characteristic early double-peak structure (see the inset in Figure 2). The two peaks are due to the two GRB jets becoming "visible" in radio at different times. The observed early radio data of J1419+3940 show no such structure, but instead show a gradual decay. Such a decay could be observed if the double-peaked structured is observed such that the three observations are carried out serendipitously at the right "phases" of the two peaks; however, such a coincidence is unlikely. A more likely explanation is that the radio emission was indeed gradually decaying during this time. Gradual decay is still consistent with an on-axis GRB afterglow. In this scenario the GRB jet is pointing toward or away from Earth from the beginning; therefore, there will be no peaks in the radio light curve. However, while an on-axis GRB afterglow could explain the gradually decaying radio data observed at 1.4 GHz, it would predict a high initial flux that is inconsistent with the early upper limit observed at 0.3 GHz. Therefore, an on-axis GRB afterglow can be ruled out.
• 2.  
  Late-time temporal decay. The long-GRB afterglow radio light curve is expected to flatten after about a year, when the bulk of the shock-accelerated electrons in the outflow become nonrelativistic (Sironi & Giannios 2013; Law et al. 2018). This results in a deep-Newtonian temporal decay of the radio flux F ∝ t⁻¹. In contrast, we found that for the neutron star merger case the bulk of shock-accelerated electrons remains ultrarelativistic until about 40 yr after the merger. The velocity of the shock front is nonrelativistic for the whole period. This corresponds to a temporal decay of the flux F ∝ t^−1.7. This latter decay seem to be the case for the observed flux of FIRST J1419+3940. Near-future follow-up observations will further differentiate these two decay types. Given the expected differences in the temporal decay of the radio flux for the neutron star and afterglow models, the predicted flux difference for an observation during the summer of 2021 is around 500 μJ for the radio frequencies considered here, which is beyond the uncertainty expected from the power-law decay fit on the 2009–2018 observations. Therefore, we found that one additional observation in the summer of 2021 with VLA can substantially constrain the deep-Newtonian explanation of the late radio light curve.

As an alternative probe of the source's possible origin, we computed the a priori detection probability of a neutron star merger with the strategy that uncovered FIRST J1419+3940. We then compared this number to the a priori probability of a long-GRB afterglow origin.

We carried out Monte Carlo simulations in which we randomly placed neutron star mergers in space and time and checked whether the same detection strategy that was used to identify FIRST J1419+3940 would identify them.

For each neutron star merger we randomly drew (i) ejecta mass within [10⁻²M_⊙, 10⁻¹M_⊙] and ejecta velocity within 0.1c–0.2 c using independent uniform distributions; (ii) circum-merger density from the distribution of reconstructed densities for short GRBs (see Figure 6, top left, in Fong et al. 2015); (iii) location drawn from a uniform volumetric distribution within d_max = 108 Mpc, which is the maximum distance in the analysis of Ofek (2017); and (iv) time of the merger within the 50 yr period prior to the time of observation. Using these simulations we determined the average duration $\langle {\rm{\Delta }}t\rangle =1.1\,\mathrm{yr}$ for which a neutron star radio flare's flux is >4 mJy at 1.4 GHz, which is the threshold of the FIRST survey (see Law et al. 2018).

We adopted a neutron star merger rate of ${{ \mathcal R }}_{\mathrm{NS}}={900}_{-790}^{+2940}$ Gpc⁻³ yr⁻¹ at ∼90% confidence level (we obtained the expected value by averaging the four different analyses discussed in Section VII.C of Abbott et al. 2019). We note here that the higher end of this range (≲10³ Gpc⁻³ yr⁻¹) has been ruled out by radio transient searches (Levinson et al. 2002; Gal-Yam et al. 2006). We additionally took into account that the FIRST survey covered about f_FIRST ∼ 25% fraction of the sky, that the Karl G. Jansky Very Large Array sky survey (VLASS) used by Ofek (2017) covered about f_VLASS ∼ 50% fraction of the sky, and that the galaxy sample used in the study by Ofek (2017) had a completeness of f_galaxy ∼ 30% out to the covered d_max = 108 Mpc. With these factors the expected number of neutron star merger radio flares detected by the survey is

$$\begin{eqnarray}&&{{ \mathcal N }}_{\mathrm{NS}}=\displaystyle \frac{4}{3}\pi {d}_{\max }^{3}\,{{ \mathcal R }}_{\mathrm{NS}}\,\langle {\rm{\Delta }}t\rangle {f}_{\mathrm{FIRST}}\,{f}_{\mathrm{VLASS}}\,{f}_{\mathrm{galaxy}},\end{eqnarray} \tag{ 7 }$$

which gives us ${{ \mathcal N }}_{\mathrm{NS}}={0.2}_{-0.17}^{+0.64}$. Nevertheless, we note that the model parameters above, such as the neutron star's environment, ejecta mass, and and ejecta speed, are poorly constrained at present. Our uncertainty of ${{ \mathcal N }}_{\mathrm{NS}}$ is therefore larger than the statistical uncertainty quoted here.

For comparison we adopted here the relevant estimate of the long-GRB afterglow detection rate here from Law et al. (2018). They estimated $\langle {\rm{\Delta }}t\rangle =2000$ days to be the observed duration of FIRST J1419+3940 instead of our Monte Carlo simulation, and used a different source rate of ${{ \mathcal R }}_{\mathrm{GRB}}\sim 60$ Gpc⁻³ yr⁻¹ (see also Wanderman & Piran 2010; Goldstein et al. 2016), but otherwise carried out the same computation as we describe above. They found ${{ \mathcal N }}_{\mathrm{GRB}}=0.06$, which is likely significantly lower than the a priori neutron star merger probability.