An extreme magneto-ionic environment associated with the fast radio burst source FRB 121102
D. Michilli, A. Seymour, J. W. T. Hessels, L. G. Spitler, V. Gajjar, A. M. Archibald, G. C. Bower, S. Chatterjee, J. M. Cordes, K. Gourdji, G. H. Heald, V. M. Kaspi, C. J. Law, C. Sobey, E. A. K. Adams, C. G. Bassa, S. Bogdanov, C. Brinkman, P. Demorest, F. Fernandez, G. Hellbourg, T. J. W. Lazio, R. S. Lynch, N. Maddox, B. Marcote, M. A. McLaughlin, Z. Paragi, S. M. Ransom, P. Scholz, A. P. V. Siemion, S. P. Tendulkar, P. Van Rooy, R. S. Wharton, D. Whitlow
References
Methods
The analyses described here were based on the PRESTO, PSRCHIVE, and DSPSR pulsar software suites, as well as custom-written Python scripts for linking utilities into reduction pipelines, fitting the data, and plotting.
We observed using the Arecibo ‘C-band’ receiver (dual linear receptors), in the frequency range GHz, and the Puerto-Rican Ultimate Pulsar Processing Instrument (PUPPI) backend recorder. The full list of observations is reported in Extended Data Table 1. We operated PUPPI in its ‘coherent search’ mode, which produced s samples and MHz frequency channels, each coherently dedispersed to pc cm-3. Coherent dedispersion within each 1.56-MHz channel means that the intra-channel dispersive smearing is s even if the burst DM is pc cm-3 higher/lower than the fiducial value of pc cm-3 used in the PUPPI recording. The raw PUPPI data also provide auto- and cross-correlations of the two linear polarizations, which can be converted to Stokes I, Q, U, and V parameters in post-processing. Before each observation, both a test scan on a known pulsar (PSR B0525+21) and a noise-diode calibration scan (for polarimetric calibration) were performed.
Dedispersed time series with DM = pc cm-3, in trial steps of pc cm-3, were searched using PRESTO’s single_pulse_search.py, which applies a matched-filter technique to look for bursts with durations between s to s (for any putative burst that only has a single peak with width s, the sensitivity will be degraded by a factor of a few, at most). The resulting DM-time-S/N events were grouped into plausible astrophysical burst candidates using a custom sifting algorithm and then a dynamic spectrum of each candidate was plotted for human inspection and grading. We found 16 bursts of astrophysical origin, and used the DSPSR package to form full-resolution, full-polarization PSRCHIVE ‘archive’ format files for each burst.
On August 26, 2017, we observed FRB 121102 using the GBT ‘C-band’ receiver ( GHz, with dual linear receptors) as part of a program of monitoring known FRB positions. Observations were conducted with the Breakthrough Listen Digital Backend, which allowed recording of baseband voltage data across the entire nominal 4-GHz bandwidth of the selected receiver. Scans of a noise-diode calibration, of the flux calibrator 3C161 and of the bright pulsar PSR B0329+54 supplemented the observations.
In post-processing, a total intensity, low-resolution filterbank data product was searched for bursts with DM pc cm-3, using trial DMs in steps of 0.1 pc cm-3 and a GPU-accelerated search package to perform the incoherent dedispersion. We detected 15 bursts with S/N . Here we present the properties of just the two brightest GBT bursts in order to confirm the large RM observed by Arecibo and to quantify its variation in time. A detailed analysis of all GBT detections is presented in Gajjar et al. (in prep.). A section of raw voltage data (1.5 s) around each detected burst was extracted and coherently dedispersed to a nominal DM of 557.91 pc cm-3 using the DSPSR package. Final PSRFITS format data products have time and spectral resolutions of s and kHz, respectively.
Data analysis
We calibrated the burst ‘archives’ using the PSRCHIVE utility pac in ‘SingleAxis’ mode. This calibration strategy uses observations of a locally generated calibration signal (pulsed noise diode) to correct the relative gain and phase difference between the two polarization channels, under the assumption that the noise source emits equal power and has zero intrinsic phase difference in the two hands. This calibration scheme does not correct for cross-coupling or leakage between the polarizations. While leakage must be present at some level, the high polarization fraction, complete lack of circular polarization, and consistency of the test pulsar observations with previous work all give us confidence that calibration issues are not a significant source of error for the RM determination. In addition, the flux density of GBT observations was calibrated using the flux calibrator.
We initially performed a brute force search for peaks in the linear polarization fraction (Extended Data Fig. 3), and discovered rad m-2 in the Arecibo data. Each burst was corrected for Faraday rotation using the best-fit RM value for that burst. Residual variations in the resulting PA() were used to refine the initial values by fitting
where is the unit vector of the linear polarization, was used to fit the whole sample of bursts together, imposing a different RM per day and a different PA∞ per telescope. The results of these fits are reported in Table 1 and an example is shown in Fig. 2.
Applying the optimal RM value to each burst, we produced polarimetric profiles showing that each burst is consistent with being 100% linearly polarized after accounting for the finite widths of the PUPPI frequency channels (Fig. 1; Extended Data Fig. 2). In fact, the measured Arecibo bursts are depolarized to 98%, consistent with an uncorrected intra-channel Faraday rotation of
where is the speed of light, is the channel width, and is the central channel observing frequency. At 4.5 GHz this corresponds to 9∘, and the depolarization fraction is
We supplemented our above analysis with a combination of RM Synthesis and RMCLEAN (e.g. Extended Data Fig. 4). Ensuring the presence of minimal Faraday complexity is possible by integrating across the full bandwidth and taking advantage of a Fourier transform relation between the observed values and the Faraday spectrum (the polarized brightness as a function of RM). This approach is commonly known as RM Synthesis, and can be coupled with a deconvolution procedure (RMCLEAN) to estimate the intrinsic Faraday spectrum. While RM Synthesis and RMCLEAN can have difficulty in properly reconstructing the intrinsic Faraday spectrum under certain circumstances, the spread of clean components is a reliable indicator of spectra that contain more than a single Faraday-unresolved source.
As in previous studies, a search for periodicity in the burst arrival times remains inconclusive.
Determining the exact DMs of the bursts is complicated by their changing morphology with radio frequency. Measuring DM based on maximizing the peak S/N of the burst often leads to the blurring of burst structure and, in the case of FRB 121102, an overestimation of DM. We have thus chosen to display all bursts dedispersed to the same nominal DM from Burst #6 (Fig. 1 and Extended Data Fig. 1). Taking advantage of the narrowness of Burst #6, we estimated its optimal DM by minimizing its width at different DM trials. We measured burst widths at half the maximum by fitting von Mises functions using the PSRCHIVE routine paas (Table 1). These widths correspond to the burst envelope in the case of multi-component bursts.
Flux densities of the Arecibo bursts were estimated using the radiometer equation to calculate the equivalent RMS flux density of the noise:
where K and K Jy-1 are the system temperature and gain of the receiver, respectively, MHz is the observing bandwidth and s is the sampling time. GBT observations were instead calibrated using a flux calibrator as discussed above. Due to the complicated spectra of the bursts, we quote average values across the frequency band (Table 1).
The burst dynamic spectra in Extended Data Fig. 1 show narrow-band striations that are consistent with diffractive interstellar scintillations caused by turbulent plasma in the Milky Way. Autocorrelation functions (ACFs) of burst spectra show three features: a very narrow feature from radiometer noise, a narrow but resolved feature corresponding to the striations, and a broad feature related to the extent of the burst across the frequency band. The striation feature has a half width that varies from 2 to 5 MHz from burst to burst and is comparable to the scintillation bandwidth expected from the Milky Way in the direction of FRB 121102. The NE2001 electron density model provides an estimate s for the pulse broadening at 1 GHz. This predicts a scintillation bandwidth that ranges from 5 to 11 MHz across the 4.1 to 4.9 GHz band. We conclude that the measured ACFs and the NE2001 model prediction are consistent to within their uncertainties and that the narrow striations are due to Galactic scintillations.
A model for FRB 121102’s rotation measure (RM) and scattering measure (SM)
The measured RM implies a source frame value
We can use the previously estimated –270 pc cm-3 (in the source frame) and RMsrc to constrain the properties of the region in which the Faraday rotation occurs. In the absence of other information, we can set a constraint on the average magnetic field along the line of sight in the Faraday region with the ratio
If only a small portion of FRB 121102’s total DM is from the highly magnetized region, the field could be much higher.
The best constraint on pulse broadening comes from the measurement of the scintillation (diffraction) bandwidth of MHz at 4.5 GHz (see above). This implies a pulse broadening time at 1 GHz:
This scattering time is consistent with that expected from the Milky Way using the NE2001 model and therefore is an upper bound on any contribution from the host galaxy. Compared to scattering in the Milky Way, this upper bound is below the mean trend for any of the plausible values of DMHost, especially when the correction from spherical to plane waves is taken into account.
The ratio host-galaxy is a factor larger in the source frame but that is still far from sufficient to account for the apparent scattering deficit compared to the Galactic -DM relation. Given the apparent extreme conditions of the plasma in the host galaxy, it would not be surprising if its turbulence properties cause a scattering deficit. For example, scattering is reduced if the inner scale is comparable or larger than the Fresnel scale, either due to a large magnetic field or a high temperature.
Constraints on the properties of the Faraday region
Comparison of the magnetic field and thermal energy densities enables us to constrain the density (), electron temperature (), and length scale () of the region responsible for the observed Faraday rotation. We parametrize this relation with
where is a scaling factor, is the magnetic field strength, and is the Boltzmann constant. This assumes a 100% ionized gas of pure hydrogen with temperature equilibration between protons and electrons. Under equipartition, . In more densely magnetized regions, . Field reversals will reduce the total RM, requiring a lower value of in order to match constraints. The absence of free-free absorption at a frequency of 1 GHz sets an additional constraint on the permitted parameter space.
In Extended Data Fig. 6, we explore a range of physical environments. We consider a smaller lower limit, i.e. pc cm-3, on the dispersion measure than the previously estimated –270 pc cm-3, because not all of the DM may originate from the Faraday region. Galactic HII regions typically show rad m-2 and weak magnetic fields with , although calculations suggest it is possible for HII regions to achieve high RMs under some circumstances. Parameter space for typical HII region plasma at K is almost entirely excluded, and considering a range of possible HII regions sizes and densities shows that these are incompatible with the constraints. At higher , wide ranges of parameter space are permitted. In the case of equipartition, we have explicit unique solutions. For K, we find a density of on a length scale pc, i.e., comparable to the upper limit on the size of the persistent source. Higher temperature gas ( K) can be extended to pc. For both of these solutions, the characteristic magnetic field strength is 1 mG.
The large RM of FRB 121102 is similar to those seen toward massive black holes; notably, rad m-2 is measured toward Sgr A*, the Milky Way’s central black hole, and probes scales of Schwarzschild radii (0.001 pc). The constraints on , , and are also consistent with the environment around Sgr A* (Extended Data Fig. 6). The high RM toward the Galactic Centre magnetar PSR J17452900 (Fig. 3), rad m-2, at a projected distance of 0.1 pc from Sgr A* , is evidence for a dynamically organized magnetic field around Sgr A* that extends out to the magnetar’s distance. Notably, 4.5 years of radio monitoring of PSR J17452900 has shown a 5% decrease in the magnitude of the observed RM, while the DM remained constant at the 1% level (Desvignes et al., in prep.). This suggests large fluctuations in magnetic field strength in the Galactic Centre, on scales of roughly parsec.
While models considering the presence of only a massive black hole have been proposed, there is no observational precedent for microsecond bursts created in such environments. Rather, the FRB 121102 bursts themselves could arise from a neutron star, perhaps highly magnetized and rapidly spinning, near an accreting massive black hole. The proximity of PSR J17452900 to Sgr A* demonstrates that such a combination is possible. In this model, the black hole is responsible for the observed persistent source, whereas the bursts are created in the magnetosphere of the nearby neutron star.
Alternatively, the association of FRB 121102 with a persistent radio source has been used to argue that the radio bursts are produced by a young magnetar powering a luminous wind nebula. This model is not well motivated by Galactic examples, since the most luminous (non-magnetar powered) Galactic pulsar wind nebula is only times as luminous as the persistent source coincident with FRB 121102, and Galactic magnetars have no detectable persistent radio wind nebulae. Also, while giant flares from magnetars can produce relativistic outflows, an upper limit on the RM from one such outburst is 4 orders of magnitude below that observed for FRB 121102.
Nonetheless, under the millisecond magnetar model, the properties of the persistent source constrain the putative magnetar’s age to be between several years and several decades with a spin-down luminosity of to times higher than any local analog. Furthermore, the millisecond magnetar model predicts that the nebula magnetic field strength scales with the integrated spin-down luminosity of the magnetar. Extended Data Fig. 6 describes a range of sizes, densities, and temperatures for the Faraday-rotating medium that are consistent with Crab-like pulsar wind nebulae, known supernova remnants, and a simple model for swept-up supernova ejecta.
Data availability
The calibrated burst data are available, upon request, from the Corresponding Author.
Code availability
The code used to analyse the data is available at the following sites: PRESTO (https://github.com/scottransom/presto), PSRCHIVE (http://psrchive.sourceforge.net), DSPSR (http://dspsr.sourceforge.net).
References
Extended Data
Extended Data Table 1: List of 4.5 GHz Arecibo observations used in this study. These are a subset of all FRB 121102 observations to date.
Extended Data Table 2: Results of RM Synthesis and RMCLEAN. RMs were determined by fitting a quadratic function to the peak of the deconvolved Faraday spectrum. RM uncertainties were determined by dividing the nominal FWHM of the RM resolution element by twice the signal-to-noise ratio at the peak of the RM spectrum. RMdisp is the second moment (dispersion) of the RMCLEAN clean components discovered during the Faraday spectrum deconvolution. Upper limits indicate that the value scales with RM pixel size. A value of zero means that all clean components fell within the same pixel, and indicates a Faraday spectrum that is indistinguishable from being infinitely thin.
Extended Data Figure 1: Pulse profiles and spectra of the 16 Arecibo bursts. The bursts are de-dispersed to pc cm-3 (which minimizes the width of Burst #6) and plotted with time and frequency resolutions of s and MHz, respectively.
Extended Data Figure 2: Polarimetric properties of the 11 brightest bursts detected by Arecibo. a: linear polarization fraction of the bursts as a function of frequency. A solid line shows the theoretical depolarization due to intra-channel Faraday rotation calculated using Eqs. 3 and 4. b: PA∞ as a function of frequency. For both panels, values are averaged over consecutive channels. c: PA∞ as a function of time. A time offset is applied to each burst in order to show them consecutively. Vertical, dashed lines divide different observing sessions. All values in this figure have been corrected for the RM calculated with a global fit. Grey regions in b and c indicate the 1- uncertainty around the PA value from the global fit.
Extended Data Figure 3: Linear polarization fraction of the bursts as a function of RM. Different colours represent different observing sessions (see legend). A grey line indicates the average RM yielding the largest polarization fraction in the first observing session.
Extended Data Figure 5: RM and PA∞ values of the different bursts. Coloured, 1- error bars represent individual bursts, with central values highlighted by black dots. Horizontal grey regions are values obtained from a global fit. Values used in the figure are reported in Table 1.
Extended Data Figure 6: Physical constraints from source parameters. Parameter space for electron density () and length scale () of the Faraday region for three different temperature regimes, K. The shaded red region indicates parameter space excluded by optical depth considerations (). The solid black line gives the maximum DMHost permitted, while the shaded grey region shows the DM down to 1 . The solid blue line gives RMsrc. The shaded blue region gives the range . The intersection of grey and blue regions outside of the red region are physically permitted. The arrows indicates the upper limits on the sizes of the persistent source (left) and the star-forming region (right), respectively. The parallel dashed lines represent fits to a range of galactic and extragalactic HII regions. The parallel dotted lines represent the evolution of 1 and 10 M⊙ of ejecta evolving up to 1000 years at a velocity of in the blast-wave phase following a supernova. The filled downwards triangle and diamond are for the supernova remnants Cas A and SN 1987A, respectively. The filled circle represents the mean density and diameter of the Crab Nebula, whereas the filled square represents the characteristic density and length scale of a dense filament in the Crab Nebula. The star indicates the density at the Bondi radius of Sgr A*.