Sub-second periodicity in a fast radio burst

Fast radio bursts (FRBs) are millisecond-duration flashes of radio waves that are visible at distances of billions of light-years. The nature of their progenitors and their emission mechanism remain open astrophysical questions. Here we report the detection of the multi-component FRB 20191221A and the identification of a periodic separation of 216.8(1) ms between its components with a significance of 6.5 sigmas. The long (~3 s) duration and nine or more components forming the pulse profile make this source an outlier in the FRB population. Such short periodicity provides strong evidence for a neutron-star origin of the event. Moreover, our detection favours emission arising from the neutron-star magnetosphere, as opposed to emission regions located further away from the star, as predicted by some models.

Fast radio bursts (FRBs) are millisecond-duration flashes of radio waves that are visible at distances of billions of light-years. 1 The nature of their progenitors and their emission mechanism remain open astrophysical questions. 2 Here we report the detection of the multi-component FRB 20191221A and the identification of a periodic separation of 216.8(1) ms between its components with a significance of 6.5σ. The long (∼ 3 s) duration and nine or more components forming the pulse profile make this source an outlier in the FRB population. Such short periodicity provides strong evidence for a neutron-star origin of the event. Moreover, our detection favours emission arising from the neutron-star magnetosphere, 3,4 as opposed to emission regions located further away from the star, as predicted by some models. 5 Operating on the Canadian Hydrogen Intensity Mapping Experiment (CHIME), CHIME/FRB 6 is an ongoing experiment to find and study a large number 7 of FRBs. CHIME is a cylindrical North-South oriented transit radio interferometer observing in the 400-800-MHz range. Upon detection of an FRB, the so-called intensity data, i.e. the total intensity of the signal as a function of time and frequency, are stored. Additionally, channelized complex voltages (referred to as baseband data) with full polarisation information are stored for a subset of FRBs (see Methods).
In fewer than 0.5% of the events detected by CHIME/FRB, five or more separate components are visible in the pulse profiles 7 obtained by summing all frequency channels of an intensity dataset after correcting for the effects of the dispersion measure (DM). Particularly striking is FRB 20191221A, with a total duration of roughly three seconds and at least nine overlapping components ( Figure 1 and Table 1). No other FRB candidate observed by CHIME/FRB contains a comparable or greater number of sub-components. Since the detection pipeline of CHIME/FRB is not optimized to find bursts longer than 128 ms, it is possible that some events with comparable duration eluded detection. However, FRB 20191221A was identified using the individual peaks in its profile. Therefore, we expect the fraction of missed long-duration events to be small if single peaks can be identified in their pulse profiles. A detailed analysis of the  Figure 1: Radio signal from FRB 20191221A. a, Waterfall plot of the signal intensity (colourcoded) as a function of time and frequency. Frequency channels missing or masked due to radio frequency interference are replaced with off-burst median values and are indicated in red. Effects of dispersion have been removed, and data have been averaged to 3.125 MHz frequency resolution and to 7.86432 ms time resolution. b, In black, the pulse profile obtained by averaging the frequency channels of the waterfall plot where signal is visible. The Gaussian function convolved with an exponential used to model nine peaks in the profile is plotted in red, and peak locations are highlighted by vertical lines. c, Residuals of the fit, with a red horizontal line at zero residual.  Table 1: Properties of FRB 20191221A. Uncertainties are reported at 1-σ confidence level. The arrival time is that of the brightest sub-burst at the Solar System's barycentre and infinite frequency. The DM is calculated to maximize the peak S/N in the timeseries. The scattering timescale is referenced to the centre of the band, i.e. ∼ 600 MHz. Fluence is for the full band-averaged profile, and peak flux is the maximum in the profile (with 1-ms time resolution). completeness of the CHIME/FRB search pipeline has been presented elsewhere. 7 Significant peaks are visible in the power spectrum of FRB 20191221A obtained by performing a Fast Fourier Transform (FFT) on its pulse profile (Figure 2), indicating a possible periodicity in the times of arrival (ToAs) of single components of the pulse profile. To confirm this, individual sub-components have been fitted with a Gaussian function convolved with a single exponential to account for scattering, a pulse broadening caused by the propagation of the radio waves in turbulent plasma. 8 The resulting ToAs have been used to perform a timing analysis around the initial period derived from the power spectra. A simple timing model with two free parameters, a period and an arbitrary overall phase, have been fitted to the data. The b, Residuals of a timing analysis assuming that the peaks forming the FRB profile are separated by integer numbers times a period P = 216.8 ms. 1-σ error bars are often smaller than the symbol size. The horizontal dotted line indicates a phase of zero around which residuals have been rotated. c, Study of the statistical significance of the measured periodicity by using a periodicity-sensitive scoreŜ described in the Methods. The grey histogram has been obtained with an ensemble of simulations, whereas the value measured for the FRB is represented with a vertical dotted line. The corresponding probability of obtaining such a periodicity by chance is indicated on the plot.
fit residuals are shown in Figure 2 and the period is reported in Table 1. A comparison of the   residuals with a straight line, expected for a perfectly  Millisecond to second periodicities may suggest that the bursts are generated by Galactic radio pulsars that have been misidentified as extragalactic. However, FRB 20191221A has a measured DM ∼ 4 times larger than the maximum value expected by models of the Milky Way electron content. 9, 10 We searched for evidence of ionized 11,12 or star-forming 13 regions in the direction of the FRB that could account for the excess DM but found none. We conclude that the source is likely extragalactic.
No additional bursts have been detected from FRB 20191221A up to March 10th, 2021 above an S/N of 9 at a position consistent within ∆RA = 2.2 deg cos −1 (Dec) and ∆Dec = 1 deg of these three sources using an algorithm built on a density-based spatial clustering of applications with noise (DBSCAN). 14 The nominal DM range threshold for the clustering was set to 13 pc cm −3 , corresponding to the largest DM uncertainty in the real-time pipeline of CHIME/FRB. 6 CHIME/FRB continues to monitor the sky location of these FRBs daily for possible additional bursts.
Multi-component bursts that are associated with repeating sources of FRBs often exhibit downward-drifting subbursts, 15 however, all of the components forming FRB 20191221A show a similar spectrum. This is visible in ∼ 5% of the FRB population detected by CHIME/FRB. 16 It must be noted that the spectrum of FRB 20191221A is affected by the telescope response, mainly due to the detection at a location offset from the centre of a formed beam, which produces strong bandpass effects (see Methods). 17 Our modeling of the pulse profile of FRB 20191221A, visible in Fig. 1 and described in detail in the Methods, shows that its single components have a relatively narrow width of 4(1) ms on average, even though they overlap due to an extreme scattering timescale τ s = 340(10) ms at 600 MHz. Although this estimate may be affected by unresolved features in the profile mimicking an exponential decay, it is clear that the FRB emission experienced strong scattering, significantly in excess of the Galactic contribution expected given its sky position, 9 pointing to propagation through a turbulent plasma. The pulse width corrected for the pulse broadening corresponds to a duty cycle of ∼ 1.8%, consistent with Galactic radio pulsars. 18 It is worth noting that all of the emission from this FRB is consistent with single components overlapping due to the large scattering and no envelope of emission is required in our fit, whose residuals are consistent with noise, as is visible in Fig. 1.
Leading theories for the origin of FRBs are related to magnetars 2,19,20 and are divided into models where the emission is either generated in the star's magnetosphere or triggered in plasma regions by a flare of the star. The detection of periodicity is naturally explained by the first class of magnetar models, and it has been extensively observed in Galactic neutron stars, albeit with orders of magnitude lower luminosities. 18,21 By contrast, the second class of models does not necessarily predict a millisecond modulation in the emitted signal. 5 The periodic structures in the bursts could be explained by a rotating neutron star with beamed emission similar to Galactic radio pulsars where, for an unknown reason, a train of single pulses has an abnormally high luminosity for a short period of time. One unlikely scenario to explain the larger radio luminosity of the FRB compared to Galactic sources is that it may represent the observation of a gravitationally micro-lensed extragalactic pulsar. Alternatively, the single components of periodic FRBs could be generated by the magnetospheric interaction of merging neutron stars. These possibilities are explored in the Methods, along with suggested tests of these models. In the meantime, CHIME/FRB is continuing to detect hundreds of FRBs, which should allow for additional periodicities to be detected in the near future.

CHIME/FRB sensitivity beams
The CHIME radio telescope is a transit interferometer formed by 1024 dual-polarisation antennas observing in the 400-800 MHz range. The field of view of single antennas summed together incoherently is defined as the primary beam of the telescope, whose FWHM sensitivity spans ∼ 110 degrees in the N-S direction and 1.3 (2.5) degrees the E-W direction at the top (bottom) of the observing bandwidth. 6 The telescope antennas can be added coherently to produce a formed (or synthesized) beam in one or more directions and increase the telescope sensitivity towards those directions. Formed beams have a size between ∼ 0.3 and 0.7 degrees, depending on the observing frequency and zenith angle. 6 Therefore, their sensitivity and bandpass vary spatially more rapidly than those of the primary beam.
In the real-time search for FRBs, 1024 beams are formed on the sky within the primary beam via an FFT (FFT beams). 27 The total intensity measured by these FFT beams as a function of time and frequency is referred to as intensity data. Intensity data have a time resolution of 0.98304 ms and are divided into 16, 384 channels. Additionally, for a part of the detected FRBs, ∼ 100 ms of baseband complex voltages are also stored. 6 The baseband data have a time resolution of 2.56 µs, are divided into 1024 frequency channels and contain full polarisation information. Thanks to the phase information available in baseband data, synthesized beams can be formed to virtually any position on the sky within the telescope's field of view. 28

Sources of interest
We initially identified FRB 20191221A as an interesting source due to the number of peaks in its pulse profile. Only intensity data was stored for this source. We then visually inspected the sample of events detected by CHIME/FRB for other similar sources. We found none with comparable characteristics and the closest for number of components were found to be FRBs Extended Data Figure Table 1. c, d, Residuals of a timing analysis assuming that the peaks forming the FRB profile are separated by integer numbers times these periods, respectively. 1-σ error bars are often smaller than the symbol sizes. Horizontal pink lines indicate a phase of zero around which residuals have been rotated. e, f, Study of the statistical significance of the measured periodicity by using the periodicitysensitive scoreŜ. The grey histograms have been obtained with an ensemble of simulations, whereas the value measured for each FRB is represented with a vertical pink line. The corresponding probability of obtaining such a periodicity by chance is indicated on the plots. Both intensity and baseband data are available for FRBs 20210206A and 20210213A. This allowed us to localize the sources with sufficient precision to form beams in these directions and, therefore, limit the effect of their bandpass. 28 Therefore, the lack of emission above ∼ 500 MHz seen for FRB 20210213A is astrophysical. The DM of FRB 20210213A is ∼ 10 times the expected Galactic contribution corresponding to an estimated redshift of z ∼ 0.4. 34 FRB 20210206A, located at a Galactic latitude b = −2.5 deg, has a DM ∼ 1.5 times larger than the Galactic contribution, placing it at a supposed redshift z ∼ 0.1 but with the caveat that its extragalactic nature is less certain given the larger uncertainty of models at low latitudes. 35 As opposed to FRB 20191221A, FRBs 20210206A and 20210213A show narrower widths and shorter scattering timescales (see Extended Data Table 1) typical of the FRB population. 16 The baseband data also allowed us to study the polarisation properties of FRBs 20210206A and 20210213A. A rotation measure RM = +193.6(1) rad m −2 has been mea-  (4), 496.7(2) * For circumpolar sources (δ > +70 • ), the two entries correspond to exposure in the upper and lower transit, respectively.
Extended Data Table 1: Properties of FRBs 20210206A and 20210213A. All quantities have been calculated with baseband data except for fluences and fluxes. Uncertainties are reported at 1-σ confidence level. The arrival time is that of the brightest sub-burst at the Solar System's barycentre and infinite frequency. The DM is calculated to maximize the peak S/N in the timeseries. The scattering timescale is referenced to the centre of the band, i.e. ∼ 600 MHz. Fluence is for the full band-averaged profile, and peak flux is the maximum in the profile (with 1-ms time resolution). sured for FRB 20210206A, suggesting a significant extragalactic contribution. On the other hand, FRB 20210213A appears to be unpolarised. Possible reasons for this are discussed in the following, together with a detailed description of the polarisation analysis.

Properties of the bursts
To localize an FRB with CHIME/FRB intensity data, we fit the spectra of the burst detected in different FFT beams with a model of the CHIME/FRB beams and an underlying burst spectrum using a Markov chain Monte Carlo (MCMC) method. 36 The model of the CHIME/FRB beams contains a description of both the synthesized 27,37 and primary 7 beams and the underlying spectrum is modeled as a Gaussian. The free parameters are therefore width, mean, and amplitude of the underlying Gaussian model spectrum, along with the sky position. We use a flat prior on the position of the event that spans 5 • to either side of meridian in E-W, while in N-S the prior spans the extent of the beams that detected the event. The position and uncertainties are derived from the 2D posterior distribution (in 'x', the E-W coordinate, and 'y', the N-S coordinate) marginalized over the parameters of the Gaussian spectral model. Since FRB 20191221A did not have baseband data available, the position reported in Table 1 is derived from the intensity localization described in the present section. The posterior probability distribution is double-peaked in the E-W direction and so two positions are reported for the event.
A detailed description of the algorithm used to obtain the sky position of FRBs by using baseband data has been presented elsewhere. 28 In summary, a grid of partially overlapping beams is produced around the intensity localization and a total S/N value is calculated in each beam. The resulting intensity map of the signal is fitted with a mathematical model describing the telescope response to determine the source position. The localization and its uncertainties have been calibrated with a sample of sources with a known position to account for any unmodelled systematics. 28 Flux and fluence calculations are determined using the intensity data for each burst with the same method presented in previous CHIME/FRB papers. 17,32,[38][39][40][41] In summary, meridian transits of steady sources with known spectra are used to sample the conversion between beamformer units and Janskys as a function of frequency across the N-S extent of the primary beam. For each burst, the beamformer to Jansky conversion closest in zenith angle (assuming N-S symmetry) is applied to the intensity data to obtain a dynamic spectrum in physical units roughly corrected for This is what we report in Table 1 for FRB 20191221A. However, for bursts that have a baseband localization, we can achieve more realistic fluence results by using the beam model to scale between the location of the calibrator at the time of transit and the location of each FRB. This is what we report in Extended Data Table 1 for FRBs 20210206A and 20210213A.
The exposure of the CHIME/FRB system to the sources reported in this work was determined for the interval from August 28, 2018 to March 1, 2021. For each source, the exposure is calculated by summing the duration of daily transits across the FWHM region of the synthesized beams at 600 MHz. Two of the three sources, having declinations > +70 • , transit through the primary beam twice per day. These sources have their upper and lower transit exposures calculated separately as the beam response for the two transits is different. While calculating the total exposure, we include daily transits for which the CHIME/FRB detection pipeline was fully operational, which is determined using recorded system metrics. Transits occurring on days when the detection pipeline was being tested or upgraded are not included. Additionally, system sensitivity varies on a day-to-day basis due to daily gain calibration as well as changes in the detection pipeline and the RFI environment. We evaluate the variation in sensitivity for each sidereal day using observations of 120 Galactic pulsars with the CHIME/FRB system. 39,41 For each source, transits for which the sensitivity varied by more than 10% from the median in the aforementioned observing interval are excised from the total exposure. On average, the observing time corresponding to the excised transits amounted to 4% of the exposure for each source. The uncertainty in the source declination is a source of error in the measurement of the exposure. The source declination dictates where the transit path cuts across a synthesized beam with the transit duration being maximum if the path crosses the beam centre and zero if the path lies between two beams. 40 We estimate the resulting uncertainty in the exposure by generating a uniform grid of positions within the 68% confidence localization region for each source. The reported exposure in Table 1 and Extended Data Table 1 is the average for these sky positions with the error corresponding to the standard deviation.

Modeling of pulse profiles
FRB 20191221A is composed of multiple peaks overlapping due to a large broadening caused by scattering. To properly calculate the periodic separation among its components, it is important to avoid human bias in selecting significant peaks in the pulse profile. We reduced human bias in the following way. First, we smoothed the pulse profile using the Savitzky-Golay filter as implemented in SciPy. 1 The filter requires two input parameters, the window length, and the polynomial order. We explored the space of the two parameters up to a window of 600 ms and a polynomial order of 12. For each combination, separate peaks in the smoothed pulse profiles were identified from local maxima if the peaks were wider than 3 bins to avoid noise spikes.
We grouped together different combinations of parameters yielding the same number of peaks, noting their location in the profile and excluding values lower than 7 and higher than 14 peaks as visually implausible. For each value of the peak number between 7 and 14, the different peak locations resulting from each of the initial parameter combinations were grouped with a kernel density estimator, varying its window until we obtain the expected number of peaks. We apply this method to obtain seven combinations of initial peak positions composed of 7 to 14 1 https://scipy.org/ Extended Data Figure 3: Reduced chi-square test as a function of the number of components used to model the profile of FRB 20191221A. The vertical line highlights the chosen number of components, while the horizontal line is placed at the χ 2 red value for 9 components. The minimum χ 2 red variation that can be measured confidently with our data is estimated with Eq. 2 and plotted as en error bar for each number of peaks. peaks. We model each peak using the exponentially modified Gaussian (EMG) described by 7,42 where erfc is the complementary error function, A is the signal amplitude, µ and w are the Gaussian mean and width, respectively, and τ is the scattering timescale. For each number of peaks, the pulse profile is modelled with a sum of one EMG function per peak through an MCMC sampling using emcee 36 with wide uniform priors for all parameters. The scattering timescale is set to be the same for all peaks, while the other parameters are allowed to vary.
The resulting χ 2 red values, presented in Extended Data Figure 3, are used to select the number of peaks yielding the best model. We estimate the minimum variation in the χ 2 red value that we can  (5)  7 1952(1) · · · · · · 8 2171(2) · · · · · · 9 2604(2) · · · · · · Extended Data Table 2: List of times of arrival (ToAs) for the peaks forming each event. The ToA of each component is reported in milliseconds relative to the first peak. 1-σ uncertainties on the last digit are indicated in parenthesis.
measure confidently given the number of free parameters N as As visible in Extended Data Figure 3, the model with 9 components has the smallest χ 2 red that deviates significantly from the previous values. Therefore, we choose this as the best model to reproduce the data and use its parameters in our analysis. Choosing a different number of peaks in the profile leads to values of the significance of the periodicity that are still very high, especially for a number of peaks between 9 and 12.
For FRBs 20210206A and 20210213A, the components forming these two events do not overlap. Therefore, it was not necessary to perform the method described above. Instead, we directly used the locations of the peaks in the smoothed profiles as initial conditions for the MCMC sampling using the same model described above. The scattering timescale and average width of the resulting models are presented in Table 1 and Extended Data  Table 3: Statistical significance of the FRB periodicities. For each FRB, we report the number of ToAs measured from the profile, the number of gaps and trials considered in theŜ-periodicity analysis, and the resulting probability and significance. The values in the last four columns are derived using the Rayleigh (Z 2 1 ) test and show the period obtained in the analysis, the resulting value of the test, and the false alarm probability and significance of the periodicities.

Timing analysis
The periodicities of the three FRBs were initially investigated through a power spectrum of the pulse profile. This has been computed as the absolute square of the Fast Fourier Transform (FFT) of the pulse profiles. The FFT has been calculated with the implementation offered by the SciPy module. The resulting power spectra presented in Fig. 2 and Extended Data Fig. 2 show clear peaks for the three sources. The periods corresponding to the most prominent peaks in the power spectra have been refined through an additional timing analysis. We ran a leastsquares fit as implemented in SciPy using a simple model with the period and an arbitrary overall phase as the only free parameters, calculating residuals as the modulo of ToAs (reported in Extended Data Table 2) and the period. The results of these fits are the values reported in Table 1 and Extended Data Table 1.

Significance of the periodicity:Ŝ score
The following steps were used to estimate the significance of the periodicity calculated for each of the three sources and reported in Table 1 and Extended Data Table 1. 1. For each event, we use the ToAs reported in Extended Data Table 2 to compute a statistiĉ S which is sensitive to periodicity. The construction ofŜ is described below.
2. To assign a statistical significance, we use a frequentist approach. We evaluate the statistic Extended Data Figure 4: Times of arrival (ToAs) of the components of FRB 20191221A as a function of their measured cycle. The cycle is defined in Eq. 3. The periodicity appears clearly as the points fall nearly on the straight gray line, which highlights the trend expected for a period of 216.8 ms. Vertical lines mark gaps where no pulse is observed within one period.
S on simulated events and rank the 'data' valueŜ data within the ensemble of simulated valuesŜ sim , obtaining a p-value. The results are summarized in Extended Data Table 3. 3. For FRB 20191221A, there is an extra step. With 10 8 simulations, we find that none of theŜ sim values exceedŜ data . This shows directly that the level of periodicity in this 9-component event is extraordinarily unlikely to occur by random chance. To assign a p-value, we fit an analytic model PDF to the tail of theŜ sim distribution and integrate the model PDF.
In the rest of this section, we describe the details in the above steps.
To motivate the definitions which follow, consider Extended Data Figure 4. Each point is one pulse in the 9-component event, with the ToA t i on the y-axis. The values on the x-axis are the 9-element integer-valued vector n i = (0, 2, 3,5,7,8,9,10,12) (3) Figure 4 as the points falling nearly on a straight line. Note that there are four 'gaps' in this example, i.e. pulse periods where no pulse is observed (either because it is physically absent or buried in the noise). Each gap is represented by consecutive entries in n i which differ by 2 (rather than 1). Based on this picture, we construct statistics as follows. For a fixed choice of gap vector n i , we fit the points (n i , t i ) in Extended Data Figure 4 to a straight line

Periodicity appears in Extended Data
where this equation defines the residuals r i , and the parameters (f , T 0 ) have been chosen to minimize i r 2 i . Note that this fitting procedure weights all ToAs equally, and does not use statistical errors on ToAs. We define a statistiĉ The maximum in the equation is taken over all trial gap vectors n i with ≤ G gaps, where G is an input parameter to the pipeline. More formally, we take the maximum over integer-valued vectors n i such that n i = 0, n i−1 < n i , and n p < n 1 + p + G, where p = len(n) is the number of peaks. The number of such vectors is By default, we choose G = N gaps + 1, where N gaps is the number of gaps that are empirically seen in each event. This is a conservative choice, which deliberately ensures that N trials is a few times larger than the value obtained with G = N gaps (see Extended Data Table 3).
Our procedure for simulating ToAs has two parameters: a mean spacingd, and a dimensionless 'exclusion' parameter 0 ≤ χ ≤ 1 defined by χ =

Minimum allowed spacing between pulses
Mean spacingd between pulses .
We simulate a sequence t 1 < · · · < t p of ToAs by independently randomly generating arrival We also tried an exponential distribution, but find that it gives lower p-values (higher significance). To be conservative, we use the uniform distribution throughout. When we assign a statistical significance by rankingŜ data within a histogram ofŜ sim values, we find that the value ofd does not affect the statistical significance, while the statistical significance decreases as χ is increased in the simulations. However, the statistical significance is not strongly dependent on χ for a reasonable range of the parameter. Therefore, we choose χ = 0.2 to represent our analysis.
For FRB 20191221A, there is an extra step. With 10 8 simulations, we find that none of thê S sim values exceedŜ data . Therefore, we fit the 10 5 largest values S i of theŜ sim -distribution to an analytic PDF of the form by maximizing the likelihood i P (S i |a, q). The analytic PDF p(S) is an excellent visual fit to the tail of theŜ sim distribution. More quantitatively, a KS test shows no statistical difference between p(S) and the simulations (p-value 0.58). To assign a bottom-line p-value, we integrate the analytic PDF from S =Ŝ data to S = ∞. This gives a p-value of 6.7 × 10 −11 , corresponding to Gaussian significance 6.5σ.
The preceding derivation of the statisticŜ was heuristic, based on intuition from Extended Data Figure 4. However,Ŝ can also be interpreted as a likelihood ratio statistic. This provides a systematic derivation, and also shows thatŜ is near-optimal. Let H 0 be the null hypothesis that the ToAs t i are Gaussian distributed with mean T and variance σ 2 . In this model, the conditional likelihood of obtaining ToAs t i given model parameters (T , σ) is Let H 1 be the alternate hypothesis that the t i are given by the linear regression in Eq. (4), where the residuals r i are Gaussian with variance σ 2 . The conditional likelihood of obtaining t i given Then a short calculation shows that theŜ statistic is the log-likelihood ratio of the two models, after maximizing all model parameterŝ As a further test of our pipeline, we applied theŜ statistic to CHIME/Pulsar 43 observations of the bright pulsar PSR B1919+21. If we select three pulses with one gap, periodicity is detected at the ∼ 3σ level. If we select four pulses or more with either one or two gaps, periodicity is detected at > 4σ.

Rayleigh (Z 2 1 ) periodicity analysis.
We have shown in the previous section that theŜ statistic is nearly optimal for determining the statistical significance of the periodicity observed in each FRB. However, we also calculated the significance through a different statistic that is commonly used in studies of periodicities of high-energy pulsars, the Rayleigh (Z 2 1 ) statistic. 44,45 In general, the Z 2 n test statistic is defined as where N is the number of peaks, n is the number of harmonics, and φ j is the phase of each ToA, t j . The phase of each ToA is determined using φ j = νt j , where ν is the modulation frequency. In the following analysis, we use n = 1 harmonics in Eq. (14), which corresponds to the Rayleigh test, and the ToAs listed in Extended Data Table 2. All of the ToAs were weighted equally in our analysis.
We performed a blind search for periodicity using the Z 2 1 test statistic defined in Eq. (14) and the ToAs listed in Extended Data Table 2 Table 3.
where ν 0 is the putative periodicity determined from the ToAs measured from each event.
In this analysis, the tail-fitting procedure described above is not used. Instead, the FAPs are calculated based on the number of Monte Carlo simulations that have a maximum Z 2 1 test statistic which exceeds the Z 2 1 value obtained using the measured ToAs. We find that the periodicity observed from FRB 20191221A has a significance of 6.2σ using this method. The equivalent significance of the periodicities observed from FRBs 20210206A and 20210213A are both < 1σ using the Rayleigh (Z 2 1 ) statistic (see Extended Data Table 3).

Polarisation analysis
Full polarisation information is stored for FRBs 20210206A and 20210213A. The polarisation analysis follows a similar procedure to that previously applied to other CHIME-detected FRBs. 41,46 In particular, an initial RM estimate is made by applying RM-synthesis 47,48 to the Extended Data Figure  For FRB 20210206A, RM-synthesis results in an unambiguous RM detection. This detection is refined by applying a Stokes QU-fitting routine that directly fits for the modulation between Stokes Q and U from Faraday rotation as well as the modulation between U and V introduced by an instrumental delay between the two linear polarisations. Optimal values are determined numerically through Nested Sampling, a Monte Carlo method that seeks to optimize the likelihood function given a model and data. Further details on the CHIME/FRB polarisation analysis pipeline are presented elsewhere. 49 The instrumental delay, once fitted for, can be used to produce a delay-corrected spectrum. The RMs produced from these two independent methods agree within the measurement uncertainties. An RM = +193.6(1) rad m −2 is given by QU-fitting and it is used to produce the polarized burst profile shown in Extended Data Figure 5.
We ues beyond this range, bandwidth depolarisation from Faraday rotation within a single frequency channel becomes significant at the native channelization of CHIME/FRB baseband data. 49 We developed an algorithm 49 that employs a phase-coherent method of correcting for bandwidth depolarisation in data for which the electric field phase is retained. 51 52 We note that ionospheric RM has not been corrected for in our analysis, but it does not exceed a few rad m −2 .

Model of gravitational lensing
As a possible explanation to the magnification that would be necessary to observe a radio pulsar located in another galaxy, we explore the observability of a pulsar that is gravitationally microlensed. In this model, pulses from a radio pulsar in a binary system are magnified in intensity by the gravitational field of its companion. The magnification would be modulated in time such that the signal from the pulsar would be convolved with a bell-shaped curve. The pulse profiles detected for the three FRBs presented here are qualitatively consistent with this morphology.
We consider a binary system with a pulsar of mass M s = 1 M and explore the parameter space of lensing masses, system alignment, and orbital separations to explain the properties of the three FRBs.
We use a test pulsar at a distance of 1 Gpc (approximately the DM-inferred distance of FRB 20210213A) emitting periodic pulses with a luminosity of 1 Jy kpc 2 , comparable to Galactic radio pulsars. Without any magnification, such a pulsar would be observed on Earth with a peak flux of ∼ 10 −11 Jy. Therefore, a magnification µ 10 11 is needed to explain the fluxes of ∼ 1 Jy observed for the three FRBs here reported.
We constrain the parameter space by first requiring the lens mass M lens to be able to magnify the pulsar 53 where G is the gravitational constant and f is the observed frequency.
The orbital separation of the binary system D ls and its alignment are constrained by requiring that we only observe the magnification curve for the duration of the event ∆t. The minimum alignment angle β of the binary system can be calculated as 53 Extended Data Figure 6: Parameters of a binary system producing a radio signal compatible with the FRBs presented here through gravitational lensing. The system, located at 1 Gpc, contains a 1 M pulsar emitting 1 Jy pulses which are lensed by its companion. The allowed parameter space is shown with brighter colors as a function of the minimum alignment angle (color-coded), companion mass and separation of the binary system.
where the Einstein angle θ E is defined as D os is the distance between the observer and the source, D ol is the distance between the observer and the lens, and D ls is the orbital separation that we want to constrain. β is further constrained by β ω∆t, where ω is the orbital angular velocity of the lens, defined as and ∆t is the duration of time in which the pulsar is magnified enough to be detected above the noise floor. We use Eqs. 16,17, and 19 to constrain the properties of the binary system and plot the parameter space shown in Extended Data Figure 6.
The orbital inclination of the lensing system would need to be aligned to the line of sight within 10 −17 − 10 −18 arcseconds, depending on the lens mass. The small required alignment angle implies a low probability to observe such an event despite the large trial factor provided by the number of galaxies within a Gpc distance. In addition, in this scenario, we would expect to detect a larger number of FRBs from gravitationally lensed pulsars in the nearby universe since they would require a lower magnification. Therefore, the absence of multi-peaked, periodic FRBs with a small DM excess 7 implies that it is unlikely that we have detected a gravitationally micro-lensed pulsar of average luminosity at a distance of 1 Gpc.

Model of merging compact objects
A different model that could produce periodic FRBs considers merging neutron stars that emit the FRB signal. One possible interaction of merging neutron stars to produce periodic FRBs is through a unipolar inductor process where the companion orbiting through the magnetic field acts as a conductor driving a current loop. The latter accelerates electrons and positrons to emit curvature radiation [54][55][56][57] in orbital frequencies ranging from few Hz to kHz, corresponding to orbital separations of 10-1000 km in the binary neutron star case. Another proposed mechanism to extract energy is through the magnetic braking and spin-orbital synchronization of merging binary neutron stars. 58,59 The coherent emission is hypothesized to arise from the magnetosphere in a manner roughly similar to isolated pulsars as the rotation rate of one of the neutron stars rapidly increases or decreases to synchronize with the binary rotation. In such a case, FRB emission may show multiple peaks corresponding to a favourable orbital phase for a range of orbital periods.
The loss of angular momentum and energy through gravitational wave radiation causes the compact object binary orbits to decay with a predictable relation between the orbital angular frequency ω and time t. We consider the equation 60 Here, m = m 1 + m 2 and η = m 1 m 2 /m 2 are the total mass and reduced mass, respectively, of two components of mass m 1 and m 2 . We have assumed that spin-orbit and spin-spin coupling are negligible. Since the orbital angular frequency is related to the observed pulse period by it follows thatṖ For a given set of trial masses, (m 1 , m 2 ), we integrated numerically Eq. (21) twice to calculate the orbital phase φ, i.e. we go fromω(ω, t) → ω(t) → φ(t). The initial orbital period was chosen to be larger than the periods measured for these FRBs and the system was evolved until ω = 2π × 10 6 rad s −1 , very close to the final merger. We numerically inverted φ(t) to get t(φ), the time of passage of the components through a specific orbital phase. We use this to fit the ToAs through the same procedure used to study the significance of the periodicity. We modified Eq. (4) as where φ 0 is an arbitrary initial phase and the integer-valued n i vector is defined in Eq. An eventual future detection of repeating bursts from these FRB sources would disprove this model.