Propagation effects at low frequencies seen in the LOFAR long-term monitoring of the periodically active FRB 20180916B

LOFAR (LOw Frequency ARray) has previously detected bursts from the periodically active, repeating fast radio burst (FRB) source FRB 20180916B down to unprecedentedly low radio frequencies of 110 MHz. Here we present 11 new bursts in 223 more hours of continued monitoring of FRB 20180916B in the 110-188 MHz band with LOFAR. We place new constraints on the source's activity window $w = 4.3^{+0.7}_{-0.2}$ day, and phase centre $\phi_{\mathrm{c}}^{\mathrm{LOFAR}} = 0.67^{+0.03}_{-0.02}$ in its 16.33-day activity cycle, strengthening the evidence for its frequency-dependent activity cycle. Propagation effects like Faraday rotation and scattering are especially pronounced at low frequencies and constrain properties of FRB 20180916B's local environment. We track variations in scattering and time-frequency drift rates, and find no evidence for trends in time or activity phase. Faraday rotation measure (RM) variations seen between June 2021 and August 2022 show a fractional change $>$50% with hints of flattening of the gradient of the previously reported secular trend seen at 600 MHz. The frequency-dependent window of activity at LOFAR appears stable despite the significant changes in RM, leading us to deduce that these two effects have different causes. Depolarization of and within individual bursts towards lower radio frequencies is quantified using LOFAR's large fractional bandwidth, with some bursts showing no detectable polarization. However, the degree of depolarization seems uncorrelated to the scattering timescales, allowing us to evaluate different depolarization models. We discuss these results in the context of models that invoke rotation, precession, or binary orbital motion to explain the periodic activity of FRB 20180916B.


INTRODUCTION
Fast radio bursts (FRBs; see Petroff et al. 2022 for a review) are astrophysical transients that produce millisecond-duration coherent radio emission.Their characteristic large dispersive delays are the result of their extragalactic distances, which have been confirmed in some cases by a robust association with a host galaxy (e.g., Chatterjee et al. 2017;Bannister et al. 2019;Marcote et al. 2020;Kirsten et al. 2022;Niu et al. 2022;Ravi et al. 2022).They are many orders-of-magnitude more luminous than the individual pulses from Galactic pulsars, including even the giant pulses from the young Crab pulsar.Since the serendipitous discovery of the 'Lorimer burst' (FRB 20010724A;Lorimer et al. 2007), over 600 FRB sources have been catalogued, most by the Canadian Hydrogen Intensity Mapping Experiment's FRB detection system (CHIME/FRB; CHIME/FRB Collaboration et al. 2021Collaboration et al. , 2023)).In this sample, most have been one-off events (apparent non-repeating FRBs); ≈ 4% of FRBs have been seen to repeat (Spitler et al. 2016;Pleunis et al. 2021b; CHIME/FRB Collabora-★ E-mail: a.gopinath@uva.nltion et al. 2023) -FRB 20180916B with a periodic activity cycle of 16.33 ± 0.12 day (CHIME/FRB Collaboration et al. 2020;Pleunis et al. 2021b), and FRB 20121102A with a possible ≈ 160-day periodicity (Rajwade et al. 2020;Cruces et al. 2021).It has been argued that the apparent one-off events are less active repeating sources that have simply not been observed to repeat yet (Caleb et al. 2019;Ravi 2019;James et al. 2020;CHIME/FRB Collaboration et al. 2023).However, burst morphological and spectral differences between repeating and apparently non-repeating FRBs point to the possible existence of at least two distinct populations of FRBs (Hessels et al. 2019;Pleunis et al. 2021b;CHIME/FRB Collaboration et al. 2023).Alternatively, this dichotomy may indicate a different burst mechanism in the same type of progenitor source (Kirsten et al. 2023), or may be a consequence of beaming direction (at least in the case of burst duration; Connor et al. 2020).
We can further explore the origins of FRBs through studies of their host galaxy environments, changing burst activity (Li et al. 2021a;Hewitt et al. 2022;Nimmo et al. 2023), and evolution of burst properties such as dispersion measure (DM) and Faraday rotation measure (RM; Michilli et al. 2018b;Hilmarsson et al. 2021b).This can help establish differences, if any, in the hosts (Bhandari et al. 2022) and local environments of FRBs (Kirsten et al. 2022).
FRB 20180916B was localized to the vicinity of a star-forming region in a spiral host galaxy (Marcote et al. 2020;Tendulkar et al. 2021).At a luminosity distance of ∼ 149 Mpc ( = 0.0337 ± 0.0002), it is offset by about 250 pc from the peak brightness of the nearest star-forming knot in its host galaxy.Assuming that it was born in this region and has subsequently moved away from it, the inferred age of ∼ 10 5−7 yr is inconsistent with that of a young magnetar, and instead led Tendulkar et al. (2021) to suggest a high-mass X-ray or gammaray binary progenitor model (though an OB runaway scenario cannot be ruled out).The firmly established activity period (CHIME/FRB Collaboration et al. 2020) aids multi-frequency follow-up observations of the source because it provides some predictability to use in the planning of observations.As a result, FRB 20180916B has been detected at over 4 octaves in radio frequency, from 110 MHz (Pleunis et al. 2021a, hereafter P21) up to 5400 MHz (Chawla et al. 2020;Pilia et al. 2020;Marthi et al. 2020;Marcote et al. 2020;Nimmo et al. 2021;Pastor-Marazuela et al. 2021;Sand et al. 2022;Bethapudi et al. 2022).The Low-Frequency Array (LOFAR) High Band Antenna (HBA) detection of 18 bursts (P21) strongly constrains free-free absorption in the local environment of the source.These detections also demonstrated that FRB 20180916B's periodic activity is systematically delayed towards lower frequencies, in addition to detecting subtle, but measurable, RM variations of 2 − 4 rad m −2 (P21, 2 − 3% fractional variations).
The LOFAR telescope (van Haarlem et al. 2013) operates at the lowest radio frequencies () detectable from the Earth's surface (just above the ionospheric cut-off).At these low frequencies, propagation effects such as scattering (∝  −4 ), dispersive delay (∝  −2 ), and Faraday rotation (∝  −2 ) are much more pronounced.The temporal broadening due to increased scattering can lead to a reduction in the signal-to-noise ratio (S/N) of the bursts, making detection more challenging.Conversely, this also means that these effects can be measured and tracked for changes at very high precision with LOFAR -which offers ∼ 0.1 rad m −2 precision for tracking subtle, but still significant1 , RM variations of a few rad m −2 (P21).
The spectro-temporal properties of the bursts give insight into the local environment and emission mechanism.Like most repeaters, FRB 20180916B's bursts often show a ∼ 10 − 30% fractional bandwidth of emission and sub-burst drifting structure, argued to be caused by a radius-to-frequency mapping in the magnetosphere of a NS (Hessels et al. 2019;Wang et al. 2019;Lyutikov 2020).The bursts have been observed to be ∼ 100% linearly polarized at high frequencies (Nimmo et al. 2021), with no detectable circular polarization, and depolarizing only at the lowest frequencies (P21).This depolarization has been observed for other repeating FRBs as well (Feng et al. 2022), including the other possibly periodically repeating FRB 20121102A (Plavin et al. 2022).In addition, the flat polarization position angles (PPAs) that remain consistent even between bursts (Nimmo et al. 2021;Michilli et al. 2018b) disfavour rotational and precession models that expect PPA variations as a function of both the rotational and precession phases.
Recently, Mckinven et al. (2023a) observed a secular change in the RM of FRB 20180916B using CHIME/FRB, further demonstrating a progenitor residing in a dynamic magnetized environment.This large fractional RM change (of > 40% over a 9-month duration) came after measuring roughly stable RMs with comparatively small variations (< 10%) for nearly three years since its discovery in 2018 (CHIME/FRB Collaboration et al. 2019;P21).
In this paper, we present results from a monitoring campaign of FRB 20180916B with LOFAR at 110 − 188 MHz spanning 78 activity cycles spread across nearly three years.We examine whether the chromatic shift in activity window (P21; Pastor-Marazuela et al. 2021;Bethapudi et al. 2022) remains constant in time, look for changes in burst rate, and compare our bursts to those detected by CHIME/FRB at 400 − 800 MHz in the same time period.We also track scattering, Faraday rotation measure, and depolarization over multi-year timescales.Section 2 describes the observational setup and the burst search; Section 3 presents the analysis and results of our observations, whose implications are discussed in Section 4.

OBSERVATIONS AND BURST SEARCH
LOFAR has observed FRB 20180916B between 2019 June 6 and 2022 November 17, with observations of the source during 51 of the 78 16.33-day activity cycles in this time.Results of the observations between 2019 June 6 to 2020 August 6 (13 activity cycles), with 18 bursts detected, were reported in P21 (half of the same LOFAR bursts were reported in Pastor-Marazuela et al. 2021).
For the work presented here, we followed two observing strategies, totaling to 223 h of exposure on source, spread across 38 activity cycles: (1) Observations concentrated during the expected activity window as reported in P21, in order to maximize our chances of detecting high-S/N bursts to measure scattering timescales and RM variations.We followed this strategy between 2021 June to 2021 November, and 2022 June to 2022 November (89 h). ( 2) Observations distributed nearly uniformly throughout the 16.33-day activity cycle, in order to establish a well-bounded activity window and duty cycle.We followed this strategy between 2020 December to 2021 May, and 2021 December to 2022 May (134 h).
For each observation, signals from the HBAs of up to 24 LO-FAR Core stations were coherently added in the COBALT2.0 beamformer (Broekema et al. 2018) to form a single tied-array beam pointing towards the source.We used a pointing direction of  J2000 = 01 h 58 m 00.s 7502,  J2000 = +65 • 43 ′ 00.′′ 3152 (2.3-mas uncertainty), from Marcote et al. (2020), for all the observations.The ∼ 3 ′ Full-Width-Half-Maximum (FWHM) of the tied-array beam covered the uncertainty region of the FRB localization.Nyquist sampled complex voltage data for two linear polarizations were recorded for 400 subbands of 195.3125 kHz each, providing 78.125 MHz of bandwidth between 110 − 188 MHz, with typical observation durations of 1 h.
The complex voltage data were analysed using the pipeline described in more detail in P21; we briefly repeat the main aspects here.We channelized the complex voltage data to a frequency resolution of 3.05 kHz (64 channels per subband), after which the two orthogonal polarizations were squared and summed to form Stokes I total intensity data in filterbank format (Lorimer 2011).These data were averaged to a time resolution of 983.04 μs (after averaging by a factor of 3), and subsequently averaged in frequency by a factor of 16 using digifil (van Straten & Bailes 2011), and using incoherent dedispersion to a DM of 348.772 pc cm −3 , the best-fit DM from Nimmo et al. (2021).The temporal smearing due to incoherent dedispersion in a single channel at the given DM is between 1.3 ms (at 188 MHz) and 6.7 ms (at 110 MHz).Given the expected DM smearing, intrinsic burst widths (FWHM > 30 ms), and scattering timescales (∼ 50 ms) measured for this source in P21, we decide that a time resolution of 983.04 μs is sufficient to detect bursts from FRB 20180916B.This  approach to creating the dedispersed filterbank files is identical to that used in P21, such that the results can be compared.Note that the data were coherently dedispersed for subsequent measurement and analysis of burst properties, as presented in Section 3.3.
We used rfifind from the PRESTO software suite (Ransom 2001) to identify radio frequency interference (RFI); we replaced the values of these frequency channels and time integrations with random noise with the appropriate local statistics.From these RFI cleaned files, dedispersed time-series filterbanks were generated using the graphics processing unit (GPU)-accelerated Dedisp library (Barsdell et al. 2012) for DMs ranging from −20 to +20 pc cm −3 around the nominal DM of FRB 20180916B, in steps of 0.1 pc cm −3 .FRB 20180916B and other repeating FRBs often show bursts with fractional bandwidths of ∼ 10 − 30% (Gourdji et al. 2019;Pleunis et al. 2021b;Kumar et al. 2021).To increase the sensitivity to such bursts, dedispersed time-series were also created for three overlapping halves of the band (110 − 149, 130 − 169, and 149 − 188 MHz), and seven overlapping quarters of the band (110 − 129, 120 − 139 MHz, etc.).
Candidate bursts were generated by cross-correlating the dedispersed time-series with top-hat functions with widths up to 150 ms, using a GPU-accelerated version of PRESTO's single_pulse_search.py3 .The diagnostic plots from single_pulse_search.py showing the S/N of bursts versus DM and time were inspected by eye for bursts, as a fail-safe.The events identified by single_pulse_search.py at separate DMs were grouped using SinglePulseSearcher4 (Michilli et al. 2018a) to eliminate redundant burst candidates obtained from dedispersing the same burst at slightly different DMs.This was done by choosing the highest-S/N burst in each group containing more than 4 candidates with S/N > 5, and rejecting the rest.The remaining candidate bursts were assessed using the FETCH5 (Agarwal et al. 2020) deep-learning-based classifier.FETCH is trained to classify candidates based on their frequency-time image (i.e., dynamic spectrum) and their DM-time image (i.e., bow-tie plot).Since the S/N drops steeply with changing DM at the low radio frequencies of our observations (Cordes & McLaughlin 2003), the characteristic DM-time bow-tie variation would not be visible at large DM deviations.Hence, we modify the DM range for the DM-time image to go from [0.9DM, 1.1DM] instead of the default [0, 2DM] used in FETCH.All FETCH candidate bursts with a single_pulse_search.py reported S/N > 6 -and with a mean probability of being astrophysical > 50% from the 11 available FETCH deep-learning models -were visually inspected to remove possible contamination by residual RFI.
In our observations, we detect 11 new bursts, bringing the total number of bursts detected from FRB 20180916B by LOFAR to 29.The bursts have been labeled B1 through B29 (of which B1−B18 were reported in P21), in chronological order of barycentric arrival time.The bandpass-corrected dynamic spectra and time series profiles of bursts B19−B29 are plotted in Figure 1 and the corresponding burst parameters can be found in (Table 1).

Bursting activity
Figure 2 shows the barycentric times of arrival, total bursts per 16.33day activity window, and burst rates of all 29 LOFAR bursts along with the CHIME/FRB bursts detected since the first LOFAR observations of FRB 20180916B began (2019 June 6), plotted against the activity phase.Activity phases are determined by setting MJD 58369.40 at phase  0 (i.e,  = 0.0), such that  = 0.5 is the mean of the folded phases of the CHIME/FRB bursts in P21.We find that the 11 new bursts fall within the previously observed 4.1-day (25% of the 16.33day cycle) LOFAR activity window (P21) between activity phases 0.53 <  < 0.79.
In order to examine if the properties of the activity window remain stable, accounting for the varying exposure of the LOFAR observations, we repeat the Markov-Chain Monte Carlo (MCMC) analysis from P21 (their Section 3.5, Figure 10) for bursts from different time ranges.We briefly repeat the description of the analysis here, which uses emcee by Foreman-Mackey et al. (2013) for the MCMC analysis.For each pair of the central window phase  LOFAR c and activity window width ,  burst arrival times are randomly drawn in the windows of allowed activity in the specified epoch, assuming a uniform burst rate across activity cycles and across the allowed activity phases.The number of bursts  for a given epoch is based on the known average burst rate in the window of observed activity.This analysis assumes that there is no significant derivative of the activity period, and random values for  LOFAR c and  are drawn from uniform distributions with 0.0 ≤  LOFAR c ≤ 1.0 and 4 ≤  ≤ 10 day.We thus arrive at a probability distribution function of obtaining bursts that fall within the observed phases, given the number of realizations out of 10,000 draws in which all drawn bursts fall within the observed activity phases.The posterior distributions were determined by using 20,000 steps of 32 walkers until the chains converged.The best fit parameters were obtained from the posterior distributions after discarding a burn-in phase of 1,000 steps and thinning the chains by 5 steps.
These simulations to fit for  LOFAR c and  in the 110 − 188 MHz band are done for the following four time ranges and the distributions with fits shown in Figure 3: (a) From MJD 58640 to MJD 59087, with 18 detected bursts.The best fit parameters are  LOFAR c = 0.72 +0.07 −0.04 and  = 5.02 +2.34 −0.78 (restating the window properties from P21); (b) From MJD 58640 to MJD 59854, with 29 detected bursts (this includes all LOFAR observations of and bursts from FRB 20180916B to date).We find that the most probable activity window width is  = 4.32 +0.67  −0.24 day, which is more tightly constrained than previously found in P21.The window is centered at a phase of  LOFAR c = 0.67 +0.03 0.02 , which is earlier than the previously obtained value by nearly a day (although still consistent within errors).This is still significantly offset from the CHIME/FRB window centered at a phase 0. and  = 6.66 +2.19  −1.90 .The implications of the activity window properties remaining consistent before and after the start of the secular RM trend, i.e., cases (c) and (d), are discussed later in Section 4.
Burst rates are calculated across time -for each 16.33-day activity period (only considering exposures within the observed LOFAR activity window), and across activity phase -including all cycles during which we had observations.Assuming that the bursts detected within the activity window follow a Poisson process, we calculate the 1 errors following Gehrels (1986).The obtained rates are plotted in panels 1c and 2c of Figure 2, with upper limits for epochs with no detected bursts.We find that the burst rates are largely consistent between activity cycles within 1 (panel 1c, Figure 2).By distributing our observations across all activity phases (as described in Section 2), we are able to place upper limits on the LOFAR burst rate throughout the activity cycle, outside of the expected LOFAR activity window.The average burst rate and the 3- Poissonian errors during the activity window we obtain in this work, i.e., case (b), is 0.17 +0.28  −0.09 h −1 , which differs from the upper limit on burst rate outside the expected window of 0.006 h −1 , by over 3 (see also panel 2c of Figure 2).

Scattering timescale and sub-burst drift rate variations
Scattering timescale fits can be covariant with other frequency-time properties of the bursts such as DM, sub-burst drift rate (Hessels et al. 2019), and the intrinsic width of the underlying burst (as it was See Section 3.2 for fitting method used.b  low and  high are determined by fitting a Gaussian profile to the spectrum of the burst and taking the ±3 limits of the fitted profile.If the  low or  high thus obtained are beyond the edges of the observing band, the edge of the observed band occupied by the burst is quoted.c Profile integrated S/N of the time series profile, after summing over the frequency range occupied by the burst.
Note that this value is different from the peak S/N and our detection metric single_pulse_search.py-based S/N.d Fluence determined using only the frequency channels occupied by the burst.
The procedure followed for fluence calculation and fluence error estimation is the same as described in P21.
± 27 a Scattering timescales are only quoted for single-component bursts with S/N > 10.
The errors quoted here are only statistical from the MCMC fit.We incorporate deviations in scattering from MCMC fits to simulated bursts to account for covariance with drift rates as an additional source of error in Figure 4. See Section 3.2 for a full description.b Drift rates using the ACF method are only quoted for bursts with S/N > 10. c Drift rates using the time of arrival method are only quoted for bursts with S/N > 20.
The S/N limit here is higher than for the ACF method since we are more limited by S/N by dividing individual bursts into subbands to calculate the times of arrival at each subband.d Observed RM, uncorrected for ionospheric contribution.Error in RM calculated as We assume here that bursts from this source are intrinsically made up of multiple sub-bursts that drift towards lower frequencies as time progresses.In the limit of low S/N, the individual sub-bursts blend together to produce one 'blob' of emission (see, e.g., Figure 7 in Gourdji et al. 2019).When measuring DMs in this low S/N limit, one is biased high due to any dedispersion algorithm (optimizing for burst S/N or burst structure) preferring to superimpose the sub-bursts that are drifting downwards in frequency.This is seen in the middle panel of Figure 3  Given the absence of evidence for significant DM variations, we choose to dedisperse all bursts to the same structure-optimized DM = 348.772pc cm −3 from microsecond burst structure measured by Nimmo et al. (2021).Using a DM value optimized for microsecond structure circumvents the issue of covariance of the DM with spectrotemporal changes of the burst profile, particularly exacerbated by the obscuring of structure by scattering at lower frequencies.Upon dedispersing the dynamic spectra, we can measure the times-ofarrivals of the burst at different parts of the burst bandwidth, to check for residual quadratic delays from dispersion. Figure A1 shows these delays versus frequency for the brightest bursts (S/N > 20).Since the times-of-arrivals progressively move to later times at lower frequencies for all the bursts, any variations in DM must all be in the same positive direction compared to DM = 348.772pc cm −3 , which we find to be unlikely.In the case of stochastic DM variations, we would expect that the DM randomly fluctuates around the used DM, rather than in one direction preferentially.We see that deviations in DM can be no greater than +0.25 pc cm −3 , although we note that most of the delays in individual bursts appear non-quadratic and thus interpret it as sub-burst drifting instead.The delays being larger for bursts at the bottom of the band can be explained by the expected linear scaling of drift rate with frequency, as we discuss in more detail below.
We are unable to perform a joint-fit for the scattering timescales and drift rates due to insufficient S/N of the bursts in the frequencytime dynamic spectra.Hence, we fit for these two quantities separately and calculate the uncertainties due to drift-scattering covariance using simulated bursts.Downward drifting sub-bursts towards later times and lower frequencies have been seen to be characteristic of repeating FRBs (Hessels et al. 2019; Fonseca et al. 2020).The two distinct components visible in the dynamic spectrum of B25 show that this 'sad trombone effect' is manifest in bursts from FRB 20180916B even at LOFAR frequencies.In the case of bursts that do not show distinct sub-bursts, yet show delays towards lower frequencies in the dynamic spectra post dedispersion, we assume that the downward drifting sub-structures have been obscured by low S/N, or shorter waiting times, or large scattering timescale, or any combination of the above.To measure the linear drift rate of the bursts (   or ), we fit a two-dimensional Gaussian ellipsoid to the two-dimensional auto-correlation function (ACF) of the dedispersed dynamic spectra of the bursts (see Figure A2 for the 2D ACF fits obtained).The 2D-Gaussian ellipsoid is parameterized as follows, after applying a counterclockwise rotational transformation by an angle : Here,   and   are the standard deviations of the Gaussian ellipsoid along orthogonal axes after the rotational transform is performed.The inclination of this ellipsoid cot(), gives us the drift rate in MHz ms −1 .We use the lmfit Python package to perform a non-linear least-squares fit for the free parameters A,   ,   and .The uncertainty on  is calculated by taking the square-root of the variance from the fit, summed with the square-root of the reduced- 2 of the residuals from the fit as an additional source of error.
Previously, the magnitude of the drift rate has been seen to roughly linearly decrease with decreasing frequency for FRB 20180916B as well as FRB 20121102A (Chamma et al. 2021;Hessels et al. 2019).For FRB 20180916B, we are able to recover this relation of decreasing drift rate magnitudes with frequency, even within the 78-MHz LOFAR bandwidth (Figure 5).We obtain the relation  = (−0.02± 0.01)  from a least-squares fit of a line to the mean value of the drift rates within the LOFAR band (without weighting the drift rates by their errors).Considering measurements from higher frequencies (Chawla et al. 2020;Chamma et al. 2021;Pastor-Marazuela et al. 2021;Sand et al. 2022), the drift rate scales as  = (−0.027± 0.005)  for this source between 110 MHz and 1520 MHz, which agrees with the  = −0.02 relation obtained by Sand et al. (2022) for drift rates measured between 400 − 1500 MHz.
As a second method to measure the drift rate, we use the slope of a linear fit to the times-of-arrival at different frequencies (see Figure A1).The drift rates measured this way are also plotted in Figure 5, and quoted in Tables 2 and B1.They are found to be consistent within errors for most bursts where the drift rates were measured using the two-dimensional ACFs, and they follow an inverse linear trend with frequency with a slope of −0.04 ± 0.03.
Inhomogeneities in the medium through which the bursts propagate results in a geometric time delay between rays from the burst received at the telescope, manifesting as a temporal broadening of the pulse.This pulse broadening due to multi-path scattering can be described by a one-sided exponential function in the case of an infinitely extended, thin screen.The resulting time () variant burst profile I() is a convolution of the pulse broadening function with the underlying emitted burst profile: Here,  s is the scattering timescale and Γ(,   ) is the underlying burst profile before scattering, which we assume to be a Gaussian with standard deviation   .Scattering timescale varies with frequency as  s ∝   , where generally scattering index  is between −4 (Williamson 1974;Lang 1971, for a thin-screen scattering model with Gaussian inhomogeneities) and −4.4 (Lee & Jokipii 1975, Kolmogorov spectrum).However, it is possible that scattering scales with frequency with a different .To test this we measure scattering timescales at different subbands of individual bursts by performing a scipy least-squares fit of a model of the burst as given in Equation 2to the different subbands (Figure A4), and measure varying indices between −4.5 ≤  ≤ −1.7.Given the large errors on the indices (see Figure A4), we choose to scale the scattering timescales measured below to an  = −4.0 to reference them to 150 MHz, and include the deviation in scattering timescale if scaled using  ∈ [−4.4,−3.6] as a part of the error.The  fit for burst B6 has the smallest error associated and agrees with  = −4.0.
Scattering timescales were measured by fitting the total intensity profiles of each burst with a symmetric Gaussian envelope ( Γ(,   )) convolved with an exponential tail first by a least-squares fit using the scipy Python package.The obtained values were used as the initial values for a least-squares fit using the MCMC method (with the emcee Python package; Foreman-Mackey et al. 2013).The results of the fitting algorithm can be compared with the burst profile in Figure A5, which also shows the residuals after subtracting the model fit from the data.One Gaussian component was assumed for all bursts except B25 where the downward-drifting substructures are distinguishable in the dynamic spectrum (Figure 1).We assume that scattering time scales with frequency as  −4 , and only measure it for bursts with S/N > 10.Since our scattering fits to bursts B7 and B14 (see Figure A5) are affected by the presence of possible secondary peaks in their frequency-integrated time profiles, we consider the measured value for these bursts to only be an upper limit.The measured scattering timescales (Tables 2 and B1), referenced to 150 MHz, span a range between ∼ 30 ms and ∼ 70 ms (excluding B7 and B14), as seen in Figure 4.
To examine the effect of drifting on the scattering timescales obtained through such a fitting process, we fit scattering timescales to the time profiles of simulated bursts with non-zero drift rates in their dynamic spectrum, using the same fitting process.The dynamic spectrum of each of the simulated bursts was modeled as a single tilted 2D Gaussian, as described in Equation 1 (since most of the detected bursts do not show distinct components), after which the channels containing the burst had their Gaussian profiles convolved with an exponential scattering tail as described in Equation 2. The simulated bursts had the following input properties: • Standard deviation of the Gaussian intrinsic burst width in time   = 16 ms, the median value from the instrinsic widths of B1−B29.
• Standard deviation of the Gaussian profile extent of the burst in frequency   = 6 MHz, the median value from the frequency extents of B1−B29.
• Scattering timescales   referenced to 150 MHz, sampled uniformly in the range [0 ms, 150 ms], at the central frequency of the burst.
• Time series S/N (signal summed over the extent of the burst and divided by the square root of the width of the burst in time) of 20.
• Equal fluences, by normalizing the time series by the area under the burst.
Figure A3 shows the fractional difference between the fitted scattering timescale (  (fit) ) and simulated scattering timescale (  (input) ), for varying input drift rates.We treat this fractional difference as the covariance of the fitted scattering timescales with drift rate.The measured drift rates for bursts B1−B29 are then matched with the input drift rates of the simulated bursts, and the measured scattering timescales are matched with the fitted scattering timescales from the simulation.An additional source of one-sided-error (based on whether the fit is over/under-estimated in relation to the input scattering in the simulated bursts) on the measured scattering timescales (  (meas) ) of bursts B1−B29, given the measured drift rate is derived as : •   (meas)

Fractional polarization and RM variations
LOFAR has dual-linear polarized feeds that are stationary on the ground and aligned along the South-East to North-West and South-West to North-East directions.Depending on the elevation of the source at the time of observation, the electromagnetic waves from the burst will be projected on to the stationary telescope feeds at an angle, degradation of signal intensity and distortion of the polarization signal (Noutsos et al. 2015).We apply a Jones calibration matrix calculated from a beam-model of LOFAR using the dreamBeam6 package to mitigate this effect.While most bursts remained mostly unaffected by this correction (owing to high elevation of the source during detection by the telescope), it reduced the fraction of linear polarization detected for a few of the bursts by up to 12%.Owing to the millisecond-duration of our bursts, effects of changing projection effects for the dipole length as the source is tracked in the sky will be negligible.LOFAR real-time calibration corrects for cable and geometric delays during beam-forming.We use ionFR 7 (Sotomayor-Beltran et al. 2013) to calculate the ionospheric contribution to the measured RM (RM iono in Table 2).For polarization analysis, the complex voltage data was coherently dedispersed to a DM= 348.772 pc cm −3 using dspsr (van Straten & Bailes 2011) at a time resolution of 983.04 μs and frequency resolution of 12.2 kHz.At this frequency resolution, depolarization due to uncorrected intra-channel Faraday rotation at the lowest-observed frequency in the band (109.88MHz) is less than 1% at the maximum absolute RM value reported for FRB 20180916B by P21 (and < 12% up to an RM= 1000 rad m −2 ).The linear hands give us the correlated polarization products: XX, XY, YX and YY that are used to calculate the 4 Stokes parameters , ,  and , with  2 =  2 +  2 +  2 .Linear polarization is computed as  meas = √︁  2 +  2 .The true linear polarization is obtained after de-biasing the linear polarization  meas based on the prescription in Everett & Weisberg (2001), where: Faraday rotation occurs as the electromagnetic wave passes 7 https://github.com/csobey/ionFRthrough a magnetized plasma.It has the effect of rotating the plane of polarization of a linearly polarized wave, with the degree of rotation scaling as the inverse square of the radio frequency.For a cold, magnetized plasma, the change in PPA due to Faraday rotation (Ψ) can be expressed as a function of the wavelength () as follows: where PPA Ψ = 1 2 tan −1   (periodic as  in phase); and RM quantifies the slope of the PPA versus  2 relation.
We measure RM by implementing an RM-synthesis (Burn 1966;Brentjens & de Bruyn 2005) algorithm.The values obtained using our code were verified using the RM-tools software8 (Purcell et al. 2020).Burn (1966) define the Faraday Dispersion Function (FDF) as the Fourier transform of the complex linear polarization vector per unit RM.The measured linear polarization vector is modified/multiplied by the 'RM Spread Function' (RMSF) as a result of the finite window of bandwidth of the telescope.This window determines the resolution obtainable in RM-space, and is given by the FWHM of the RMSF.For the LOFAR setup used in this work, the RM resolution is ≈ 3.8  2 max − 2 min = 0.77 rad m −2 .In the case of FRBs, a Faraday-thin FDF is expected to first order, since we expect a near-perfect point source behind a Faraday screen (Michilli et al. 2018b).In this scenario, all the polarized flux is Faraday rotated by the same value -i.e., the FDF can be approximated to a delta function at the value of a single RM, convolved with the RMSF.We obtain the RM value from the single peak of the Faraday-thin FDF.By definition, Faraday rotation affects Stokes  and , being seen as = (-0.027±0.005) c LOFAR (ACF method, this work) LOFAR (subband least-squares method, this work) GBT, CHIME/FRB and GMRT (Chawla+2020 and Sand+2022) CHIME/FRB (Chamma+2021) GBT ( Sand+2022) Apertif (Pastor-Marazuela+2021) Figure 5. Drift rates versus the centre frequency of the bursts in the LOFAR band (inset: orange circle markers (ACF method), black circle markers (time of arrival method)), compared to the measured rates at higher frequencies -GBT, CHIME/FRB, and GMRT 300 − 500 MHz (Chawla et al. 2020;Sand et al. 2022): dark blue triangle markers; CHIME/FRB 400 − 800 MHz (Chamma et al. 2021): blue square markers; GBT 600 − 1000 MHz (Sand et al. 2022): yellow diamonds; Apertif 1220 − 1520 MHz (Pastor-Marazuela et al. 2021): green pentagonal markers.The 1- errors of the drift rates are included.We measure drift rates for all detected FRB 20180916B LOFAR bursts B1−B29, i.e., including bursts reported in P21.We obtain the relation  = −0.02± 0.01  from a least-squares fit of a line to the drift rates in the LOFAR band.The linear decrease in drift rate (using the ACF method) with the observed frequency measured within the LOFAR band is consistent within errors with the best-fit linear relation of  = (−0.027± 0.005)   for this source between 110 MHz and 1520 MHz.Neither of these fits weight the drift rates by their individual uncertainties.oscillations in the  and  spectra due to the rotation in the complex linear polarization vector ( + ) as a function of  2 .We obtain the total linear polarization and the PPA of the pre-Faraday-rotated burst by 'derotating' the  and  spectra (integrated over the burst duration) by the determined RM as shown below: where Ψ 0 is the PPA at infinite frequency.RM measurement was performed using only the corresponding frequency channels of the  and  spectra occupied by the Stokes  spectrum of the burst in order to increase S/N.The lack of polarization in B26 and B29 could in part be due to their very low S/N (< 8).For B24 and B25 (note that B25 is at the top of the LOFAR band), we obtain no RM measurement despite higher S/N and lower scattering.The two visible components in B25's dynamic spectrum were also searched separately for a peak in the FDF.For these bursts that appear completely depolarized in the LOFAR band, we place conservative upper limits on / after derotating the  and  spectra by an 'expected' RM value using the best-fit slope of 0.197 rad m −2 day −1 from Mckinven et al. (2023a), during their observed secular RM trend (since the true RM, assuming the intrinsic emission was polarized, is not known in these cases).
P21 reported on the polarimetric properties of 3 out of the 18 bursts they presented, since 14 of the bursts (B5 − B18) had only totalintensity data available.In this work, we detect linear polarization from 4 out of 11 bursts and measure their RMs (bursts B19, B20, B23, and B27).We detect no significant circular polarization in any of these bursts.The remaining 7 bursts appear unpolarized.Polarization profiles of linear polarization, , and circular polarization, Stokes , are plotted in the mid-panel of each sub-figure in Figure 1 along with their total-intensity (Stokes ) profile.The polarization profiles shown in Figure 1, like the total-intensity profiles, were only summed over the frequencies where the burst is present in the band.In the top panels for each burst are the corresponding PPAs.
The measured RMs, assuming a  2 scaling of the PPAs (and the corresponding ionospheric contribution to the RM, RM iono ), along with the linear polarization fraction measured at this RM can be found in Table 2.We plot the measured RMs as a function of time and activity phase in Figure 6.Until B22 on 2021 March 07, bursts from FRB 20180916B showed subtle but detectable RM variations in the LOFAR band of ∼ 2 − 3 rad m −2 .For burst B23, we measure an RM = −105.6rad m −2 , i.e., an absolute change of 10.9 rad m −2 and a ∼ 8% fractional change in the RM compared to the previous LOFAR measurement from B20.When over-plotted with RM measurements at 400 − 800 MHz from Mckinven et al. (2023a), this apparently abrupt change in RM agrees with the secular trend in RMs observed in bursts detected with CHIME/FRB.The most recent RM measurement with LOFAR for burst B27, with RM= −54.4 rad m −2 , points towards a decrease in the RM gradient with time.
For LOFAR bursts from FRB 20180916B that had sufficient S/N (including B1−B3, published in P21), we were able to divide the bursts into halves/quarters in the frequency band they occupied, and measure the linear polarization fraction in each of these frequency chunks.An example of the Stokes , , and Stokes  profiles for burst B1 in quarters of its total occupied bandwidth is shown in Figure 7 (top).This allows for quantifying the variation of polarization fraction with frequency within the same burst.Figure 7 (bottom) has the split-bandwidth linear polarization fraction measurements, as well as the full-bandwidth values plotted against central frequency.It shows that FRB 20180916B bursts exhibit depolarization towards lower frequencies not only between bursts as reported in P21, but also within the individual bursts.It is also clear from Figure 7 (bottom) that the frequency of depolarization is not constant between bursts in the LOFAR band.
We further check for RM jumps in time across the pulse profile (as seen in some pulsar profiles in Noutsos et al. 2009) but do not find any evidence for this.However, we do find small variations in RM for the same burst at different parts of its frequency band.Any additional dependence of RM on wavelength/frequency, when fit for RM ∝  2 , can be seen as a deviation from this canonical  2 scaling.We find hints of deviation from the expected  2 scaling for one of the bursts; further analysis of this will be explored in a future paper.

DISCUSSION
We see that the frequency-dependent activity window of FRB 20180916B at LOFAR remains stable in time.With our improved constraints on the activity window at 110−188 MHz, we find that the width of this window (4.3 +0.7 −0.2 day) is likely smaller than the 5.2-day window in the 400 − 800 MHz range of CHIME/FRB.This is unlike what Pastor-Marazuela et al. (2021) find from observations at higher frequencies, where the activity window at the lower frequencies of the 600 MHz CHIME/FRB band is wider than that at the higher frequencies of 1370 MHz with Apertif.Bethapudi et al. (2022) detected the highest-frequency bursts from FRB 20180916B to date, by extrapolating the chromatic activity to frequencies of 4 − 6 GHz to pick observing windows; however, more observations distributed across all activity phases are required to constrain the width of the window at these frequencies.Along with the new LOFAR window constraints we present here, this suggests that while the systematic delay in activity spans nearly 6 octaves in frequency, a systematic increase in activity width towards lower frequencies does not hold true for the entire frequency range of detections to date.

Implications of the observed propagation effects
We find that the drifting rate of sub-bursts (the 'sad trombone' effect) does not depend on the activity phase, which means its cause should be decoupled from the mechanism causing the 16.33-day periodicity.The drifting could be caused by intrinsic emission processes, like radius-to-frequency mapping of the emission region (Hessels et al. 2019).Marcote et al. (2020) find a scintillation bandwidth of 59 ± 13 kHz at 1.7 GHz, corresponding to a scattering timescale of 0.015 ms at 1 GHz, after correcting for plane-wave scattering of an extragalactic burst (Ocker et al. 2021).This value is found to be consistent with the expected line-of-sight scattering from the interstellar medium (ISM) in the Milky Way disk based on the NE2001 model (Cordes & Lazio 2002).At 150 MHz this spans a range of scattering timescales between 14 ms (scattering index  = −3.6)and 63 ms (scattering index  = −4.4).The total scattering observed will be a linear sum of the line-of-sight contributions: , where  s,mw ,  ,s,halo ,  ,host ,  ,frb are scattering contributions from the ISM of the Milky Way, the Milky Way halo, the ISM of the host, and from the circum-source medium of the FRB, respectively. ,s,halo is < 4 ms at 150 MHz based on the estimate from Ocker et al. (2021).We find that  s,mw and  ,s,halo can account for nearly all of the scattering we measure in the bursts at 150 MHz, with minimal contribution from  ,host .We find hints of a few ms-timescale variations on weeks to months timescales that could arise from the local environment of the FRB.We do not find strong evidence for longer timescale variations on years timescales, that could possibly arise from the relative motion of the observer and the Milky Way and host galaxy ISM.
Up to factor-two variations in scattering time over a few minutes to days timescales have been reported for another repeating FRB, FRB 20190520B (Ocker et al. 2023).Such variations are also observed in the Crab pulsar, accompanied by DM variations, and are attributed to arise from the discrete structures ('filaments') of the Crab nebula (McKee et al. 2018;Driessen et al. 2019).
Rotation measure variations (> 5% fractional), both drastic and moderate, have been previously seen in various other repeating FRBs: FRB 20121102A (Michilli et al. 2018b;Hilmarsson et al. 2021b;Plavin et al. 2022), FRB 20180301 (Luo et al. 2020), FRB 20201124A (Hilmarsson et al. 2021a;Xu et al. 2022;Jiang et al. 2022), FRB 20190520B (Anna-Thomas et al. 2022), and 12 additional repeating FRBs as shown in Mckinven et al. (2023b).FRB 20121102A and FRB 20190520B exhibit extreme RM variations > 10 4 rad m −2 , and both have a compact persistent radio source associated with them (Chatterjee et al. 2017;Marcote et al. 2017;Niu et al. 2022).Out of these, the two periodically active repeaters FRB 20121102A and FRB 20180916B reside in vastly different types of galaxies, while both being close to star-forming regions.Until MJD 59355, FRB 20180916B was in fact one of the only repeaters to exhibit relatively minor RM variations (2 − 3 rad m −2 ; P21) in comparison, despite these variations being higher than typically observed towards Galactic pulsars.While it is not possible to comment on the RM evolution of non-repeating sources, repeating FRBs seem to favour residing in dynamic magneto-ionic environments.
The phase centre and width of the activity window remain consistent between bursting activity observed before and after MJD 59355 (although more detections are required to obtain better constraints after MJD 59355), when the first burst with an RM change > 10 rad m −2 was detected by Mckinven et al. (2023a).We do not detect any obvious change in the chromatic activity window at LO-FAR accompanying the large RM variations following relatively stable RM measurements over ∼2.4 years.Consistent values between CHIME/FRB and LOFAR RM measurements (including the latest RM measurement reported in this work, as conveyed to us by private correspondence with Ryan Mckinven) show that the RM variations are unlikely to be dependent on the frequency band in which  the bursts are detected.This implies that the emission at different frequencies comes from the same Faraday depth, thus placing the Faraday screen(s) where the RM variations originate, in front of all the emission region(s) at different frequencies.
In the absence of any significant correlation between the small (< 0.8 pc cm −3 ) DM variations with the RM variations, Mckinven et al. (2023a) attribute the observed RM evolution to changes in the strength/direction of the magnetic field in the local medium producing Faraday rotation (and not changes in the electron density of the Faraday screen).They calculate limits on the rate of change in the line-of-sight component of the magnetic field between a fraction of a G/day and ∼ 10 5 G/day, depending on the assumed DM contribution from the Faraday-active medium (the greater the assumed DM contribution, the lower the required change in the magnetic field to account for a given change in RM).
Bursts from FRB 20180916B typically have lower fractional linear polarization towards lower frequencies, even within individual bursts (Figure 7).However, there is a large scatter in the fractional polarization between different bursts, even at the same radio frequency (Figure 7).For instance, consider bursts B1, B3, B19, B23, B27 which are centered within a 10-MHz range between 157−167 MHz.Despite being close in emission frequency, these bursts show variations in linear polarization fractions in the range ∼ 34% to ∼ 60%.B25, B28 and B29, also centered at the top of the band show no detectable polarization at all.This is comparable to the range of polarization fractions observed across the entire 80-MHz LOFAR band.Beniamini et al. (2022) discuss depolarization effects on FRBs caused by multi-path propagation through a magnetized scattering screen.They expect a a power-law dependence of the polarization fraction, decreasing with frequency.Feng et al. (2022) found that a model of exponential decrease in polarization with frequency due to a scattering screen can explain the decrease in / at lower frequencies for five repeating FRB sources.The analysis used by Feng et al. (2022) applies a model originally developed by Burn (1966) and extended in a series of works since then (e.g., Sokoloff et al. 1998;Brentjens & de Bruyn 2005).These works focus on depolarization of incoherent synchrotron radiation.Their results, and in particular the exponential cutoff of the polarization below a certain frequency threshold differ qualitatively and quantitatively from the depolarization that will be suffered by coherent radiation passing through a scattering screen, as explored in Beniamini et al. (2022).The latter is the appropriate case for FRBs, whose enormous brightness temperature implies that the radiation must be coherent.
It is possible that the large scatter in the polarization fractions of bursts centered around similar frequencies, when bursts from multiple activity cycles are considered, arise due to different lines-of-sight through a Faraday active scattering screen in the local environment of the source.There might still be more constrained trends within each activity cycle that get washed out when bursts from multiple activity cycles are considered.On the contrary, we do not find a significant correlation between the scattering timescales of the bursts (at the central frequency of the band occupied by the burst) and linear polarization fractions of the bursts (see Figure A6).
In a scattering model of depolarization, the depolarization is caused by an external Faraday screen separate from the emitting region.Depolarization can also occur in a Faraday-active medium co-spatial with the emitting region.The plane of polarization undergoes differential Faraday rotation in the emission region itself, causing depolarization when summed over the extended region.The fractional depolarization (f int depol ) due to this 'internal depolarization' model (Gardner & Whiteoak 1966;Burn 1966;Sokoloff et al. 1998)  .Top: Total intensity (black) and polarization profiles (linear: red; circular: blue) for a single burst (B1), at different radio-frequency ranges of the total spectral range occupied by the burst.The band occupied by the burst was divided into four quarters labelled Q1, Q2, Q3, Q4 from the lowest to highest frequency, respectively.The profiles for each quarter have been normalized by the peak of the total intensity, for visual purposes.The respective linear polarization fraction measured in these quarters can be found in the figure below.Bottom: Linear polarization fractions versus the central frequencies of the bursts.Additionally, B1 was split in bandwidth into quarters, and B3, B4, B19 into halves; and the linear polarization fractions of these parts plotted separately.We see that the linear polarization decreases with decreasing frequency even within the spectral extent of a single burst.Overplotted is the internal depolarization model (see Section 4.1 and Equation 7).Additionally, we see that bursts B25, B28, and B29 are completely depolarized despite occupying the top-half of the frequency band.Upper limits are calculated for the bursts with no detected polarization (B21, B22, B24, B25, B28, and B29) at an interpolated RM value based on the best-fit slope for the RM secular trend from Mckinven et al. (2023a).can be described as: for the simple case of a uniform slab with a regular magnetic field and Faraday depth R emit , where the differential Faraday rotation Δ due to the Faraday depth of the slab is R emit  2 .By eye, we find that R emit of the emitting region would occupy a range between ∼0.03 and ∼0.45 rad m −2 (Figure 7).The best-fit value of 0.18 rad m −2 at the 600-MHz central frequency of the CHIME/FRB band (Mckinven et al. 2023a) also lies within this range.Note that, in this case, no correlation between the temporal scattering (caused by external foreground screens) and depolarization fraction is expected, which agrees with our findings (Figure A6).Plavin et al. (2022) found that internal depolarization by an emitting region of Faraday depth 150 rad m −2 could explain the depolarization seen for bursts from the possibly periodically active, repeating source FRB 20121102A from 1.7−5 GHz.Given that the emission at 150 MHz and 600 MHz are likely from different emission regions (since they are band-limited and happen at different activity phases), the RM evolution is likely still caused by a different Faraday screen along our line-of-sight to the different emission regions.

A self-consistent model for the observed propagation effects
FRB 20180916B is periodically active for a total 8.8-day window every ∼16.3 days, first at higher frequencies and a few days later at lower frequencies (from 6 GHz down to 110 MHz).Explanations for the periodic activity of FRB 20180916B have so far come in three classes: a binary orbit, precession, or the rotation of a compact object.We will briefly review these three hypotheses and discuss them in context with observational constraints on the source.Binary orbit: Tendulkar et al. ( 2021) favour a high-mass X-ray binary model (with a late O-type or B-type companion) given FRB 20180916B's 16-day periodic activity and inferred age of O (10 6 yr) based on the observed offset from a star-forming knot in its host galaxy.They suggest that the modulation of activity window with frequency in such a scenario would come from free-free absorption in the swept back wind of the companion (Ioka & Zhang 2020;Wada et al. 2021).Other proposed binary models include a Be/X-ray binary system (Li et al. 2021b) with a NS and Be-star, and an eccentric binary system consisting of a magnetar and early-type star (Barkov & Popov 2022).In Li et al. (2021b), interactions between the circumstellar disk of a companion Be-star and the NS produce starquakes on the NS resulting in FRBs, while free-free absorption in the disk causes a chromatic activity window.If the RM evolution is the result of the periastron passage of the NS through such a disk, then the orbital period must be greater than the 3-year timescale of RM evolution observed for the source thus far.Such a system would need an alternate explanation for the observed 16.3-day activity period.A well-known example of a similar system is the Galactic binary containing PSR B1259−63 and a Be star.In the 50 days around the periastron in its 1237-day orbit, changes in DM (∼ 6−8 pc cm −3 ), RM (∼6000 rad m −2 ), and an increase in scattering (such that the pulsar gets eclipsed by the companion disk for 35 days) were observed (Johnston et al. 1992(Johnston et al. , 1996(Johnston et al. , 2005)).These drastic changes in the propagation effects undergone by the pulses are thought to arise from the stellar disk of the Be-star companion.Alternatively, an eccentric binary system where the emission frequency depends on the orbital separation, along with viewing angle limitations from beaming geometry, could explain the asymmetric chromatic activity windows seen for FRB 20180916B.
Slow rotation: Beniamini et al. (2020) proposed that FRB 20180916B is an ultra-long-period magnetar that rotates every ∼16.3 days.In this scenario, the emission is associated with a polar cap beam with a high-to-low-frequency chromaticity from radius-tofrequency mapping (see also the simulations in Li et al. 2021b).Any DM variations should not vary periodically with the rotational phase.The recent discoveries of radio emission from sources with periodicities of 76 s (Caleb et al. 2022) and 1091 s (Hurley-Walker et al. 2022) and a 6.7-h magnetar candidate (De Luca et al. 2006) add credence to the existence of a population of ultra-long-period magnetars, that most surveys are biased against detecting (Beniamini et al. 2023).However, a 2 week period is still a stretch from the minutes-hours rotational periods that have so far been observed to exist.
Precession: Precession of the FRB-emitting NS/magnetar may arise either from a freely precessing, non-spherical star deformed by its internal magnetic field (Levin et al. 2020;Zanazzi & Lai 2020); orbitally induced forced precession of the NS in a NS -black hole (BH) binary (Yang & Zou 2020); forced precession produced by a fall-back disk (Tong et al. 2020); or, a tidal-force-induced precession of an emitting jet (Chen et al. 2021).Li et al. (2021b) show that the chromatic windowing in the activity of FRB 20180916B can be reconciled with a slowly rotating or freely precessing magnetar with a curvature radiation emission model (with radius-to-frequency mapping), with asymmetric emission about the magnetic dipole axis.The forced precession scenarios require the activity windows at all frequencies to be centered around the same precessional phase, rendering them implausible.Precession models can be tested by monitoring for a time derivative of the activity period (which would also be the precessional period; Wei et al. 2022).Our observed flat PPAs within each burst are consistent with measurements in the same frequency band in the past (P21) as well as at higher frequencies (Nimmo et al. 2021;Pastor-Marazuela et al. 2021), which are accommodated in both free and forced precession models by Li et al. (2021b).We do not calculate the absolute value of the PPAs, which if compared with those at higher frequencies would further constrain precession and slow rotation models that predict an evolution of the PPAs with activity phase.
The activity of the FRB 20180916B in its 16.33-day periodic window is frequency dependent, and we see in this work that the chromaticity remains stable with time in the LOFAR band (Section 3.1).However, this chromaticity does not seem to extend to the RM variations, with lower frequency LOFAR measurements of RM (of bursts B23 and B28) agreeing with the RM trend observed at higher frequencies by CHIME/FRB.In the context of binary models where the activity is modulated by the stellar wind or the disk of a companion, we infer that the RM variations must arise from a medium distinct from the one responsible for the chromatic activity windows.This strengthens the same conclusion drawn by Mckinven et al. (2023a) based on the fact that the RM variations occur on much longer timescales than the 16.3-day periodic activity.Such a magnetized screen that affects bursts occurring at different parts of the activity phase must exist between the source and us along our lineof-sight.This screen can be conceived to be a supernova remnant or a pulsar wind nebula (like the Crab nebula), although observations of the field around FRB 20180916B do not detect any persistent radio counterpart on milliarcsecond to arcsecond scales (Marcote et al. 2020).In case of an expanding supernova remnant, the |RM| would decrease (mostly) monotonically over much longer timescales (Yang et al. 2023) than the source has been currently observed for.

CONCLUSIONS
FRB 20180916B, with its 16.33-day periodic activity is one of the most prolific repeating FRBs.This allowed us to track its activity, as well as the observable properties of the bursts in the 110 − 188 MHz LOFAR band with high cadence, over years timescales.We have tracked frequency-dependent variations in the spectro-temporal and polarimetric properties of the burst, some of them caused by propagation effects in the intervening medium, and find no correlations with the chromatic activity of the source.Continued monitoring of the source with LOFAR and other telescopes will chart the trajectory of the observed RM evolution, which transitioned from being stochastic to a secular decrease.This will provide strong constraints on any self-consistent model of this prolific repeating FRB source.  .Difference between the measured scattering timescale at 150 MHz and the input scattering timescale, for simulated bursts with varying input drift rates and scattering timescales.See Section 3.2 for a full description of the simulated bursts.This is used to place systematic uncertainties on the scattering timescale fits for the detected bursts, accounting for covariance with the drift rate of the burst.

FrequencyFigure 1 .
Figure 1.Dynamic spectra (bottom panels of sub-figures), total intensity and polarimetric profiles (mid panels of sub-figures) and PPAs (top panels of sub-figures, for bursts with detected linear polarization) -for the 11 new bursts (B19−B29) detected from FRB 20180916B during the LOFAR campaign reported here.The PPAs for each burst have been rotated by an arbitrary angle such that they are centered around 0 deg).All bursts are dedispersed to a DM = 348.772pc cm −3 , the best-fit DM from Nimmo et al. (2021) who used burst micro-structure to determine this value.See Section 3.2 in the main text for a detailed justification of why we used this single DM value.The spectra have a time and frequency resolution of 7.864 ms and 0.781 MHz, respectively.All the time series profiles -total-intensity Stokes I (in black), linear polarization L (in red), and circular polarization Stokes V (in blue) -are obtained by summing only over the part of the band containing the burst.The white horizontal sections of the dynamic spectra represent masked RFI bands.The central frequency and 3- frequency extent of the bursts (upon fitting a Gaussian profile to the burst spectrum) are shown in the panels on the right of the dynamic spectrum in each sub-figure.
5 +0.058 −0.058 (CHIME/FRB Collaboration et al. 2020), lending support to the chromatic modulation of activity; (c) From MJD 58640.4 to MJD 59355, during the stochastic RM variations of 2 − 3 rad m −2 .The best fit window parameters are  LOFAR c = 0.68 +0.03 −0.02 and  = 4.44 +0.83 −0.33 ; (d) From MJD 59355 to MJD 59854, during the observed secular trend in RM in the 400 − 800-MHz CHIME/FRB band (Mckinven et al. 2023a, see Section 3.3).The best-fit window parameters are  LOFAR c Sand et al. (2022) constrain DM variations to < 1 pc cm −3 at 800 MHz for bursts in January 2020 using structure maximization, whileMckinven et al. (2023a) find < 0.8 pc cm −3 variations during a 9-month period in 2021 (that overlaps with our observing period).Using a different approach to find structure optimized DM, Lin et al. (2022) claim DM variations of ∼1.5 pc cm −3 .However, the optimized DM using their DM-power algorithm on the highest S/N burst in their sample agrees with the microsecond, structure-optimized DM fromNimmo et al. (2021) to within 0.005 pc cm −3 .There is a ninemonth difference between these two consistent DM measurements.Pastor-Marazuela et al. (2021) find a structure optimizing DM = 348.75± 0.12 pc cm −3 for APERTIF bursts at L-band, and find that this limits the rate of change of DM to less than 0.05 pc cm −3 per year, with bursts separated by up to ∼ 17 months from theNimmo et al. (2021) measurement.

Figure 3 .
Figure 3. Corner plot from the MCMC analysis showing constraints on the width of the LOFAR activity window , and the activity phase at which it is centered  LOFAR c, for different time ranges.These constraints are derived from accounting for the telescope exposure across all activity phases in the 16.33-day period.With ∼ 2.5 times the total exposure time since P21, and 49.8% of this spread outside the observed activity window to date, we obtain stronger constraints on the width of the window.In blue: Constraints on  and  LOFAR c using only the first 18 bursts, and exposures up until MJD 59088, as published in P21 (their Figure10).In orange: Constraints on  and  LOFAR c using all bursts and observations up to MJD 59854 (top-left sub-figure).In black: Constraints obtained for exposures before MJD 59355, at which the first burst with a ∼ 10% change in RM is observed.In red: Constraints obtained for exposures after (and including) MJD 59355, at which the first burst with a ∼ 10% fractional change in RM is observed.This shows that properties of the activity window at LOFAR have not shifted since the observed secular trend in RM began.

Figure 4 .
Figure 4. Variation of scattering timescales of the bursts with time and activity phase, obtained by fitting a Gaussian convolved with an exponential tail to the time profiles of the bursts.Only values for bursts with a S/N > 10 are marked in the plot.All scattering times have been referenced to 150 MHz, assuming   ∝  −4 .The error bars are a quadrature sum of the 1- statistical errors from the MCMC fitting, the errors derived from Figure A3 for the covariance between scattering and drift rate (as described in Section 3.2), and difference in scattering timescale if the values were scaled using   ∝  [−3.6≤ ≤ −4.4] when referenced to 150 MHz.The data points are coloured by their corresponding fitted drift rates.We include measurements for all detected FRB 20180916B LOFAR bursts B1−B29, i.e., including bursts reported in P21.The black dashed line represents the expected line-of-sight scattering from the Milky Way disk based on the NE2001 model, after accounting for plane-wave scattering of extragalactic bursts, and scaling the expected value at 1 GHz using   ∝  −4 .The gray shaded region shows the Milky Way ISM expected scattering contribution (NE2001 model) for frequency scaling indices in the range [ −3.6 ≤  ≤ −4.4].

Figure 6 .
Figure 6.RM versus time (Left) and versus activity phase (Right).Until B19, the RM variations were only ∼ 2 − 3 rad m −2 .The LOFAR RM measurements at 110 − 188 MHz agree with the RM trend observed by Mckinven et al. (2023a) at 400 − 800 MHz using CHIME/FRB.The LOFAR data points are also coloured by their corresponding linear polarization fractions, while unpolarized bursts are denoted by vertical dashed black lines at their respective arrival times and phases.Note that the linear polarization fractions vary with the central frequency of the bursts and with time (see Figure7).
Figure7.Top: Total intensity (black) and polarization profiles (linear: red; circular: blue) for a single burst (B1), at different radio-frequency ranges of the total spectral range occupied by the burst.The band occupied by the burst was divided into four quarters labelled Q1, Q2, Q3, Q4 from the lowest to highest frequency, respectively.The profiles for each quarter have been normalized by the peak of the total intensity, for visual purposes.The respective linear polarization fraction measured in these quarters can be found in the figure below.Bottom: Linear polarization fractions versus the central frequencies of the bursts.Additionally, B1 was split in bandwidth into quarters, and B3, B4, B19 into halves; and the linear polarization fractions of these parts plotted separately.We see that the linear polarization decreases with decreasing frequency even within the spectral extent of a single burst.Overplotted is the internal depolarization model (see Section 4.1 and Equation7).Additionally, we see that bursts B25, B28, and B29 are completely depolarized despite occupying the top-half of the frequency band.Upper limits are calculated for the bursts with no detected polarization (B21, B22, B24, B25, B28, and B29) at an interpolated RM value based on the best-fit slope for the RM secular trend fromMckinven et al. (2023a).

Figure A3
Figure A3.Difference between the measured scattering timescale at 150 MHz and the input scattering timescale, for simulated bursts with varying input drift rates and scattering timescales.See Section 3.2 for a full description of the simulated bursts.This is used to place systematic uncertainties on the scattering timescale fits for the detected bursts, accounting for covariance with the drift rate of the burst.

FrequencyFigure A4 .
Figure A4.Scattering timescales measured at four equipartitioned subbands for each burst, versus frequency.Both the axes are logarithmic.The fitted scattering indices,  (where scattering scales with frequency as  s ∝   ), for each burst are shown in the legend along with the 1- errors.

Figure A5 .
Figure A5.Top panel of each sub-figure:The MCMC-based fit of a Gaussian, convolved with an exponential tail model, to the frequency-integrated time profiles of bursts with S/N> 10, as described in Section 3.2.The data is at a time resolution of 3.93 ms.Bottom panel of each sub-figure: Difference between the fit model and measured burst profile in black (with the orange dotted line showing the line of zero difference).The fit for B14 might be overestimating the scattering timescale due to the presence of a smaller second peak in the time series profile of the burst, arising from a potential second burst component.Note that the bursts were dedispersed to a DM= 348.772 pc cm −3 see Section 3.2 in the main text for a detailed justification of why we used this single DM value.

Figure A6 .
Figure A6.Linear polarization fraction versus the scattering timescales of bursts at 150 MHz.The Pearson correlation coefficient between the two is −0.21, indicating a very weak correlation, if any.

Table 1 .
(Nimmo et al. 2021)e Section 3 for a description of how parameters were determined.For all bursts, arrival times and burst widths are computed for DM = 348.772pccm−3(Nimmo et al. 2021).Burst width (FWHM of Gaussian), for a Gaussian function with exponential tail fitted to the time series. a

Table 2 .
Measurements of propagation effects for the new bursts presented here.See also TableB1for a re-fitting of the bursts previously published in P21.
Burst  (at 150 MHz) a Drift rate (ACF method) b Drift rate (time of arrival method)