© ESO 2020

1. Introduction

Intra-hour variability (IHV) refers to quasi-random intensity modulation that is similar to canonical interstellar scintillation, but on vastly shorter timescales. It is an extreme manifestation of radio-wave propagation through turbulent plasma, but the nature of the plasma clouds causing the fast variability has remained unclear so far. It is generally thought to be caused by nearby (d ∼ 10 pc) turbulent plasma (Bignall et al. 2007; de Bruyn & Macquart 2015). Because the scattering strength of plasma is weaker if the plasma clouds are at small distances, the plasma densities inferred from IHV observations are orders of magnitude higher than that in the ambient ionised interstellar gas. The formation and survival of such over-pressured clouds is at the heart of the mystery. No multi-wavelength counterparts to the plasma have yet been found.

A recent development is the statistical association of three known IHVs with nearby hot stars by Walker et al. (2017) suggested that the scattering is caused by the ionised sheath of tiny self-gravitating molecular clumps similar to the cometary knots seen in the Helix nebula, for example. If true, this would be extremely interesting as these clouds would constitute a new population of objects within the interstellar medium. Alternatively, Pen & King (2012) and Pen & Levin (2014) have refined the model of Goldreich & Sridhar (2006), which suggests that the extreme scattering phenomena are caused by thin plasma sheets (which have a basis in magnetohydrodynamics, MHD, theory) oriented along the sightline. The implication here being that there is nothing inherently unusual about the intervening plasma save a fortuitous orientation.

Only very few IHVs are known. More such objects need to be discovered to better understand their statistical properties and for more secure empirical associations. More cases with extreme properties need to be studied to test the limits of proposed theories. The recent spate of wide-field survey radio telescopes operating in the GHz band, such as the Australian Square Kilometre Array Pathfinder (ASKAP) and the Aperture Tile In Focus (Apertif) project on the Westerbork Synthesis Radio Telescope (WSRT), are in an excellent position to provide larger samples. Here we report on the first IHV discovered by the recently started survey with the Apertif phased-array frontends on the WSRT.

2. Discovery and follow-up

We serendipitously discovered large and rapid scintillations in the radio source J1402+5347 (RA_J2000 14^h 02^m , Dec_J2000 53° 47′ 11″) with the Apertif phased-array feed system on the WSRT in commissioning data taken on April 9, 2019. Because of the east-west orientation of the WSRT, highly variable sources show prominent artefacts in the form of radial spikes in long-exposure images, which makes them easy to identify (Fig. 1). The light curve extracted for the source (Fig. 2; details in Appendix A) shows strong flux-density variations that are about 38% of the mean value. Archival WSRT observations from 2013 of the same field at the same frequency showed no significant variability in J1402+5347. No counterpart to J1402+5347 is visible in optical images, but an object is present at the location of J1402+5347 in images produced by the WISE satellite (Cutri et al. 2012). The lack of an optical counterpart together with the colours of this object in the WISE data suggests that the background object is a dust-obscured active galactic nucleus (AGN) at fairly high redshift (z ≳ 1; Jarrett et al. 2017).

Fig. 1. Discovery image of J1402+5347 (marked by the blue lines) made from the observation on April 9, 2019. The nearby spiral galaxy M 101 lies to the north. J1402+5347 was identified as a highly variable source because of the prominent radial spikes that are centred on it.

Fig. 2. Light curves of J1402+5347 observed over a year with a monthly cadence. The observation dates are indicated. The annual modulation in the variation rate is apparent. We observe two standstills in August 2019 and November 2019, while the most rapid variations occurred in April 2019 and March 2020.

Follow-up observations taken in May 2019 also showed strong variability at about 30%. Further observations in June and July 2019 showed weaker and slower variability (Fig. 2). An annual modulation of the scintillation properties has been observed in the other known IHVs (Dennett-Thorpe & de Bruyn 2003; Walker et al. 2009; Bignall et al. 2019) and is caused by the annual change in the relative velocity of the Earth with respect to the plasma clouds that is due to the orbital motion of the Earth around the Sun. Motivated by this, we continued monitoring the source with a cadence of about one month. The intensity and rate of the variations continued to decrease and J1402+5347 virtually ceased to vary in August. Following a brief period of slow variability in September and October, another standstill was observed in November. After November, the rate and intensity of the variations increased steadily, and by March 2020, they even exceeded the high levels seen in the discovery observation of April 2019.

3. Timescale analysis

The scintillation geometry and physical properties of the plasma clouds causing the scintillation can be determined from the annual variation in the scintillation timescale. We computed the scintillation timescale for each observation using two methods that are widely used in the analysis of time-series. The first method determines the decorrelation timescale from the noise-bias-corrected ACF directly computed from the light curves (Appendix B). The temporal lag at which the ACF of the variations falls below a specified threshold is taken as the scintillation timescale. For easy comparison with previous work (Dennett-Thorpe & de Bruyn 2003; Bignall et al. 2019), we define the threshold as 1/e times the peak value. This method is more accurate as it makes minimal assumptions. However, it is only effective for epochs of rapid variability as it requires several scintels to be observed within a single observation run (Dennett-Thorpe & de Bruyn 2003). For epochs near the standstills, we used a method that is commonly called Gaussian process (GP) regression (Appendix B). This method models the observed variations as a GP with a judiciously chosen functional form for the temporal covariance function. Although this method assumes an analytic form of the covariance function, it can effectively determine the decorrelation timescale and uncertainties even in epochs near the standstills (Bignall et al. 2019). The decorrelation timescales from both methods are tabulated in Table 1 and agree with each other within the errors for all the epochs to which both methods could be applied.

The scintillation timescale of just 6.5 min in March 2020 and the similarly short timescale in the April 2019 observations represent the fastest interstellar scintillation reported to date. This is remarkable because, compared to previous studies, our observations are at lower frequency for which normally the scintillations are significantly slower.

This timescale directly relates the physical size of the scintels to the transverse speed between the Earth and the scattering plasma. A single-epoch measurement cannot simultaneously constrain the scintel size and the relative speed. However, the annual modulation in the scintillation timescale can be used to estimate both parameters. In particular, the two standstills critically constrain both the screen velocity and the geometry of the scintels. One annual standstill can be obtained when the screen velocity on the ecliptic plane equals the Earth’s orbital velocity. The second standstill requires the scintels to be anisotropic such that the standstill occurs when the relative velocity between the Earth and the plasma cloud is along the major axis of the scintels. We fit a simple model with one-dimensional scintels to the observed annual modulation, as suggested by Walker et al. (2009, see Appendix C). The annual variation in this model is fully specified by three parameters: (a) the semi-minor axis of the scintels, a_⊥, (b) the systemic velocity of the scattering plasma along the minor axis of the scintels, v_⊥, and (c) the orientation of the major axis on the plane of the sky, θ_R (measured from north to east). Figures 3 and 4 show the posterior distributions of the model parameters and the best-fitting annual modulation curve, respectively. We constrain the semi-minor axis of the scintels to be just a_⊥ = 2.17(0.08)×10⁴ km. We note that the reduced chi-squared for the best-fitting model is close to unity, which implies that the one-dimensional scintel model sufficiently and parsimoniously captures the annual variation in the data. To understand the degree of anisotropy implied by the data, we also computed the posterior distribution for a two-dimensional scintillation model (Appendix C). The two-dimensional model places a lower limit on the ratio of the long and short axes of the scintels, R, of log₁₀(R) > 2.12. As expected from the high anisotropy, the data do not constrain the space velocity of the scattering screen parallel to the long axis of the scintels.

Fig. 3. Posterior distribution of the model parameters describing the annual modulation of the scintillation rate. The panels on the diagonal are the marginalised one-dimensional posterior distributions. The dashed black lines are drawn at 0.16, 0.5, and 0.84 quartiles. The annotated text gives the best-fit parameters and ±1σ bounds implied by the posterior distributions.

Fig. 4. Scintillation rate determined through a Gaussian regression analysis (blue points with ±1σ errors); the maximum likelihood model (black curve) is overplotted. The annotated text gives the model parameters and the reduced χ2 of the fit. The phase at the vernal equinox is defined to be zero, with phase increasing with time. Magenta crosses show the inverse of the scintillation timescale measured directly from the autocorrelation function (ACF) of the light curves during epochs of sufficiently rapid variations.

4. Discussion and outlook

Incoherent synchrotron sources are generally expected to only display refractive interstellar scintillation. Refractive scintels of a point-like source will have a size equal to or larger than that of the first Fresnel zone: , where λ is the wavelength and d is the distance to the scattering screen (Goodman & Narayan 2006). If the measured semi-minor axis of the scintels is associated with the Fresnel scale, the distance to the scattering plasma is d = 0.24 pc, which is within the sphere of gravitational influence of the Sun. This would place the scattering plasma within the canonically defined Oort cloud.

If J1402+5347 is unusually compact, then it may display diffractive interstellar scintillations, as is usually seen in pulsars. In this case, the scintel size is associated with the diffractive scale. Because the diffractive scale can be much smaller than the Fresnel scale, the data admit significantly larger screen distances in the diffractive regime. Regardless, the length scale of scintels cannot be smaller than the projected size of the source at the scattering screen. We consider for simplicity a circular source of angular radius θ_src. The requirement d θ_src <  a_⊥ places a constraint on the apparent brightness temperature of the source of T_b >  10¹¹(d/pc)² K.

Walker et al. (2017) have shown a statistical astrometric association between sightlines to sources showing intra-hour variability and hot stars such as Vega, Spica, and Alhakim. They also reported that the long axis of the scintels of IHVs preferentially points towards the associated star. We searched the HIPPARCOS catalogue to verify if the sightline to J1402+5347 passes anomalously close to any known O-, B-, or A-type star (see Fig. 5). The B3 star Alkaid (η UMa, HIP 67301) at a distance of d = 32 pc is unusually close to the sightline towards J1402+5347. The volume density of B-type stars in the HIPPARCOS catalogue that lie within a 50-degree-wide cone around this sightline is about 4.3 × 10⁻⁵ pc⁻³. The impact parameter of the sightline with Alkaid is about 2.7 pc and the angle between the scintel long axis and the star-ward direction for a pairing of J1402+5347 and Alkaid is about 9 ± 3 deg. The probability of finding a B-type star by chance with the observed impact parameter or closer and with the observed angular alignment or better is about 0.2%. In addition, the space velocity components of Alkaid with respect to the Solar System barycentre are −17.8 km s⁻¹ and −2.2 km s⁻¹ along the RA and Dec axes, respectively. This is about 9 km s⁻¹ different from the space velocity component of the scattering plasma as constrained by our observations, although, as pointed out by Walker et al. (2017), a clear agreement between space velocities may be precluded because of the free expansion of the scattering plasma, and to a lesser extent, because of projection effects. Regardless, an association of the scattering plasma with Alkaid sets d = 32 pc and places a rather extreme lower limit on the apparent brightness temperature of J1402+5347 of T_b >  10¹⁴ K. This limit can be slightly ameliorated if we postulate that the radio source itself is extended precisely along the major axis of the scintels.

Fig. 5. Left panel: line of sight and transverse distance of stars in the HIPPARCOS catalogue with respect to the sightline towards J1402+5347. The stars are colour-coded by their spectral type (legend in the annotated text). Right panel: same as the left panel, but the stars are plotted in the plane of the sky. The central black square shows J1402+5347. The solid and dashed blue lines mark the orientation of the major and minor axes of the scintels, respectively.

An association with Alkaid also has critical implications for the variance in the scattering properties of the Galactic plasma at large. The angular broadening expected from the line-of-sight integrated Galactic plasma (unrelated to the scattering screen discussed here) for this sightline is estimated to be about 0.36 mas according to the widely used NE2001 model (Cordes & Lazio 2002). The brightness temperature necessary for an association with Alkaid, however, requires that the source diameter does not exceed about 8.6 μas, which is 40 times smaller than the anticipated angular broadening. An association with Alkaid therefore suggests that the bulk of the scattering in the Galactic interstellar medium must occur in highly localised clumps so that this strong sightline-to-sightline variation can be accounted for.

Whether the scattering plasma is located in the Oort cloud or is associated with Alkaid may be determined by extending theoretical analyses, such as the analysis of Goodman & Narayan (2006), to the case of highly anisotropic scintels, and by broadband observations of J1402+5347 at higher frequencies. If the observed scattering at 1.4 GHz is indeed diffractive, two effects must be seen in such observations: (a) the scintillations will de-correlate on frequency scales Δν/ν ≪ 1, and (b) the transition to weak scattering must occur at frequencies that are several times higher than 1.4 GHz, and similarly strong light-curve variations must then extend to these high frequencies. If the scattering is associated with the Oort cloud, then the observed variations are due to refractive scintillations in the weak-scattering regime. We therefore expect to observe broadband variations (Δν/ν ∼ 1) at higher frequencies with a continuously declining level of fractional variations with increasing frequency. The observations of April 9, 2019, show that the variations at the upper and lower end of the Δν/ν ≈ 15% observing band are very similar, with a Pearson correlation coefficient of 0.84 (Appendix D), suggesting refractive variations and hence that the scattering plasma is located in the Oort cloud. However, further theoretical investigation on the spectral decorrelation properties of highly anisotropic scintels are required to clearly determine the distance to the scattering screen.

We end by noting that because of the larger apparent source sizes at lower frequencies, fast intra-hour variability at decimetre wavelengths will be significantly biased towards intervening plasma clouds that are nearby. The discovery of rapid scintillation in J1402+5347 in early Apertif data bodes well for ongoing wide-field surveys in the decimetre band. Based on the discovery of J1402+5347 and the sky area covered by the Apertif survey at the time, we expect to discover approximately ten more such systems, and indeed, several new IHV sources have more recently been identified in Apertif data (Oosterloo et al., in prep.).

Acknowledgments

We dedicate this Letter to the memory of Ger de Bruyn and of J.P. Macquart. HKV thanks Mark Walker for discussions. Credit is also due to Mark Walker for pointing out the possible association with Alkaid. We thank the anonymous referee for a meticulous reading of the manuscript which helped us to improve our paper, and for rectifying some errors in notation. EAKA is supported by the WISE research programme, which is financed by the Netherlands Organisation for Scientific Research (NWO). YM acknowledges funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 617199. LCO acknowledges funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 617199. JMvdH acknowledges funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 291531 (“HIStoryNU”), JvL acknowledges funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 617199 (“ALERT”), and from Vici research programme “ARGO” with project number 639.043.815, financed by the Netherlands Organisation for Scientific Research (NWO). This work makes use of data from the Apertif system installed at the Westerbork Synthesis Radio Telescope owned by ASTRON. ASTRON, the Netherlands Institute for Radio Astronomy, is an institute of the Dutch Science Organisation (De Nederlandse Organisatie voor Wetenschappelijk Onderzoek, NWO).

References

Appendix A: Interferometric data reduction

The observations of J1402+5347 presented here were performed with the Apertif phased-array feed system that was recently installed on the WSRT. With Apertif, 40 partially overlapping beams, each with a full width at half maximum of about 35 arcmin, are formed on the sky. The data from each beam are used independently to image the region covered by each beam using standard aperture synthesis. The total area imaged by combining the images from all beams is about with a spatial resolution of 12″ × 17″ at the declination of J1402+5347. The observing band was 122 MHz wide, divided into 156 channels, centred on 1.365 GHz. J1402+5347 was detected in three overlapping beams. Each observation had a duration of 11 to 11.5 h. The calibration of the data from each beam followed standard procedures, involving excision of radio-frequency interference, flux-, bandpass-, and cross-calibration using 3C147, followed by self-calibration. The resulting images had small direction-dependent imaging errors, which were corrected for using standard peeling techniques. This applies in particular to the strong source 3C295, located SW of J1402+5347: because of the large distance from the beam centres, it showed large direction-dependent effects. The resulting images are noise limited with a noise level of about 30 μJy beam⁻¹ that varies slightly between observations.

To construct the light curves, the source model resulting from the self-calibration was subtracted from the visibilities, except for the source model for J1402+5347. This created a visibility data set that contained only emission from J1402+5347. The phase centre of this visibility data set was shifted to the position of J1402+5347 so that the phases of all visibilities were zero. Finally, for each time stamp, the visibilities of all baselines and frequencies were summed to obtain the flux density of J1402+5347 for each time stamp. After primary beam correction, this yields the light curves shown in Fig. 2.

To verify the accuracy of the light curves, the same procedure was followed for J1403+5349, which has a very similar flux density, 9 arcminutes W of J1402+5347. An example of such a comparison is shown in Fig. A.1.

Fig. A.1. Light curves derived for April 9, 2019, of the variable source J1402+5347 and the steady comparison source J1403+5349, 9 arcmin west of J1402+5347, illustrating the accuracy of the method for constructing the light curves. The horizontal axis shows the time since the start of the observation.

Appendix B: Determination of the scintillation timescale

We computed the scintillation timescale using the autocorrelation function (ACF) and Gaussian process modelling (GP). For the estimates based on the ACF, we followed a procedure that closely mimics that presented by Bignall et al. (2003). We first subtract the mean flux density from the light curves to obtain a time series, x_i, with a sampling interval Δt. For each pair ij we compute the product x_ix_j. The list of products is then binned on a grid of size Δt, and all products falling into the same bin were averaged. Thermal noise only biases the zeroth bin as it is uncorrelated between light-curve samples. This bias, equal to the variance of the noise, is estimated to be equal to half the variance of the difference between successive samples in the time series. The 1/e timescale, t_e, is estimated using a straight-line fit to a segment of the ACF containing five samples, centred on the sample that is closest to 1/e. The error on the scintillation timescale is given by , where N is the number of observed cycles of variation on timescale t_e and ϵ is a factor of about unity that depends on the exact profile of the ACF around t_e. If the total duration of observation is T (∼11.5 h in our case), then N = (T/t_e). The factor ϵ was estimated with numerical Monte Carlo techniques (Dennett-Thorpe & de Bruyn 2003; Bignall et al. 2003) to be in the range 0.25 − 0.9. We conservatively assumed ϵ = 1. The timescale and its errors thus computed are only accurate when several scintels are observed, which is not the case during the standstills. An example ACF is plotted in Fig. B.1.

Fig. B.1. Numerical ACF computation: mean subtracted light curve (top panel) and the ACF (bottom panel) for the April 9, 2019, light curve.

Our second approach was to model the temporal scintillation as a GP with a covariance function given by a damped exponential: K(Δt) = K₀exp(cΔt)cos(dΔt), where Δt is the temporal lag and K₀, c, and d are model parameters (Bignall et al. 2019). This choice for the covariance function is motivated by the shape of the ACF, which shows a steep roll-off at small lags, follows by oscillatory behaviour due to the dominant spectral mode in the light curves (see Fig. B.3) We assume that the scintillation is a stationary process, which means that the parameters remain the same at different epochs, and the variation in the scintillation timescale can be modelled by a single epoch-dependent scale factor, k, by substituting Δt → kΔt. This allows us to estimate the scintillation timescale and the associated uncertainty even for epochs close to the standstill by employing our knowledge of the covariance function from epochs of rapid variability. We determined the parameters K₀, c, and d from the April 2019 light curve alone instead of using a fit using the entire annual data set (Bignall et al. 2019). With these parameter values fixed, we determined the temporal scale factor k for each epoch with a single-parameter Gaussian regression.

All data were mean subtracted prior to regression analysis, and a noise variance term (estimated from the time differences) was added to the diagonal elements of the covariance matrix. For the April 2019 data set, we fixed K₀ to equal the total variance in the scintillations: K₀ = 3.477105 mJy², and we estimated c and d. The best-fit estimates were log(c) = − 3.80 ± 0.10 and log(d) = − 1.86 ± 0.05. The posterior distributions for these parameters were computed using the emcee software (Foreman-Mackey et al. 2013) and are plotted in Fig. B.2. The numerically computed ACF and the best-fit kernel are plotted together for comparison in Fig. B.3.

Fig. B.2. Best-fitting kernel: numerically computed ACF overplotted with the best-fit Gaussian kernel for the April 9, 2019, light curve. The decorrelation timescale of the Gaussian kernel is indicated with the dashed black line. The quasi-sinusoidal oscillations at long delay are not ideally captured by the analytic form of the kernel, However, the decay to the 1/e-value at short delays is well captured and is most relevant to estimating scintillation rates.

Fig. B.3. Gaussian process kernel parameters: posterior distribution of GP kernel parameters for the April 9, 2019, data set. The total variance K0 was fixed to equal the noise-bias corrected variance of the light curve, and only the scale parameters c and d were estimated. The dashed lines are placed at the 16, 50, and 84th quantiles of the marginalised one-dimensional distributions. The fractional error in the estimation of parameter c is about 10%.

Next, we held K₀, c and d fixed at the best-fit values quoted above and used the celerite (Foreman-Mackey et al. 2017) package to fit a single scale parameter k. The 1/e crossing timescale for the April 2019 epoch is largely set by the value of parameter c, which is determined with a formal error of about 10%. To account for this uncertainty in our error budget, we added a 10% error in quadrature to the formal errors on the timescale from the GP regression analysis.

Throughout our analysis, the likelihood computations were made using routines from the celerite package (Foreman-Mackey et al. 2017), the best-fit values were determined using scipy.optimize with the L-BFGS-B method, and the posterior distributions were sampled using the emcee package (Foreman-Mackey et al. 2013).

Appendix C: Annual modulation

The timescale of a temporal scintillation caused by anisotropic scattering is given by (Walker et al. 2009)

(C.1)

where R = a_||/a_⊥ is the anisotropy of the scintels, a_⊥ and a_|| are the length scales of scintels along the minor and major axes, respectively, and v_⊥ rel and v_|| rel are the relative velocity between the Earth and the scintels along the minor and major axes of the scintels, respectively. For highly anisotropic scintels (Walker et al. 2009), we can set R → ∞ and obtain t_s = a_⊥/v_⊥ rel. The annual modulation profile can be fit with a three-parameter model. The parameters are (i) the semi-minor axis of the scintels, a_⊥, (ii) the orientation of the major axis of the scintels in the plane of the sky, θ_R ∈ (0, 2π], measured from north towards east, and (iii) the screen velocity in the plane of the sky perpendicular to the major axis of the scintels v_⊥ ≥ 0 (defined to be along θ_R − π/2). We carried out our fits in a tangent plane perpendicular to the sightline with cardinal axes along the directions of increasing RA () and Dec (). If is the velocity of the Earth projected onto this plane, then the relative velocity perpendicular to the major axis of the scintels is given by . We fit this equation to the scintillation timescales tabulated in Table 1 (column “GP scintillation rate”). We carried out the fits to the scintillation rate instead of the scintillation timescale to avoid infinities in the model. First, we carried out a brute-force search for the maximum likelihood value on a coarse three-dimensional parameter grid spanned by 10⁸ cm < a_⊥ <  5 × 10¹⁰ cm, 0 km s⁻¹ < v_⊥ < 100 km s⁻¹, and 0 <  θ_R <  2π. The likelihood peaked at about a_⊥ ≈ 2 × 10⁹ cm, v_⊥ ≈ 22 km s⁻¹, and θ_R ≈ 0.7. We then used the emcee software (Foreman-Mackey et al. 2013) to obtain the final fit values and their formal errors. We used 200 randomly initialised walkers, each walking 2000 steps. We imposed a flat prior on the parameters in the range 2 × 10⁸ cm <  a_⊥ <  2 × 10¹⁰ cm, 15 km s⁻¹ < v_{⊥, R} < 35 km s⁻¹, and 0.1 <  θ_R <  1.1. The posterior parameter distributions were computed and plotted from the emcee samples using the corner package.

We also ran a separate MCMC simulation to sample the posterior distribution for the two-dimensional scintillation model given in Eq. (C.1). The additional parameters in the two-dimensional model and v_|| and R. We placed the prior bounds v_|| ∈ [0, 100) km s⁻¹ and log₁₀R ∈ [0, 4), with the relative velocity parallel to the scintel long axis given by . The posterior distributions are given in Fig. C.1.

Fig. C.1. Same as Fig. 3, but for a two-dimensional scintel model given in Eq. (C.1).

Appendix D: Broad-band nature of scintels

To investigate a possible frequency dependence of the intensity fluctuations, we constructed separate light curves for the upper and lower halves of the observing band for the observations on April 9, 2019. These separate light curves are shown in Fig. D.1, where the blue line represents the light curve for the frequency range 1311−1408 MHz and the red line shows the range 1408−1505 MHz. It is clear from this figure that the light curves are very similar, showing that there is no frequency structure with the observing band. The Pearson correlation coefficient between the two light curves is 0.84.

Fig. D.1. Separate light curves for the lower and upper halves of the observing band of the discovery observation on April 9, 2019. The blue line is the light curve for the frequency range 1311−1408 MHz, and the red line shows the range 1408−1505 MHz.

All Tables

All Figures

	Fig. 1. Discovery image of J1402+5347 (marked by the blue lines) made from the observation on April 9, 2019. The nearby spiral galaxy M 101 lies to the north. J1402+5347 was identified as a highly variable source because of the prominent radial spikes that are centred on it.
In the text

	Fig. 2. Light curves of J1402+5347 observed over a year with a monthly cadence. The observation dates are indicated. The annual modulation in the variation rate is apparent. We observe two standstills in August 2019 and November 2019, while the most rapid variations occurred in April 2019 and March 2020.
In the text

	Fig. 3. Posterior distribution of the model parameters describing the annual modulation of the scintillation rate. The panels on the diagonal are the marginalised one-dimensional posterior distributions. The dashed black lines are drawn at 0.16, 0.5, and 0.84 quartiles. The annotated text gives the best-fit parameters and ±1σ bounds implied by the posterior distributions.
In the text

Fig. 4. Scintillation rate determined through a Gaussian regression analysis (blue points with ±1σ errors); the maximum likelihood model (black curve) is overplotted. The annotated text gives the model parameters and the reduced χ2 of the fit. The phase at the vernal equinox is defined to be zero, with phase increasing with time. Magenta crosses show the inverse of the scintillation timescale measured directly from the autocorrelation function (ACF) of the light curves during epochs of sufficiently rapid variations.

In the text

Fig. 5. Left panel: line of sight and transverse distance of stars in the HIPPARCOS catalogue with respect to the sightline towards J1402+5347. The stars are colour-coded by their spectral type (legend in the annotated text). Right panel: same as the left panel, but the stars are plotted in the plane of the sky. The central black square shows J1402+5347. The solid and dashed blue lines mark the orientation of the major and minor axes of the scintels, respectively.

In the text

	Fig. A.1. Light curves derived for April 9, 2019, of the variable source J1402+5347 and the steady comparison source J1403+5349, 9 arcmin west of J1402+5347, illustrating the accuracy of the method for constructing the light curves. The horizontal axis shows the time since the start of the observation.
In the text

	Fig. B.1. Numerical ACF computation: mean subtracted light curve (top panel) and the ACF (bottom panel) for the April 9, 2019, light curve.
In the text

Fig. B.2. Best-fitting kernel: numerically computed ACF overplotted with the best-fit Gaussian kernel for the April 9, 2019, light curve. The decorrelation timescale of the Gaussian kernel is indicated with the dashed black line. The quasi-sinusoidal oscillations at long delay are not ideally captured by the analytic form of the kernel, However, the decay to the 1/e-value at short delays is well captured and is most relevant to estimating scintillation rates.

In the text

Fig. B.3. Gaussian process kernel parameters: posterior distribution of GP kernel parameters for the April 9, 2019, data set. The total variance K0 was fixed to equal the noise-bias corrected variance of the light curve, and only the scale parameters c and d were estimated. The dashed lines are placed at the 16, 50, and 84th quantiles of the marginalised one-dimensional distributions. The fractional error in the estimation of parameter c is about 10%.

In the text

	Fig. C.1. Same as Fig. 3, but for a two-dimensional scintel model given in Eq. (C.1).
In the text

	Fig. D.1. Separate light curves for the lower and upper halves of the observing band of the discovery observation on April 9, 2019. The blue line is the light curve for the frequency range 1311−1408 MHz, and the red line shows the range 1408−1505 MHz.
In the text