Superluminal motion of a relativistic jet in the neutron star merger GW170817

References

References

Observations, Data processing & Basic analysis

In order to establish the size and morphology of the faint radio afterglow of GW170817, we obtained Director’s Discretionary Time (program ID BM469) to observe with the High Sensitivity Array (HSA). The HSA antennas included the ten Very Long Baseline Array (VLBA) dishes, the phased Karl G. Jansky Very Large Array (VLA), and the Green Bank Telescope (GBT), although not all stations were present in all observations. The maximum baseline was typically 7,500–8,000 km.

We observed GW170817 with the HSA over four epochs between 2017 September – 2018 April. Each epoch consisted of 2–4 observations carried out over a period of up to 10 days, with approximately three hours of on-source time on GW170817 per day. The choice of the observing frequency was informed by the results from the VLA monitoring of the radio light curve, the desired angular resolution, and the ease of scheduling on the telescopes. In all epochs, a total bandwidth of 256 MHz was sampled in dual polarisation at 2-bit precision. Depending on the observing frequency, the recorded bandwidth was broken into eight 32 MHz wide bands, or two 128 MHz wide bands. A summary of the observations is given in Table 1.

The first epoch was undertaken at L band (central frequency 1550 MHz) 37 – 38 d post-merger. No fringes were seen on the GBT on one of the two observing days due to an unknown technical issue, considerably reducing overall sensitivity at this epoch. The second epoch was carried out in S band (central frequency 3200 MHz), 51 – 52 d post-merger. However, a misconfiguration of the VLA correlator on both days meant that phased VLA data was practically unusable, and hence sensitivity was severely impacted. The third epoch was observed at C band (central frequency 4540 MHz) 72 – 79 d post-merger. The fourth epoch was likewise observed at C band 227–236 d post-merger, utilising only the VLBA and VLA as the GBT was unavailable.

Each observation was structured around an 8 minute cycle as follows. We used the source J1258-2219 (a $\sim$1 Jy flat-spectrum source, separated by 2.8 degrees from GW170817) as the primary delay and gain calibrator, visiting it twice per cycle during first three epochs, and once per cycle in the fourth epoch observations. J1312-2350, a 20 mJy source separated by 0.8 degrees from GW170817, was used as a secondary phase calibrator, and was visited once per cycle in the first three epochs, and twice per cycle in the fourth epoch observations. J1258-2219 was additionally used to determine phase solutions for the VLA once per cycle. A single scan on 3C286 was included at the end of each observation to allow flux calibration of the commensally-recorded VLA interferometer data. For the C band (4.5 GHz) epochs only, we included three scans on the blazar OQ208 (B1404+286) over the course of each observation to enable polarization calibration to be determined and applied.

2 VLBI Data Processing

We followed standard data reduction procedures for HSA data using the AIPS software package . For all calibration steps that involve a sky source (fringe-fitting, leakage, and self-calibration) we used a model of the source that was iteratively refined over several passes of the entire data reduction pipeline.

The data was loaded using “FITLD” and a priori amplitude corrections were applied using “ANTAB” and “ACCOR”. We note that an issue with the VLA automatic gain control was uncovered whereby the phased VLA data exhibited large short-term amplitude variations; this could be (and was) largely mitigated by using a per-integration solution for the auto-correlation based corrections with “ACCOR”, but small residual variations which were weakly detrimental to sensitivity remained. This problem was fixed prior to the fourth observational epoch. “CLCOR” was used to correct for parallactic angle rotation and to apply the most accurate available values for Earth Orientation Parameters. TECOR was used to correct for ionospheric propagation effects, using the “igsg” model available from ftp://cddis.gsfc.nasa.gov/gps/products/ionex. We then calibrated the time-independent delays and the antenna bandpass using “FRING” and “BPASS”; in the first two epochs using a scan on the primary calibrator J1258-2219, while in the third and fourth epochs we used OQ208.

For the third epoch at 4.5 GHz only, we calibrated the cross-polar delays and instrumental polarization leakage using the tasks “FRING” and “LPCAL” and the source OQ208. This step was essential due to the large ($\sim$30%) leakage at the GBT at this frequency. “LPCAL” solves for a single leakage value per subband, while the GBT polarisation leakage varies across the 128 MHz subband; accordingly, we split each 128 MHz subband into 4$\times$32 MHz subbands to allow a coarse frequency dependence to the leakage solutions.

We solved for time dependent delays using “FRING” on the primary gain calibrator J1258-2219, followed by self-calibration on this source using “CALIB”, obtaining a single solution per subband, per scan. Finally, we improved the phase calibration using self-calibration on the secondary gain calibrator J1312-2350, deriving a single frequency-independent solution per scan.

At each stage, the solutions from the SN table were applied to the CL table using “CLCAL”. The final CL table was applied to the target using “SPLIT”. The target was then exported in UVFITS format using “FITTP” and imaged using “difmap”.

3 VLA/VLBI Interferometric data processing

We processed using VLA cross-correlated data (with the WIDAR correlator) using a custom-developed pipeline, which incorporates manual flagging, and standard interferometric data calibration techniques in CASA. The imaging was done with the CASA task clean with natural weighting, choosing an image size of 4096 pix $\times$ 4096 pix and cell size of 0.5 arcsec.

The VLA-only data gives the GW170817 flux densities of $56\pm 8$ $\mu$Jy beam${}^{-1}\,$, $54\pm 8$ $\mu$Jy beam${}^{-1}\,$and $45\pm 7$ $\mu$Jy beam${}^{-1}\,$for the three observations of the third epoch at 4.5 GHz. All three observations combined give $55\pm 5$ $\mu$Jy beam${}^{-1}\,$. For the four observations of the fourth epoch, the flux density values are $55\pm 8$ $\mu$Jy beam${}^{-1}\,$, $46\pm 8$ $\mu$Jy beam${}^{-1}\,$, $48\pm 6$ $\mu$Jy beam${}^{-1}\,$and $46\pm 6$ $\mu$Jy beam${}^{-1}\,$, while all four observations combined give $48\pm 4$ $\mu$Jy beam${}^{-1}\,$.

4 Flux comparison between the VLBI and VLA interferometric data

A comparison between the flux densities measured in the VLA-only interferometric data and those measured in the VLBI data (see Extended Data Table 1) implies that, within 1$\sigma$ uncertainties (typically 10% of the source flux density), no flux is being resolved out in the VLBI data.

5 Model fits and parameter estimations

Difmap was first used to produce a ”dirty” (un-deconvolved) image from the concatenated data from each epoch, as well as the individual observations within each epoch. In the first two epochs, there was substantial loss of sensitivity due to technical issues and the source was not detected. We place 5 upper limits of 40 $\mu$Jy beam${}^{-1}\,$(1.6 GHz, day 38) and 60 $\mu$Jy beam${}^{-1}\,$(3.2 GHz, day 52), respectively on the flux densities of GW170817, and do not consider these epochs further.

In the third and fourth epochs, a radio counterpart to GW170817 can clearly be seen in the dirty images for the concatenated datasets, and the source can also be seen (albeit at low S/N) in the individual observations. Initially, we fit the data in the visibility plane using a single circularly symmetric gaussian model component. Whilst likely an over-simplification of the true source structure, this has the advantage of being fast and simple to fit, while providing an accurate estimate of the flux centroid position. After model fitting, we read the resultant clean image into AIPS and used the task JMFIT to fit an elliptical gaussian in the image plane. Compared to model fitting, this has the advantage of providing well-constrained estimates of the uncertainty of the key parameters of interest. In the third epoch (75 days), the best-fit values of flux density and position are $58\pm 5$ $\mu$Jy beam${}^{-1}\,$and RA=13:09:48.068638(9), Dec=-23:22:53.3909(4) respectively. The uncertainties given here are purely statistical; we consider systematic contributions in the following sections. The best-fit size was a full-width half-maximum (FWHM) of 0.0 mas; i.e., the source was modeled as a point source. At day 230, the best-fit values of flux density and position were $48\pm 6$ $\mu$Jy beam${}^{-1}\,$, RA=13:09:48.068831(11) Dec= -23:22:53.3907(4) respectively, and the best-fit deconvolved size was 0.7 mas, although an unresolved source could not be excluded. The images of the source at 75 days and 230 days are shown in Extended Data Figure 1.

6 Estimating systematic contributions to flux density and position uncertainties

The absolute calibration of flux densities in VLBI maps is typically challenging due to the fact the sources compact enough to be visible at milliarcsecond resolution typically show evolution on a timescale on months to years. In cases where only a priori amplitude calibration can be performed, the accuracy of the flux density scale of a VLBI image is typically assumed to be of order 20%. In this case, we are able to use the contemporaneous VLA data to establish an absolute flux density scale, using the calibrator sources J1312-2350 and J1258-2219 (under the assumption that these sources do not have significant structure on scales larger than that resolvable by our VLBI observations). After adjusting the VLBI amplitude scale to produce the closest match to these two sources, the residual differences are typically 10% for each observation, and hence systematic uncertainties on our measured values of flux density for GW170817 are comparable to our statistical uncertainties.

Similarly, for our image centroid positions, we must consider the possibility of systematic position shifts between epochs due to calibration errors, in addition to the limiting precision attainable based on the image resolution and S/N. We neglect systematic errors due to the uncertainty in the calibrator reference position, since this would affect both epochs equally. Given the relatively close proximity of our calibrator source J1312-2350 to GW170817 (0.8 degrees), we expect any systematic errors that vary between epochs to be at most a small fraction of the synthesized beam size. Astrometric simulations suggest a typical systematic error for a single observation with the VLBA of 0.07 mas in right ascension and 0.25 mas in declination for our observing conditions (declination $-26$ degrees, angular separation 0.8 degrees). However, these simulations do not include the effect of the ionosphere, which could treble the systematic error at an observing frequency of 4.5 GHz under typical conditions. Countering this somewhat, our epochs consist of 3–4 observations spread over $\sim$7 days, and systematic errors (in particular those due to the ionosphere) are likely to be only weakly correlated over this timescale. Based on these considerations, we estimated the systematic position uncertainty to be 0.15 mas in R.A. and 0.5 mas in declination, and added this value in quadrature with the formal position fit errors at each epoch.

In order to verify this expectation, we repeated the data reduction for the third and fourth epochs after shifting the phase center of our target field to the position of the NGC 4993 low-luminosity AGN. This source is separated by 10.3 arcseconds from GW170817, and hence falls outside the field of view of the phased VLA; accordingly, the VLA was flagged before imaging. The positions obtained for the AGN have a separation of 0.05 mas in right ascension and 0.5 mas in declination (see Extended Data Figure 2). This is consistent with both their statistical uncertainties and our estimate for the systematic errors derived above. The AGN flux density is consistent with a constant value ($0.25\pm 0.02$ mJy and $0.29\pm 0.03$ mJy in the third and fourth epochs respectively, where the 1$\sigma$ uncertainties are purely statistical).

Comparison between the VLBI data and synthetic images

In order to compare the generated models with our VLBI data, we converted the simulated images (example images shown in Figure 3; for details of the simulations see the next section) into difmap models consisting of point sources at the center of each non-zero pixel in the simulated image, and performed model fitting in the visibility plane. The rotation, translation, and total flux density of the image were taken as free parameters, although we used the approximate positions and flux densities from our earlier fitting of circular gaussian components to restrict the ranges of parameter values over which we searched. For each model, we recorded the $\chi^{2}$ obtained at the best-fit values for rotation, translation, and total flux density.

Because the signal-to-noise of each individual visibility measurement is very low, determining the increase in $\chi^{2}$ that indicates a significant discrepancy between models is not straightforward. Previous authors have often relied on visual inspection of images and visibility data in order to determine model goodness-of-fit. Due to the low signal-to-noise ratio of our target image, we have taken a different approach. First, we used an image plane fit to determine the position errors in the image plane using the dataset fit with a circular gaussian component, which is a well-understood process. Second, we perturbed the position of the circular gaussian model component by up to $\pm$3$\sigma$ in right ascension and $\pm$3$\sigma$ in declination, and recorded the change in $\chi^{2}$ at offsets of 1, 2, and 3$\sigma$. A consistent increase in $\chi^{2}$ was seen regardless of the direction of the positional perturbation. Finally, we fitted other models based on the hydrodynamic simulations to the data and recorded the $\chi^{2}$ in each case. The reference positions for a given model were allowed to vary between the day 75 and day 230 datasets by up to the amount of our estimated systematic position uncertainty of 0.15 mas in R.A. and 0.5 mas in Declination. By comparison to the set of $\chi^{2}$ values obtained from the perturbed circular gaussian fits, we estimated the consistency of each hydrodynamic model with the best-fit circular gaussian model.

In addition to fitting the actual synthetic images, we first produced an estimate of the maximum source extent, by finding the largest circular and elliptical gaussian sources that produced a $\chi^{2}$ that did not deviate by more than 1 $\sigma$ from the best circular gaussian fits. For the epoch at day 75 and day 230, the largest circular gaussian source was 1.1 and 1.2 mas in diameter respectively. The best-fit elliptical gaussian converged to an unphysical one-dimensional source for each epoch, with an upper limit on the major axis of 12 mas and 9 mas for day 75 and day 230 respectively. In both cases the best-fit position angle was approximately aligned with the beam major axis and hence approximately perpendicular to direction of source motion. Tighter limits on the maximum size can be obtained if the axial ratio of the elliptical gaussian source is constrained to a physical value: for instance, in the case of the day 230 dataset, the largest source permitted with an axial ratio of 4:1 has size 3.9 mas $\times$ 0.9 mas. Hence, the source size parallel to the direction of motion is relatively well constrained.

None of the synthetic images produced a $\chi^{2}$ significantly better than a simple circular gaussian in either epoch (unsurprising, given that the source was consistent with being unresolved in both cases). Generally, we found that as the positional offset between days 75 and 230 increased, the ”best-fit” source size at day 230 also increased and was often inconsistent with the observed source compactness. This disfavoured models at low viewing angles. Conversely, models at large viewing angles were incapable of producing a sufficiently large positional offset.

The best-fitting model (narrow jet viewed at 0.35 radians, model A1 in Figures 2 and 3) was able to produce the expected positional shift between epochs: with a constant reference translation and rotation, it produced an acceptable fit to both the day 75 epoch ($\chi^{2}$ increase equivalent to a 0.9$\sigma$ position offset for the circular gaussian) and the day 230 epoch ($\chi^{2}$ increase equivalent to a 1.3$\sigma$ position offset for the circular gaussian). Among the other models, only one (model B, the very narrow jet viewed at 0.3 radians) remained consistent within 2$\sigma$ for both epochs. For all other models, the discrepancy with the best-fit circular gaussian exceeded 2$\sigma$ in one or both epochs. As can be seen in Figure 2, models A1 and B are also those that best fit the light curve.

Numerical hydrodynamic simulations

To characterize the properties of different models we carry out relativistic hydrodynamical simulations of various setups, followed by a post processing numerical calculation of their afterglow light curve and observed images at 75 and 230 days. In particular we run different type of models to see which have the potential to fit the entire data set of both the light curve and the image characteristics, i.e. the flux centroid movement and the image size constraints.

Our setup includes three components: the jet, a core of cold massive ejecta and a fast ejecta tail. Each component of the ejecta expands homologously and has a density profile of

where the normalization $\rho_{0}$ is determined by the total ejecta mass and $\alpha$ and $\beta$ which differ between models, dictate the radial and angular structures, respectively. However, our main focus was on scanning the jet’s properties such as luminosities, opening angles, injection and delay times. While some of the jets successfully break out from the ejecta if their properties allow, others may be choked inside it. We ran about ten different models, here we present four representative models that demonstrate how the different characteristics of the jet affect the observed outcome. The first two models are narrow jets which are found to fit all the observed characteristics-the gradual rise of the flux, the short plateau at the peak followed by a fast decline and the large flux centroid motion between the two image epochs. In addition we also present a wider successful jet and a choked jet. The full setup is given in Extended Data Table 2.

A full description of the hydrodynamic simulations is given in our previous work. Briefly, for each model we use three different simulations. The first one which includes the jet propagation inside the core ejecta is performed in 3D to avoid the numerical plug artifact. The second simulation includes the outflow evolution inside the tail ejecta and after breaking out of it until reaching the homologous phase. This simulation is modeled in 2D as previously we showed that after breakout the plug artifact is no longer a concern, and 2D and 3D simulations become similar. Finally, the third simulation begins when the afterglow becomes important and ends after it decays.

For the relativistic hydrodynamical simulation we use the public code PLUTO v4.0 with an HLL Riemann solver and we apply an equation of state with adiabatic index of 4/3. The setup of models ${\it A}$ and ${\it B}$ is as follows. The grid setup of the first 3D Cartesian simulation has three patches in x and y axes and two patches on the z axis. On x and y the inner patch spans from $-2\times 10^{8}\,{\rm cm}$ to $2\times 10^{8}\,{\rm cm}$ with 30 uniform cells. The outer patch is from $|2\times 10^{8}\,{\rm cm}|$ to $|3\times 10^{10}\,{\rm cm}|$ with 400 cells that are distributed logarithmically. On the z-axis the first patch is uniform from $4.5\times 10^{8}\,{\rm cm}$ to $10^{10}\,{\rm cm}$ with 200 cells followed by a logarithmic patch of 400 cells until $4\times 10^{10}\,{\rm cm}$. We convert the 3D output of the first simulation to an axisymmetric grid, which is the initial setup of the second simulation for which the setup is as follows. The first two patches on r and z axes correspond to the 3D setup. We add another patch on each axis from $3\times 10^{10}\,{\rm cm}$ ($4\times 10^{10}\,{\rm cm}$) on the r (z) axis, to $6\times 10^{11}\,{\rm cm}$ with 1200 logarithmic cells.

For the third simulation which includes two patches on each axis, we use the output of the second simulation. The first patch corresponds to the second simulation grid with 800 uniform cells until $6\times 10^{11}\times R\,{\rm cm}$ on each axis. The second patch on each axis stretches to $10^{14}\times R\,{\rm cm}$ with 6000 logarithmic cells. As the simulation is dimensionless, we use $R$ as a scaling length factor, $R$ also determines the ISM density which is set to be $\rho_{ISM}=5\times 10^{-12}\rm{gr}(R\times\,{\rm cm})^{-3}$ in simulation ${\it A}$ and $\rho_{ISM}=8\times 10^{-12}\rm{gr}(R\times\,{\rm cm})^{-3}$ in simulation ${\it B}$. Each viewing angle fit requires a different $R$. The best fits for $\theta_{obs}=0.25,0.35,0.45$ in simulation ${\it A}$ are obtained at $R=3\times 10^{5},1.7\times 10^{5},8.3\times 10^{4}$, respectively, and for $\theta_{obs}=0.3$ in simulation ${\it B}$ it is $R=5\times 10^{5}$.

The setup of simulations ${\it C}$ and ${\it D}$ was described previously (simulation ${\it D}$ is identical to the successful jet scenario, except for the engine time), and the only difference here is that for the outer patch in the third part we use a high resolution of 4000 cells rather than 2500 cells originally. The scaling of the third part of the simulation is determined by $n=4\times 10^{-2}~\rm{cm}^{-3}$ and $n=4.5\times 10^{-3}~\rm{cm}^{-3}$ in ${\it C}$ and ${\it D}$ respectively.

Finally, we verify that each of the three simulation meets the required resolution to reach convergence. We first compare the resolution of the first two simulations, from the jet launch until reaching the homologous phase, with previously-published simulations for which convergence tests have been taken. The resolution of the 3D simulation which handles the jet propagation inside the ejecta is comparable with that of the inner parts of theirs. The sequential 2D simulation has naturally a higher resolution compared with the outer parts of the 3D grid presented previously. For convergence of the third part in which the outflow interacts with the ISM, we perform another set of simulations with 2/3 the resolution aforementioned. We find that both the light curves and the images for the relevant viewing angles remain essentially unchanged with the increase in resolution.

Details of the simulation that provides the best fit to the data

Our simulation that provides the best fit to the data is of a jet with a 0.08 rad (4${}^{o}\,$) opening angle, at the time of light curve peak, that is observed at a viewing angle of $\theta_{\rm obs}=0.35$ rad (20${}^{o}\,$). In this simulation a relativistic jet is injected into the sub-relativistic merger ejecta. The jet is followed during its propagation through the ejecta, the formation of the cocoon and the breakout of the jet and the cocoon from the dynamical (sub-relativistic) ejecta. The simulation then continues to follow the interaction of the outflow (jet+cocoon) with the ISM. When this interaction starts, the jet opening angle is 0.04 rad. The cocoon dominates the observed radio emission during the first $\sim$60 days, and after this time the jet dominates. The jet expands sideways slowly during its interaction with the ISM, reaching an opening angle of 0.08 rad after $\sim$150 days at the light curve peak. On day 75, the Lorentz factor of the observed region is $\Gamma\approx 4$, which steadily drops to $\Gamma\approx 3$ by day 230.

Constraining the jet energy and the external density

The gamma-ray signal from GW170817 had an isotropic equivalent energy of $5\times 10^{47}$ erg. The afterglow suggests that this energy is not representative of the jet energy. This is consistent with models for the gamma-ray emission. Therefore, in order to constrain the jet energy and external density, we use the constraints on the geometry of the outflow together with the observed afterglow light curve to constrain the outflow energy. We use the standard afterglow model, where a narrow ultra-relativistic jet drives a blast wave into the external medium which radiates in synchrotron emission to produce the radio and X-ray afterglow. Before interacting with the external medium the jet has an initial Lorentz factor $\Gamma_{0}$. This is also the initial Lorentz factor of the blast wave that it drives, which is constant at first until the blast wave accumulates enough mass and starts decelerating. Its initial opening angle, $\theta_{j,0}$, is also constant until the Lorentz factor drops to $\sim 1/\theta_{j,0}$. At this point, if $\theta_{j,0}<0.05$ rad it starts spreading sideways rapidly until $\theta_{j,0}\sim 0.05$ rad, at which point it starts spreading sideways more slowly. We have direct constraints only of $\Gamma$ and $\theta_{j}$ near the time of the peak of the light curve. We therefore can only put a lower limit on the initial Lorentz factor, $\Gamma_{0}>4$, and an upper limit on the initial opening angle $\theta_{j,0}<0.1$ rad. Moreover, given the fast spreading of the jet if $\theta_{j,0}<0.1$ rad and $\Gamma<1/\theta$, at the time that we observe the jet its opening angle is expected to be $\theta_{j}\approx 0.05-0.1$ even if initially $\theta_{j,0}\ll 0.1$ rad and its Lorentz factor is $\Gamma_{0}\gg 4$. The Lorentz factor and the time of the peak provide a relation between the ambient medium density (assumed to be constant) and the jet isotropic equivalent energy: $E_{iso}\sim 10^{52}\frac{n}{3\times 10^{-4}{\rm cm^{-3}}}$ erg. The flux is extremely sensitive to the Lorentz factor and we can use its value at the peak to constrain the density and the fraction of the internal energy that goes to the magnetic field, $\epsilon_{B}$: $\frac{n}{3\times 10^{-4}{\rm cm^{-3}}}\left(\frac{\epsilon_{B}}{10^{-3}}\right)^{0.47}\sim\left(\frac{\Gamma}{3.5}\right)^{5.9}$, where we assume that 10% of the internal energy goes to the accelerated electrons ($\epsilon_{e}=0.1$) and that their distribution power-law index is $p=2.16$. Allowing the least constrained parameter, $\epsilon_{B}$, to vary between $10^{-2}$ and $10^{-5}$ we find that the circum-merger density is $10^{-4}-5\times 10^{-3}{\rm cm^{-3}}$ and the jet isotropic equivalent energy is $E_{iso}\sim 3\times 10^{51}-10^{53}$ erg. Since the jet opening angle at this time is 0.05–0.1 rad and it contains a significant fraction of the total energy of the relativistic outflow (jet+cocoon), we find that the energy deposited by the merger in relativistic ejecta is $10^{49}-10^{50}$ erg. The confirmation of a successful jet in GW170817 also implies high isotropy of the magnetic field.

Methods References