Long-term photometric monitoring and spectroscopy of the white dwarf pulsar AR Scorpii Purchased

ABSTRACT

1 INTRODUCTION

The binary system AR Scorpii (henceforth AR Sco) is composed of a white dwarf with an M-dwarf companion in a 3.56 h orbit (Marsh et al. 2016). In contrast to cataclysmic variable systems (CVs), which show this same binary configuration, there is no evidence that the white dwarf in AR Sco accretes mass from its companion. It is further distinguished from other similar binaries by strongly pulsed emission with a period of 1.97 min (see Fig. 1) over a broad range of wavelengths from radio (Stanway et al. 2018) to X-rays (Takata et al. 2018), which led to it being dubbed a white dwarf pulsar (Geng, Zhang & Huang 2016; Buckley et al. 2017). The physical process behind AR Sco’s pulses is, nonetheless, not the same as for traditional neutron-star radio pulsars (Katz 2017; Lyutikov et al. 2020). This is inferred from the fact that no phase coherent radio emission has been reported for AR Sco, and from the companion being within the radius of the white dwarf’s light cylinder.

Figure 1.

Light curve obtained with HiPERCAM on 2018 April 15 (see Section 2.1). HiPERCAM filters u, g, r, i, z are shown from top to bottom. Strong pulses can be seen on all bands on top of the orbital modulation.

Open in new tabDownload slide

The pulsed variability of AR Sco was first reported by Marsh et al. (2016), who detected two components of very similar frequency in AR Sco’s amplitude spectra at ≈731 and ≈738 cycles per day (c/d). The difference between the two values corresponds to the orbital frequency, suggesting that the higher frequency component represents the spin of the white dwarf, whereas the lower frequency is a re-processed beat frequency. Marsh et al. (2016) also found the maximum luminosity of AR Sco to be well in excess of the combined luminosity of the stellar components. Given the lack of accretion signatures (such as strong X-ray luminosity, broad emission lines, or flickering), they proposed this excess luminosity to be powered by the spin-down of the rapidly-rotating white dwarf. They estimated a spin frequency derivative of |$\dot{\nu }_\mathrm{spin} = (-2.86\pm 0.36) \times 10^{-17}$| Hz s⁻¹, corresponding to a rate of loss of rotational energy of |$L_{\dot{\nu }} = 1.5 \times 10^{26}$| W.

Potter & Buckley (2018a) carried out a Fourier analysis of data spanning three years (2015–2017) but could not confirm the spin-down found by Marsh et al. (2016). Instead they found a spin frequency derivative consistent with zero, but were only able to constrain it to the broad range |$-20 \times 10^{-17}\,$| to +10 × 10⁻¹⁷ Hz s⁻¹. On the other hand, Stiller et al. (2018) confirmed the occurrence of spin-down, again from just three years coverage (2016 to 2018), but found that it had a significantly larger magnitude, |$\dot{\nu }_\mathrm{spin} = (-5.14\pm 0.32) \times 10^{-17}$| Hz s⁻¹, than determined by Marsh et al. (2016). Stiller et al. (2018) argued that the apparent inconsistency between the results of Marsh et al. (2016) and Potter & Buckley (2018a) arose from an underestimate of |$\dot{\nu }_\mathrm{spin}$| in Marsh et al. (2016), likely due to sparse sampling and low-time resolution of their data. They also pointed out that their value of |$\dot{\nu }_\mathrm{spin}$| could explain the difference in the spin periods found by Marsh et al. (2016) and Potter & Buckley (2018a) being consistent with the spin period found in both works. More recently, the spin-down rate was further refined to |$\dot{\nu }_\mathrm{spin} = (-4.82\pm 0.18) \times 10^{-17}$| Hz s⁻¹by Gaibor et al. (2020). The spin-down rate measurements have thus converged with the most recent consistent with each other, but questions remain over the earlier measurements. Marsh et al. (2016) based their measurement upon data from the CRTS survey that preceded (by several years) the more recent intensive follow-up data used by Potter & Buckley (2018a), Stiller et al. (2018), and Gaibor et al. (2020). The CRTS data were very sparse, and used relatively long exposures (30 s), which rightly raises a question mark over the reliability of Marsh et al. (2016)’s measurement of spin-down (Potter & Buckley 2018a). However, could Marsh et al. (2016)’s result together with those of Stiller et al. (2018) and Gaibor et al. (2020) point towards a genuine change in the rate of spin-down in AR Sco? Such behaviour is often seen amongst close relatives of AR Sco, the intermediate polars (IPs, see e.g. fig. 2 of Patterson et al. 2020). If the same applies to AR Sco, it could suggest that (fluctuating) mass accretion has some role to play in the spin-down, despite a relative lack of any other evidence for accretion in the system (Marsh et al. 2016). A second puzzle is why Potter & Buckley (2018b) obtained so loose a constraint upon the spin derivative, some 30 times worse than Stiller et al. (2018) from data covering the same time span. Is the Fourier-based method somehow flawed in the case of AR Sco?

The orbital period of the system is also expected to change as a result of angular momentum losses caused by the interplay between magnetic braking and gravitational radiation. Stiller et al. (2018) constrained the orbital frequency derivative to |$\dot{\nu }_\mathrm{orbit} \lesssim 2 \times 10^{-18}$| Hz s⁻¹, which at their precision implied |$\dot{\nu }_\mathrm{spin} = \dot{\nu }_\mathrm{beat}$|⁠. An even tighter constraint of |$\dot{\nu }_\mathrm{orbit} \lesssim 3.8 \times 10^{-20}$| Hz s⁻¹ was found by Peterson, Littlefield & Garnavich (2019) using archival photometric plate data dating back to 1902. They pointed out that their constraint is consistent with the semi-empirical model of Knigge, Baraffe & Patterson (2011) describing angular momentum loss in CVs, which predicts an orbital frequency derivative of |$\dot{\nu }_\mathrm{orbit} \lesssim 2 \times 10^{-20}$| Hz s⁻¹for AR Sco’s orbital period. Additional studies of the orbital behaviour were performed by Littlefield et al. (2017), who found that the light curve shows a stable orbital waveform on time-scales of days, but alters slowly over time-scales on the order of years.

As well as time-series photometry, time-resolved spectra have also been obtained for AR Sco. The spectra in the visible show a blue continuum added to the red spectrum of an M5-type dwarf with strong Balmer and neutral helium lines in emission (Marsh et al. 2016). The radial velocity variation of the emission lines suggests that they originate near the surface of the M-dwarf facing the white dwarf. Garnavich et al. (2019) found the Balmer lines to show significant structure as they vary over the orbital period. Similar to what is observed for the orbital modulation of the light curve, the Balmer lines show peak flux between phases 0.4 and 0.5, and minima around phase 0.1.¹ The continuum shows a similar behaviour, but is further affected by fast variations on the beat period. Doppler tomography shows that most of the Hα emission originates on the face of the secondary star, but satellite emission features are seen, which in other systems have been attributed to long-lived prominences on the secondary star (e.g. Parsons et al. 2016). Garnavich et al. (2019) argued that the stability of these prominences requires the M-dwarf to have a magnetic field of the order of 100–500 G (Ferreira 2000), which is comparable to the white dwarf magnetic field strength at the M-dwarf (∼200 G, Takata et al. 2018), implying that interaction between the two stars’ magnetic fields is likely responsible for the energy production in AR Sco.

A pulsar model for AR Sco is supported by polarimetric observations, which detected strong linear polarization (as high as 40 per cent) modulated on both the spin and beat periods (Buckley et al. 2017). The observed broad-band spectrum is associated with synchroton radiation from relativistic electrons accelerated as the white dwarf magnetosphere sweeps past the M-dwarf, resembling the emission from neutron star pulsars, but with a different mechanism whose details remain a mystery. Most likely it relates to direct interaction between the two stars, resulting in emission from the surface or coronal loops of the M-dwarf from the magnetosphere of the white dwarf, or possibly through an associated bow shock (Geng et al. 2016; Katz 2017; Takata, Yang & Cheng 2017; Potter & Buckley 2018b; du Plessis et al. 2022). A collimated relativistic outflow, i.e. a jet has also been considered, but seems to be ruled out by radio observations (Marcote et al. 2017).

AR Sco’s nature, and the rapid spin-down of its white dwarf suggests that it could be a missing link between the two major classes of accreting magnetic white dwarfs, polars, and IPs, which are distinguished by whether or not the white dwarf’s spin is synchronous with the orbit of the binary (Katz 2017), but how AR Sco attained its current configuration is unknown. A prior stage of accretion-driven spin-up is required to explain the current spin period. However, the high-magnetic field estimated for the white dwarf in AR Sco given its measured spin-down (100–500 MG, Katz 2017; Buckley et al. 2017) presents a significant challenge as an extremely high-mass transfer rate would have been required in the past to crush the magnetosphere down to a co-rotation radius matching the 2 min spin-period. In an alternative model, Lyutikov et al. (2020) argued the magnetic field of AR Sco must be only ∼10 MG to allow its spin-up, and that this is consistent with the observed spin-down rate if the material leaving the M-dwarf is originally neutral. In that case, the material can approach the white dwarf without interacting until it is ionized by the white dwarf itself. This model requires both that the white dwarf temperature be such that the material is ionized only near the white dwarf (⁠|$T_\mathrm{eff} \approx 12\, 000{\rm K}$|⁠), and that the heated face of the M-dwarf is cool enough not to ionize the gas before it heads towards the white dwarf. Although the first condition could be fulfilled (Garnavich et al. 2021), the heated face of the M-dwarf seems to be hot enough that the material would have a high-ionization rate (Garnavich et al. 2019, 2021).

Another potential solution that reconciles the need for accretion-driven spin-up and the apparently high-magnetic field was put forward by Schreiber et al. (2021). They proposed that AR Sco only became magnetic as a result of a crystallization- and rotation-driven dynamo. In this model, AR Sco would originally be non-magnetic allowing for straightforward accretion-driven spin-up. When crystallization starts to occur in the core of the cooling white dwarf strong density stratification combined with convection creates the conditions for a dynamo, generating the magnetic field (Isern et al. 2017). If the field is strong enough, the rapid transfer of spin angular momentum into the orbit may cause the binary to detach and mass transfer to cease, leading to a system such as AR Sco.

The many open issues surrounding the modelling of AR Sco, as well as its potential central role for constraining the evolution of CVs, make precise observational constraints a valuable asset. In this work, we analyse high-speed photometry of AR Sco obtained with ULTRACAM (Dhillon et al. 2007), ULTRASPEC (Dhillon et al. 2014), and HiPERCAM (Dhillon et al. 2021) spanning seven years (2015–2022), and time-resolved spectra obtained with the Intermediate-dispersion Spectrograph and Imaging System² (ISIS) at the William Herschel Telescope (WHT) and with X-shooter at the Very Large Telescope (VLT). We further constrain the spin-down of the white dwarf using two different techniques, pulse-time arrival modelling, and Fourier analysis. We probe for evidence of orbital period change by modelling spin and beat frequencies independently in the Fourier analysis. Based on our photometry, we discuss the possibility of using AR Sco as a cosmic clock for timing observations. Using our spectroscopic data, we perform modulation Doppler tomography (Steeghs 2003), extending the work of Garnavich et al. (2019). We also measure the rotational velocity of the M-dwarf as function of orbital phase to probe for variability that could confirm that the M-dwarf is Roche lobe filling.

2 OBSERVATIONS

2.1 Photometry

We started a monitoring campaign of AR Sco shortly after its discovery in 2015. We used three different high-speed photometers: ULTRACAM at the 3.5-m ESO New Technology Telescope (NTT), ULTRASPEC mounted at the 2.4-m Thai National Telescope (TNT), and HiPERCAM at the 10.4-m Gran Telescopio Canarias (GTC) (Dhillon et al. 2007, 2014, 2021). We observed it on a total of 68 different nights over a time span of almost seven years (2015–2022). Our cadence varied between 0.1 and 9.1 s with an average of 3.3 s. This includes the readout time, which varied with the configuration and readout mode, but was never longer than 15 ms. With ULTRASPEC, we typically used a g filter, with the exception of three nights when a red-blocking KG5 filter (effectively u + g + r) or an r filter was installed. For ULTRACAM and HiPERCAM dichroic beam splitters allow simultaneous observations in three or five different filters, respectively. Filters u, g, r or u, g, i were used for ULTRACAM until 2017 May 14, after which they were replaced with similar filters whose cut-on/off wavelengths match those of the typical Sloan Digital Sky Survey (SDSS) filters, but with a higher throughput, being sometimes referred to as u_s, g_s, r_s, i_s, z_s or ‘super’ SDSS. These ‘super’ filters were also used for the HiPERCAM observations. For simplicity, we omit the s subscript in the text, but details of each of our runs and filters used can be found in Table A1.

All data sets were reduced with the dedicated HiPERCAM data reduction pipeline.³ We performed bias subtraction and flat-field correction using skyflats taken during twilight. We also used the pipeline to carry out differential aperture photometry with a variable aperture size set to scale with the seeing estimated from a point-spread function (PSF) fit. We used the same comparison star for all observations with the exception of our two HiPERCAM runs for which the chosen comparison star was not in the field of view. Our default comparison star was Gaia EDR3 6050299715251225344 (G = 13.3, G_BP − G_RP = 2.3), whereas for HiPERCAM we used Gaia EDR3 6050298302203861888 (G = 15.2, G_BP − G_RP = 1.9). For reference, AR Sco itself has reported values of G = 15.0 and G_BP − G_RP = 1.4.

To extend our baseline to earlier years, we retrieved data for AR Sco from the Catalina Real-Time Transient Survey (CRTS, Drake et al. 2009). These data span from 2005 August 01 to 2013 July 08. The typical CRTS exposure time of 30 s is short enough that the spin and beat periods are not completely smeared out.

The observing times of all data sets were corrected to Barycentric Julian Date (BJD) in the Barycentric Dynamical Time (TDB) reference system.

2.2 Spectroscopy

We obtained spectroscopic observations of AR Sco in 2016 both with X-shooter at the 8.2-m VLT (Vernet et al. 2011) and with ISIS at the 4.2-m WHT.

The X-shooter observations were carried out using a 1 arcsec slit for the UVB arm (300–559.5 nm, R = 5400), and 0.9 arcsec slits for the VIS (559.5–1024 nm, R = 8900) and NIR (1024–2480 nm, R = 5600) arms. In the UVB arm, we limited our exposure time to 16.1 s to allow for sampling of the beat/spin variability and obtained 689 exposures (readout time was 14 s). A fraction of the spectra that were taken continuously (the observations were briefly interrupted by clouds) is shown in Fig. 2. For the VIS arm, we obtained 42 exposures of 354 s (≃3 times the beat period; readout time was 58 s). The same exposure time was used for the NIR arm, but the difference in the way readout is carried out allowed 47 frames to be taken (with a readout time of 8.2 s). The spectra were reduced using the standard procedures within the ESO reflex reduction pipeline,⁴ including telluric line removal using molecfit (Kausch et al. 2015; Smette et al. 2015).

Figure 2.

Trailed spectra taken with the X-shooter UVB arm. The spectra have been normalized by the continuum prior to plotting. Yellow indicates stronger flux. Many emission lines, which generally trace the heated face of the M-dwarf can be seen.

Open in new tabDownload slide

ISIS observations were obtained using the R600B and R600R gratings in the blue and red arms, respectively. The central wavelength was set to 4500 Å for the blue arm, and to 6200 Å for the red arm. We used a 1.0 arcsec slit and set the exposure time to 7.3 s (blue) and 8.0 s (red), whereas the readout time was ≈70 s. 1078 and 1081 consecutive exposures were obtained for blue and red arms, respectively. All spectra were de-biased and flat-fielded using the starlink⁵ packages kappa, figaro, and convert. Optimal spectral extraction was carried out using pamela and wavelength calibration was performed with molly (Marsh 1989). Spectra taken with the red arm are shown as a trail in Fig. 3.

Figure 3.

Trailed spectra taken with the ISIS red arm. The spectra have been normalized by the continuum prior to plotting. The panel to the right zooms into Hα, which shows some substructures that will be discussed in Section 3.3.

Open in new tabDownload slide

3 DATA ANALYSIS

3.1 Pulse arrival time estimates

We have estimated individual pulse arrival times for each of our observing runs with ULTRASPEC, ULTRACAM, or HiPERCAM. As the pulse timing depends on the filter used (see Gaibor et al. 2020 and Section 4.2), we have only used data taken with the g filter for this purpose. Our method is similar to that previously employed by Stiller et al. (2018) and Gaibor et al. (2020), who used a Gaussian function to model the beat pulse peak. We first estimated the location of each peak using previously determined ephemeris for the beat period. Next, we cross-correlated the data around each peak with a Gaussian function with a fixed width set to minimize the residuals when fitting the ephemeris. We experimented with standard deviation (σ) values for the Gaussian ranging from 5 to 20 s, finding an optimal value of 10.1 s. To determine the maximum of the cross-correlation function, which is the estimate for the location of the peak, we determined the value at which the derivative of the cross-correlation changes sign from positive to negative using Newton–Raphson iteration. The initial ephemeris was also used to estimate the cycle number of each peak measurement. We repeated the procedure recalculating the cycle numbers using the derived ephemeris in order to get self-consistent values; this was repeated until the assumed and fitted ephemeris showed no significant change.

Data points affected by clouds were manually flagged in each run, and for each peak we computed the resulting number of bad points. Peak measurements were only accepted if the number of good points exceeded the number of bad points by at least a factor of three. Peaks were also rejected if the difference between the measurement and predicted value based on the approximate ephemeris was more than 10 per cent of the beat period. This method resulted in 3195 good pulse measurements (table included as Supplementary Material). However, the beat pulses are known to show an orbital phase dependence (Stiller et al. 2018), which needs to be corrected for. To account for this orbital modulation, we performed an initial fit to the derived pulse times, and calculated the observed minus calculated values (O–C) as a function of orbital phase ϕ_orb. We then modelled the O–C behaviour with a Fourier series, and subtracted the modelled O–C orbital behaviour from each obtained pulse arrival time⁶ (see Fig 4). Note that the values supplied as supplementary material do not include this correction, as combining them with more data can change the required number of Fourier components and their coefficients.

Figure 4.

The grey dots in the top panel show the residuals as a function of orbital phase when fitting the pulse arrival times with only a second order polynomial. The red dashed line is a Fourier fit to the residuals, which show clear dependence on the orbital phase. The bottom panel shows the residuals when an orbital component modelled by a Fourier series is included in the fit. The black points show the root mean square binned to 50 points, fitted by a Fourier series (red dashed line). There are still phases with a large amount of scattering marked by the regions shaded in red, which were excluded from subsequent fits.

Open in new tabDownload slide

This same procedure cannot be applied to the CRTS data, since the pulses are not individually resolved at the CRTS cadence. We have instead determined the maxima by fitting a cosine to the data. We first phase-folded the data to the orbital ephemeris of Marsh et al. (2016), and modelled the orbital variability with a Fourier series. The orbital contribution was then subtracted from each measurement. Given the long-time span of the CRTS observations, the spin-down can have a measurable effect in the data. To account for this, we subtracted a quadratic term obtained from fitting the ULTRASPEC/ULTRACAM/HiPERCAM data from the measured times. We then fitted cosines with a period fixed to the beat period to the CRTS data divided into two data sets at the midpoint of observations (2008 September 17). Data taken at phases showing excess scatter (see next paragraph) were excluded. This resulted in two measurements for the maximum of the beat pulse extending our baseline back to 2005 with each representing the mean pulse time from three years of CRTS coverage.

Figure 5.

The y-axis shows the residuals (observed minus calculated, O–C) after the linear part of the ephemeris is subtracted, and the orbital modulation is corrected. The alternate grey/white columns mark different years, indicated at the top. The two left-most points correspond to the estimates from the CRTS data; the black bar at the bottom indicates the time coverage of these measurements (more details in the text). Measurements from 2015 onward correspond to ULTRASPEC/ULTRACAM/HiPERCAM data, and are shown in more detail in the inset. Individual measurements are shown in grey, and black points show the data binned every 20 points. The dashed line is our quadratic fit.

Open in new tabDownload slide

3.2 Fourier analysis

In addition to using our photometric data to perform pulse arrival time estimates, we also carried out a Fourier analysis. As in the previous section, we focused on observations taken with the g filter. To minimize systematic effects, we excluded any data points affected by clouds, and only utilized data employing our default comparison star.

To identify the main harmonics and combinations contributing to the observed waveform, we calculated the Fourier transform of the data, simultaneously modelled the identified dominant contributions with a Fourier series, subtracted the model from the light curve, and then re-calculated the Fourier transform of the residuals. The modelling was done by performing a least-squares fit to the data with orbital contributions modelled as sine waves, whereas the spin and beat contributions were modelled as cosines (we differentiate between sines and cosines to make it clear what phase 0 implies). Timing parameters were common to all data, whereas amplitudes and phases were fitted individually to each run. We started by subtracting the fundamental frequency and the first harmonic for the orbital (Ω), spin (ω), and beat (ω − Ω) frequencies. This revealed the next main contribution to be the second harmonic of the beat frequency 2(ω − Ω), followed by the combinations 2ω − Ω, 4ω − 2Ω, ω − 2Ω, the third harmonic of the beat frequency 4(ω − Ω), ω + Ω, and 2ω − 3Ω. At this point, the dominant peak became again the beat period, which is not perfectly modelled using this procedure due to the varying amplitude of the peak (see Stiller et al. 2018, and Fig. 6). Fig. 7 illustrates this pre-whitening process.

Figure 6.

Sample runs from different years. The orbital phase is accurate for the 2015 run, but we have slightly displaced the runs from other years to align the pulses in order to illustrate their changing amplitude.

Open in new tabDownload slide

Figure 7.

Fourier transforms around the orbital frequency (top panel), spin and beat fundamental frequency (middle-left), first (middle-right), second (bottom-left), and third (bottom-right) harmonics. The original Fourier transform is shown in grey with the peaks that were initially subtracted indicated in the same colour as vertical lines. Each subsequent iteration of our pre-whitening procedure is shown in the same manner, i.e. the Fourier transform and the peaks subtracted before the next iteration in the same colour. In black, we show the Fourier transform after our final model is subtracted. The black dashed line is placed at an amplitude of 0.01 (1 per cent) to facilitate comparison between the middle and bottom panels. Though some combinations remain visible in the amplitude spectrum, the main contribution is again the beat frequency. Given the sparse sampling over many years, many aliases are seen, in particular one-day aliases at this scale.

Open in new tabDownload slide

Once we identified the main contributions to the observed light curve, we carried out a fit of the light curve including all of these contributions using Markov-Chain Monte Carlo (MCMC) implemented with emcee (Foreman-Mackey et al. 2013). We initially performed a fit with amplitudes and phases of all contributions allowed to vary freely, but each set to be constant throughout the whole span of our data. This led to a poor fit of the light curve, as the amplitudes of the peaks are not consistent over time (as shown in Fig. 6). We then opted for initially fitting each of our runs individually with the timing parameters fixed to the values obtained from a least-squares fit to all of the data, but amplitudes and phases free. The obtained amplitudes and phases were then fixed, and we carried out the MCMC fit with only the timing parameters, |$T_{0}^{\mathrm{orb}}, T_{0}^{\mathrm{beat}}, \nu _\mathrm{beat}, \nu _\mathrm{spin}, \dot{\nu }_\mathrm{beat}, \dot{\nu }_\mathrm{spin}$|⁠, as free. Since our runs typically do not cover a full orbit, the orbital period itself is not well constrained from sinusoidal fits to our data, but it can be indirectly obtained by modelling the beat and spin frequencies separately. We therefore allowed the spin and beat periods to vary freely, and calculated the orbital period from their difference. The linear term corresponding to T₀ was allowed to be free for the orbital and beat periods, but constrained to be within one period of our initial guess. For the spin period, as well as spin/orbit combinations, T₀ was set to be equal to the beat T₀, but a phase shift was allowed. A phase shift was also allowed for all harmonics. In addition to the periods themselves, we also included the period derivative as a free parameter for both spin and beat in our model. We have not constrained them to be equal in order to probe for a significant orbital period derivative. Data points whose residuals were more than three standard deviations away from an initial least-squares fit were excluded from the MCMC fit to avoid occasional strong peaks, which typically have low uncertainties and could bias our fit. Uncertainties were inflated by a factor of eight to achieve reduced χ² ≈ 1 (though we report the χ² prior to this renormalization below). Additionally, we excluded data taken towards the end of our last run, 2022 April 26 (MJD > 59695.2125), because we noticed that the inclusion of these data led to an increase in χ² of at least a factor of two, leading to a poor solution. The reason behind this seems to be an inversion in the strength of consecutive peaks (see Fig. 8), possibly caused by a change in which pole generates the strongest synchrotron emission power. This is also eventually observed in other runs, but has a particularly damaging effect to the fit when occurring in the last portion of data, which plays an important role in constraining the frequency derivatives. With the current data, we see no evidence for a periodical behaviour of such inversion, which seems to be stochastic. Future observations will reveal whether this behaviour persists after our last run.

Figure 8.

The inversion of pulse strength observed in our last run. The data are shown as black dots with our Fourier model shown in red. Data beyond the grey dashed line was excluded from the Fourier fit. The typical behaviour of AR Sco consists of alternating strong (s) and weak (w) pulses, whose intensity can get very similar in particular for certain orbital phases. However, after MJD = 59695.2125 there was a particularly strong pulse where a weak pulse was expected. After some pulses with the normal behaviour, the inversion becomes persistent, which caused a large increase in χ².

Open in new tabDownload slide

When running the MCMC fit, we at first applied a prior forcing the difference between spin and beat derivatives to be consistent with the upper limit determined by Peterson et al. (2019), that is |$\dot{\nu }_\mathrm{spin} - \dot{\nu }_\mathrm{beat} \lesssim 3.8 \times 10^{-20}$| Hz s⁻¹. We obtained a solution with χ² = 18415210 (for 153527 data points), whose obtained parameters are given in Table 1. Next, we removed the prior and let |$\dot{\nu }_\mathrm{spin}$| and |$\dot{\nu }_\mathrm{beat}$| vary freely. Without the prior, we obtain a solution with a lower χ² of 18381704 for the same number of points. Fig. 9 shows a section of the data with this best-fit model, and Fig. B1 shows a corner plot of the resulting timing parameters. The obtained beat period and spin and beat frequency derivatives are consistent at the 99 per cent confidence level for the two solutions, though the value of the spin period and of |$T_{0}^{\mathrm{beat}}$| shift significantly. Notably, in both cases we obtain |$\dot{\nu }_\mathrm{beat}$| consistent with our pulse time arrival analysis. The main distinction between the two fits is that, without a prior, the difference between our spin and beat period derivatives implies a significant |$\dot{\nu }_\mathrm{orbit} = (8.43\pm 0.37) \times 10^{-19}$| Hz s⁻¹, higher than the upper limit suggested by Peterson et al. (2019), whereas the value is naturally lower, |$\dot{\nu }_\mathrm{orbit} = (3.68\pm 0.17) \times 10^{-20}$| Hz s⁻¹, when the upper-limit prior is applied. We further discuss this result in Section 4.1.

Figure 9.

The observed data are shown in the top panel as black dots connected by grey lines. The light curve model is shown in red. In the bottom panel, we illustrate the different contributions, which are from top to bottom: orbital (Ω and 2Ω, in green), beat (ω − Ω, 2ω − 2Ω, 3ω − 3Ω, and 4ω − 4Ω, in blue), spin (ω, 2ω, and 3ω, in purple), and other combinations (2ω − Ω, 4ω − 2Ω, ω − 2Ω, 2ω − 3Ω, and ω + Ω, in cyan).

Open in new tabDownload slide

3.3 Doppler tomography

Similarly to Garnavich et al. (2019), we have applied the Doppler tomography technique to our obtained spectra. This technique consists of mapping the observed line profiles at different orbital phases into velocity space, thereby allowing the line emission distribution in the binary to be mapped (see Marsh 2001, for a review). We started by analysing the Hα line using our WHT spectra. Garnavich et al. (2019), previously reported the existence of long-lived prominences (see their fig. 7), whose existence, they argued, requires the M-dwarf to have a strong magnetic field (100–500 G). Owing to the short integration time of our exposures, we were able to search for variations of these prominences with beat period, which could indicate that they are modulated by the interaction between the magnetic field of the white dwarf and the companion, since the beat phase tracks the orientation of the white dwarf with respect to the binary orbit. As can be seen in Fig. 10, these features seem to persist regardless of beat phase, suggesting that they are indeed inherent to the M-dwarf. In addition, they seem to appear only for the Hα line, whereas the emission for other lines is mostly concentrated on the face of the M-dwarf facing the white dwarf, though there are some substructures that vary from line to line (see Fig. 11). The detection of He ii is particularly important, as it suggests a high temperature on the M-dwarf’s heated face that would ionize material before it leaves the star, which invalidates Lyutikov et al. (2020)’s model (as previously pointed out by Garnavich et al. 2021). The location of this stream also suggests interaction close to the M-star, in agreement with the model suggest by Katz (2017) rather than to the white dwarf, as proposed by Lyutikov et al. (2020).

Figure 10.

Doppler tomography of the Hα line for different intervals of beat phase, indicated in the bottom left corner of each panel. The white cross marks the position of the white dwarf, the red line indicates the Roche lobe of the companion, and the red dot its centre of mass. The side lobes seen at all phases have been previously reported by Garnavich et al. (2019). An animated version of this figure is available online.

Open in new tabDownload slide

Figure 11.

Doppler tomography for the CaII H&K lines (left-hand panel), and for HeII 4686 Å (right-hand panel). Both emissions originate mainly from the irradiated Roche lobe of the M-dwarf, though the HeII line shows an excess near the trailing face of the star, as also seen by Garnavich et al. (2019).

Open in new tabDownload slide

It is important to note that persistent Hα prominences have been detected in other non-accreting systems, such as QS Vir (Parsons et al. 2016) or SDSS J1021+1744, in which they have been attributed to long-lived material near the L5 point (Irawati et al. 2016). Prominences are also shown by some CVs during low states, as reported for instance for AM Her (Kafka et al. 2008) and BB Dor (Schmidtobreick et al. 2012). Similar to what we see here, these CVs also show strong prominences mainly in Hα with emission in other lines appearing more concentrated near the heated face of the white dwarf companion. Kafka et al. (2008) and Schmidtobreick et al. (2012) noted that the prominences do not coincide with the location of the outer Lagrangian points, L4 and L5, but Schmidtobreick et al. (2012) hypothesized that the magnetic field could alter the Roche geometry and result in equilibrium points where the prominences are observed. A magnetic field has also been suggested by Parsons et al. (2016) to explain the prominences seen in QS Vir, and could also be responsible for the features seen in AR Sco.

One of the assumptions of traditional Doppler tomography is that the flux of any point in the binary system is constant over time. However, emission that is modulated with time is often observed in interacting systems (e.g. Papadaki, Boffin & Steeghs 2008; Calvelo et al. 2009; Somero et al. 2012). Steeghs (2003) extended the Doppler tomography method to allow for emission modulated on the orbital period. We applied this modulated Doppler tomography to the WHT spectra Hα line (Fig. 12) and to the Hβ (Fig. 13) line using the X-shooter spectra, which have a better spectral resolution than our WHT spectra.

Figure 12.

Modulation Doppler tomography maps of the Hα line following Steeghs (2003). The top left-hand panel shows the average amplitude (non-modulated), similar to Fig. 10, but including all beat phases. The top right-hand panel shows the total amplitude that is modulated on the orbital period. The bottom left-hand panel shows modulations on a sine (peaking at orbital phase 0.5), and the bottom right-hand panel shows modulations on a cosine (peaking at orbital phase 0.0). The modulated amplitude shows complex orbital-phase dependent structures, as further illustrated in Fig. 14. We include here what would be the trajectory of the accretion stream, if it existed.

Open in new tabDownload slide

Figure 13.

Same as Fig. 12, but for Hβ. Similar structures to those in Hα can be seen, but with much smaller amplitude.

Open in new tabDownload slide

Both lines show components with a significant orbital period modulation. Unlike the average amplitudes, which are somewhat dissimilar for Hα and Hβ given that only the former shows strong prominences, the modulated amplitude shows a similar behaviour for both lines. The sine component, which peaks at orbital phase 0.5 shows two main components, one of which has a similar origin as one of the observed Hα prominences, whereas the other is located near the trailing face of the M dwarf, much like seen for the HeII 4686 Å line. The cosine component (peaking at phase 0.0) is mainly along the line connecting the two stars, with an extension towards the same Hα prominence seen in the sine map. This phase-dependent behaviour is further illustrated in Fig. 14. The extension seems to be located between what would be the Keplerian and stream trajectories, if there were any mass transfer between the two stars. This is similar to what is seen in, for example, U Geminorum (Marsh et al. 1990). In AR Sco’s case, there is no evidence for accretion; yet, the existence of persistent prominences suggests the concentration of matter in the region where the stream would be. The modulated emission could thus result from a shock in this region resulting from interaction between the two stars’ magnetic fields as the system rotates.

Figure 14.

The left-hand panel shows the total modulated amplitude for the Hα line as a vector map. The direction of the arrow indicates the orbital phase with the arrow pointing upwards for phase zero and going clockwise from there (i.e. phase 0.25 to the right, and phase 0.5 downwards, and so on). The length of the arrows is proportional to the amplitude. The panel on the right is the same as the top right-hand panel in Fig. 12, but it is repeated here with the same scale as the vector map to facilitate comparison.

Open in new tabDownload slide

3.4 Rotational velocity of the M-dwarf

Evidence of ellipsoidal variation in AR Sco’s light curve suggests that the M-dwarf is filling its Roche lobe, or very nearly so (Marsh et al. 2016). For any star with a high-Roche lobe filling factor, the radius along the line of sight can change significantly as it rotates, which translates into variable projected rotational velocity, v_rot sin i, given that v_rot = 2πR_⋆/P_⋆ (where R_⋆ and P_⋆ are the radius and rotation period of the star). Our X-shooter spectra include several lines originating in the atmosphere of the M-dwarf that allow us to probe for this v_rot sin i variability.

Following Parsons et al. (2018), we selected three sets of lines for the estimates: the KI line at 7699 Å, the NaI lines at 8183/8193 Å, and the KI line at 1.252 μm. We have used the M4.5 template from Parsons et al. (2018), which we found to yield a better fit to these lines than the M5 template. Our approach was to find the value of v_rotsin i that minimized the χ² between the observed line and the template. We normalized and continuum-subtracted both observed spectra and template prior to fitting. We also corrected the observed spectra to the rest-frame of the M-dwarf.

In order to obtain model values for comparison, we utilized the lprofile code, which is part of the lcurve package (Copperwheat et al. 2010). We calculated phase-resolved synthetic line profiles assuming the system parameters derived by Marsh et al. (2016) and an orbital inclination of 60° (Potter & Buckley 2018b; du Plessis et al. 2019). The exposure length and the spectral resolution were fixed at values matching our X-shooter spectra. We estimated the full-width at half maximum (FWHM) of each line by fitting a Gaussian profile, and did the same for the 1.252 μm line (which is the most isolated of the analysed lines). We then interpolated the observed (close to linear) relationship between FWHM and v_rot sin i to estimate the v_rot sin i value of each synthetic profile. We carried out this estimate both for a completely Roche lobe filling M-dwarf, and for radii of 90, 80, 70, and 60 per cent of the Roche lobe radius (the latter corresponds approximately to the radius of an isolated M-dwarf of mass 0.3M_⊙). These radii are a linear measure of the distance from the centre of mass of the M-dwarf to L1, that is a linear measure of its radius relative to the Roche lobe, not a volume-based measurement. Our results are shown in Fig. 15.

Figure 15.

The top panel shows v_rot sin i measurements, whereas the radial velocity is shown in the middle panel. Green circles are measurements from the NaI 8183/8193 Å line, red triangles from KI 7699 Å, and blue squares from KI 12522 Å. The orbital phase was calculated using the ephemeris of Marsh et al. (2016). The dashed line in the radial velocity panel is a sinusoidal fit to the data, which was used to correct the wavelengths to rest frame before measuring v_rot sin i. The bottom panel shows residuals to this fit. The obtained parameters agree to within less than one-sigma with those reported by Marsh et al. (2016). In the v_rotsin i panel, the lines are for radii of 100, 90, 80, 70, and 60 per cent of the Roche radius (from top to bottom).

Open in new tabDownload slide

Though a lot of scatter is present in our measurements, in particular around orbital phase 0.5 when the analysed lines from the M-dwarf become much shallower, it can be clearly seen that a low degree of Roche-filling is incompatible with the estimated values. Whereas the reduced χ² for the 1.252 μm line is similar for models with full Roche lobe filling and for a radius of 90 per cent of the Roche lobe radius (4.5 and 6.0, respectively), the value more than doubles for a radius of 80 per cent (13.6), and increases by at least a factor of 5 for 70 and 60 per cent (34 and 62, respectively). Therefore, a radius of less than 80 per cent of the Roche lobe radius is very unlikely for the M-dwarf.

4 DISCUSSION

4.1 The rates of period change

Our pulse timing analysis builds on previous work (Marsh et al. 2016; Stiller et al. 2018; Gaibor et al. 2020) and further confirms the occurrence of spin-down, now confirmed at the 50-σ level through this method. Moreover, we have shown that the spin-down rate has been constant to within the observational uncertainties over a period of 17 yr, for which there are data. Applying a quadratic term correction based on other measurements to the CRTS times, as well as excluding data taken at phases showing scatter, were crucial for obtaining a good fit to the CRTS data. Marsh et al. (2016) did not take into account the large scatter at some orbital phases, which could explain the discrepancy with more recent measurements.

For the first time, we have also succeeded in obtaining an estimate of the beat frequency derivative via Fourier analysis, obtaining a value in good agreement with results from pulse arrival timing. In order to obtain a good fit to the data, amplitudes had to be fit individually in each run given that the pulse strength is clearly variable with time. In addition, data from our most recent observation for which the strength of consecutive pulses was not consistent with the typical behaviour had to be excluded. Possibly these measures have allowed us to obtain a Fourier fit consistent with the pulse timing analysis when previous attempts have failed.

One puzzling result of our Fourier analysis was a significant difference between our derived beat and spin frequency derivatives when they were left to vary freely, which implies a measurable orbital frequency derivative. Our value of |$\dot{\nu }_\mathrm{orbit} = (8.43\pm 0.37) \times 10^{-19}$| Hz s⁻¹ is in agreement with the upper limit set by Stiller et al. (2018), but is more than 20σ above the upper limit derived by Peterson et al. (2019).

It is worth noting that Peterson et al. (2019) remarked that they too obtained a non-null |$\dot{\nu }_\mathrm{orbit}$| if the timing uncertainties in their data were not accounted for (which would be unrealistic). As their data included photographic plate measurements, that is indeed a great source of uncertainty, and they point out that it can be up to 11 min for their data. Since the inclusion of uncertainties changes what is obtained for |$\dot{\nu }_\mathrm{orbit}$|⁠, that implies that the significance of |$\dot{\nu }_\mathrm{orbit}$| depends on the assumed uncertainties. Peterson et al. (2019)’s conclusion could, for example, be altered if their uncertainties were overestimated, that is if their timing measurements were more precise than they assumed. The timing uncertainties in our measurements were taken into account, though they are negligible.

Peterson et al. (2019) support their own empirical estimate by pointing out that the models of Knigge et al. (2011) predict |$\dot{\nu }_\mathrm{orbit} \lesssim 2 \times 10^{-20}$| Hz s⁻¹ for AR Sco’s orbital period. However, this prediction does not take into account that AR Sco is highly magnetic, which affects angular momentum loss (e.g. Kolb 1995; Cohen, Drake & Kashyap 2012; Belloni et al. 2020). The direction in which angular momentum loss is affected, however, is a subject of much discussion. According to Liebert & Stockman (1985) and King (1985), the white dwarf magnetic field would lead to enhanced magnetic breaking. They argued that the wind from the donor can be temporarily trapped by the white dwarf’s magnetic field lines, thus gaining angular momentum and subsequently carrying more angular momentum out of the system. Li, Wu & Wickramasinghe (1994), on the other hand, proposed the opposite: that the donor’s wind would remain trapped in the system due to the white dwarf’s magnetic field, leading to reduced angular momentum loss. This model was supported by modelling of the polar population done by Belloni et al. (2020).

Given that we find a higher |$\dot{\nu }_\mathrm{orbit}$| than the upper limit suggested by the non-magnetic semi-empirical model of Knigge et al. (2011), this would require AR Sco to be losing more angular momentum than a non-magnetic system, or enhanced magnetic braking. We estimate |$\dot{J}_\mathrm{magnetic} \gtrsim 40 \dot{J}_\mathrm{non-magnetic}$| given the ratio between our obtained value and the estimate using Knigge et al. (2011)’s model. Considering the uncertainties surrounding angular momentum loss in magnetic cataclysmic variables, particularly in the case of a unique system such as AR Sco, it is hard to establish whether this is feasible or not.

The orbital period change could instead happen at fixed angular momentum if it is triggered by structural changes in the secondary, which can be explained by magnetic activity. Activity can affect the distribution of angular momentum within the star, changing its rotational oblateness (i.e. the gravitational quadrupole moment), which translates into a change in the gravitational field. That in turn causes the orbit and speed of the stars to change, altering the orbital period of the system without affecting its angular momentum (the so-called Applegate mechanism, Applegate 1992; Lanza, Rodono & Rosner 1998; Lanza & Rodonò 1999). This model was proposed to explain the orbital period changes often observed for eclipsing systems (e.g. Parsons et al. 2010; Bours et al. 2016), which in many cases are well above the rate expected from magnetic breaking. In addition, the rate of period change caused by the Applegate mechanism is not constant, and can even change in signal (i.e. the period can go from decreasing to increasing and vice-versa). This could potentially explain the discrepancy between our derived value and the limit obtained by Peterson et al. (2019). This explanation of course requires the M-dwarf to be magnetically active, which seems to be the case at least for AR Sco.

4.2 AR Sco as a cosmic clock

The regularity and strength of the beat pulsations observed for AR Sco, as well as the precise ephemeris that can be obtained for the beat period, raise the question of whether AR Sco could be used for calibrating the timing accuracy of high-speed observations. This is potentially useful given that it is fairly bright, varies at a speed that it accessible to many instruments, and there is not a need to target particular eclipses. Whereas the ephemerides are precise to less than a second, one important issue is their colour dependency. As already pointed out by Gaibor et al. (2020), the times of the primary pulses depend on the choice of filter, and their difference varies with orbital phase. This is illustrated in Fig. 16 for a section of one of our HiPERCAM runs.

Figure 16.

A section of one of our HiPERCAM runs showing data obtained with the u, g, r, i, and z filters (from top to bottom). The vertical lines show the pulse arrival times obtained for each filter. In the two insets, it becomes clear that there is a varying level of difference in the pulse arrival times.

Open in new tabDownload slide

It is clear, therefore, that the use of AR Sco as a timing source depends on the desired accuracy. If the orbital phases around 0.05 and 0.5 are avoided, an accuracy of the order of one second could be obtained fairly easily (see top panel of Fig. 17), although a few pulses times may require combining if signal-to-noise is low. Around phases 0.35 and 0.85, the arrival times for different filters converge, so a better accuracy can be obtained, and one could perhaps hope to calibrate to a few tenths of a second. However, for large telescopes and high signal-to-noise, AR Sco is not fully reliable as a timing fiducial, and systematic errors from the source itself will limit precision. An overall important result is that, if the timing calibration is being done with a specific filter, preferably an ephemeris obtained with the same filter should be used. Nevertheless, AR Sco could prove useful in some instances, for example as a quick check against the large errors that can sometimes plague timing work, such as whether the start, middle or end of the exposure has been used, and whether barycentric corrections have been implemented correctly.

Figure 17.

Difference between pulse arrival times for u, r, i, z compared to g. The bottom panel shows the individual pulses for our two HiPERCAM runs, as well as one ULTRACAM run added to completely cover the orbital phases. The top panel shows the average every five pulses to illustrate the overall behaviour.

Open in new tabDownload slide

4.3 The behaviour of the companion

Our Doppler tomography confirms the existence of prominences around the M-dwarf, previously identified by Garnavich et al. (2019). We have further showed that those persist throughout the beat cycle, which suggests that they are not triggered by interaction between the M-dwarf and the white dwarf. Instead they are likely inherent to the M-dwarf, and remain stable due to its magnetic field. Further evidence for this is the fact that similar behaviour is observed for other non-accreting systems (e.g. QS Vir, Parsons et al. 2016), and even for CVs during low-states (e.g. BB Dor, Schmidtobreick et al. 2012), and is thus not unique to AR Sco. Likely the prominences are explained by the magnetic field altering the Roche geometry and displacing the location of the L4 and L5 Lagrange points.

By performing modulated Doppler tomography, we have identified the occurrence of emission modulated on the orbital period, in particular along the line connecting the two stars, but also remarkably in the region where an accretion stream would be located if there were mass transfer. A possible cause for this is a shock wave in the long-lived prominence that is located in this region, triggered by interaction between the two stars’ magnetic fields.

In addition, by measuring the v_rotsin i for different orbital phases, we provide evidence that the companion is likely filling at least 80 per cent of its Roche lobe.

5 SUMMARY & CONCLUSIONS

We performed an extensive study of the white dwarf pulsar AR Sco using precise photometric data spanning almost seven years, complemented by CRTS data that extends our baseline to 2005. We refined the beat period ephemeris using pulse time arrival analysis, detecting a non-zero frequency derivative at the 50σ level. We also confirm this value via Fourier analysis. We performed two different fits, first constraining the difference between |$\dot{\nu }_\mathrm{spin}$| and |$\dot{\nu }_\mathrm{beat}$| to be no more than 3.8 × 10⁻²⁰ Hz s⁻¹, the |$\dot{\nu }_\mathrm{orbit}$| upper limit found by Peterson et al. (2019), and next leaving both to vary unconstrained. Both cases lead to significant values of |$\dot{\nu }_\mathrm{orbit}$|⁠: in the first case (3.68 ± 0.17) × 10⁻²⁰ Hz s⁻¹, and (8.43 ± 0.37) × 10⁻¹⁹ Hz s⁻¹when no prior is applied. This could suggest that AR Sco is undergoing orbital period changes in excess of what is predicted by magnetic breaking and gravitational wave losses, similar to what is observed for many eclipsing systems and often attributed to changes in the quadrupole moment of the M-dwarf due to magnetic activity. However, since the χ² values obtained for our two fits are comparable and they lead to very different |$\dot{\nu }_\mathrm{orbit}$| estimates, we believe that there is no strong evidence for measurable orbital period change at the moment, though continuous monitoring might change this picture. There is also no reason to believe |$\dot{\nu }_\mathrm{orbit}$| to be constant, as its value can depend on the magnetic activity of the M-dwarf. Large changes have been observed for post-common envelope systems (Parsons et al. 2010; Bours et al. 2016).

We additionally carry out an analysis of time-resolved spectra obtained with both ISIS/WHT and X-shooter/VLT. Modulated Doppler tomography performed here for the first time shows the occurrence of emission modulated on the orbital period in the region where an accretion stream would be. This region also shows long-lived prominences that appear inherent to the M-dwarf, given their lack of dependence on beat phase. The modulated emission could thus be the result of shock interaction between the two stars’ magnetic fields in this region as the system rotates. Our spectra also reveal evidence of tidal deformation detected via the variation on the observed projected rotational velocity, which indicates that the M-dwarf must be at least 80 per cent Roche lobe filling.

AR Sco remains a unique and puzzling system. If our findings of a possibly significant |$\dot{\nu }_\mathrm{orbit}$| are confirmed, they could require the system to be losing angular momentum at a higher rate than non-magnetic systems, which could provide clues to explain the origin of AR Sco. Alternatively, this could provide a probe for the mechanism proposed by Applegate (1992), which allows for orbital period changes at fixed angular momentum. The detection of modulated emission can shed light on the interaction between the M-dwarf and the white dwarf, which remains not fully understood.

SUPPORTING INFORMATION

suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We thank D. Steeghs for advice on how to best visualize modulated Doppler maps. IP and TRM acknowledge support from the UK’s Science and Technology Facilities Council (STFC), grant ST/T000406/1. AA acknowledges funding support from the NSRF via the Program Management Unit for Human Resources & Institutional Development Research and Innovation, grant B05F640046. SGP acknowledges the support of a STFC Ernest Rutherford Fellowship. This work has made use of data obtained at the Thai National Observatory on Doi Inthanon, operated by NARIT, and of observations made with the Gran Telescopio Canarias (GTC), installed at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias in the island of La Palma. The design and construction of HiPERCAM was funded by the European Research Council under the EU Seventh Framework Programme (FP/2007-2013) under ERC-2013-ADG Grant Agreement no. 340040 (HiPERCAM). VSD and ULTRACAM/ULTRASPEC/HiPERCAM operations are supported by STFC grant ST/V000853/1. This works is also partially based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 097.D-0685(B), and on observations made with WHT operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. The ISIS spectroscopy was obtained as part of W/2016A/26.

DATA AVAILABILITY

All data analysed in this work can be made available upon reasonable request to the authors.

Footnotes

REFERENCES

APPENDIX A: LOG OF PHOTOMOMETRIC OBSERVATIONS

APPENDIX B: CORNER PLOT FOR THE MCMC FIT

Figure B1.

Corner plot resulting from the unconstrained MCMC fit to the data. The median values quoted in the text have been subtracted for both T₀ and period values. We show |$\dot{\nu }_\mathrm{beat}$| and |$\dot{\nu }_\mathrm{beat}$| with the same scale in both axes to facilitate comparison.

Open in new tabDownload slide