Optical, X-ray, and $\gamma$-ray observations of the candidate transitional millisecond pulsar 4FGL J0427.8-6704

We present an optical, X-ray, and $\gamma$-ray study of 1SXPS J042749.2-670434, an eclipsing X-ray binary which has an associated $\gamma$-ray counterpart, 4FGL J0427.8-6704. This association has led to the source being classified as a transitional millisecond pulsar (tMSP) in an accreting state. We analyse 10.5 years of Fermi LAT data, and detect a $\gamma$-ray eclipse at the same phase as optical and X-ray eclipses at the>5$\sigma$ level, a significant improvement on the 2.8$\sigma$level of the previous detection. The confirmation of this eclipse solidifies the association between the X-ray source and the $\gamma$-ray source, strengthening the tMSP classification. However, analysis of several optical data sets and an X-ray observation do not reveal a change in the source's median brightness over long timescales or a bi-modality on short timescales. Instead, the light curve is dominated by flickering which has a correlation time of 2.6 min alongside a potential quasi-periodic oscillation at $\sim$21 min. The mass of the primary and secondary star are constrained to be $M_1=1.43^{+0.33}_{-0.19}$ M$_{\odot}$ and $M_2=0.3^{+0.17}_{-0.12}$ M$_{\odot}$ through modelling of the optical light curve. While this is still consistent with a white dwarf primary, we favour the transitional millisecond pulsar in a low accretion state classification due to the significance of the $\gamma$-ray eclipse detection.


INTRODUCTION
"Redbacks" are binary star systems which have a neutron star (NS) primary and a low-mass, near-main sequence companion. The neutron stars in these systems are detectable at radio and γ-ray wavelengths as millisecond pulsars (MSPs). In recent years, 3 "redback" systems have become increasingly important in understanding the evolution of MSPs in binary systems: PSR J1023+0038 (Archibald et al. 2009;Stappers et al. 2013), IGR J18245-2452 (Papitto et al. 2013) and PSR J1227-4853 (Bassa et al. 2014). These three systems have been observed to transition between a radio loud state, where the pulsar is detectable at radio wavelengths and there is no evidence for active accretion from the secondary, and a low-mass X-ray binary (LMXB) state, where emission from the radio pulsar is quenched and material Email: kennedy.mark@manchester.ac.uk flows from the secondary through the inner Lagrange point towards the NS primary, where it builds an accretion disc.
This paper focuses on one of the candidate tMSP systems, 1SXPS J042749.2-670434, which is a binary system with an 8.8 hour orbital period. The X-ray source has been associated with the bright γ-ray source 3FGL J0427.9-6704. Since the publication of the original γ-ray association by S16, the Fermi LAT 8-year Source Catalog has been released (4FGL; The Fermi-LAT collaboration 2019), and the γ-ray source has been renamed 4FGL J0427.8-6704. This as-sociation has been made due to the presence of deep X-ray and optical eclipses in the light curve of 1SXPS J042749.2-670434 source which potentially coincide with an eclipse of the γ-ray source (S16), with the γ-ray eclipse only detected at the 2.8σ level. Hereafter 1SXPS J042749.2-670434 and 4FGL J0427.8-6704 are assumed to be the same source, and are collectively referred to as J0427. Due to the lack of an accurate mass measurement, questions remain over whether the primary star in J0427 is a white dwarf (making the system a cataclysmic variable) or a neutron star/black hole (making the system a low-mass X-ray binary).
Based on the presence of γ-ray emission from the binary, and by estimating the primary mass using optical spectroscopy and photometry of the secondary star, S16 have suggested that the primary star is likely an MSP, and that the system is a tMSP in the accreting state. However, at the time of publication, there have been no dedicated radio observations of J0427 reported in the literature, while a positive radio detection of the source would strengthen the classification of the system as a tMSP. Mudding the waters further is the conclusion by S16 that the mass constraints on the primary are not strong, and that a white dwarf primary could still be possible, potentially making this system a member of the cataclysmic variable (CV) family of interacting binaries, which have white dwarf primaries. If this were true, the main issue would then be explaining the γray emission from the source, as no known CV has a GeV γ-ray counterpart.
Here we present high time-resolution optical photometry of J0427 taken in u s , g s , and i s filters simultaneously, X-ray data taken using XMM-Newton, optical data taken using TESS over 11 months, and 10.5 years of Fermi data in an attempt to strengthen the detection of a γ-ray eclipse.
The optical photometry reveals rapid flickering associated with an accretion disc, and tantalising hints of a periodicity close to ∼20 min, which we investigate using a combination of traditional timing analysis and Gaussian process modelling. The underlying orbital modulation coupled with the recent distance estimate to the system measured by Gaia suggests that the primary is a neutron star, but the classification of J0427 as a tMSP in an accreting state still hinges on the observed γ-ray emission from the source, with no other evidence to support the classification as a tMSP over a regular low-mass X-ray binary.
The data were reduced using the ULTRACAM pipeline as described in (Dhillon et al. 2007). The magnitudes were calibrated using previously determined zero points in each filter. Our band calibration is good to ∼ 0.1 mag in the i s and g s bands based on the measured magnitudes of two nearby comparison stars. Unfortunately there are no measured SDSS u s magnitudes for any star in the field, but based on the accuracy of the g s and i s calibration we estimate the maximum band calibration uncertainty in u s to be ∼ 0.1 mag.
Throughout the remainder of the text, any mention of u s , g s , or i s refers to data taken with ULTRACAM.

TESS
J0427 has been observed by CCD 4 of the Transiting Exoplanet Survey Satellite (TESS ) during each Sector of Cycle 1 up until the submission of this paper (S001-S011). TESS records full frame images every 30 mins during a Sector, with observations of each Sector lasting for 28 days. The data are taken through a wide filter which covers 6000-10000Å (λ c = 7865Å). The TESS cut out images around J0427 were downloaded for all available Sectors. Extraction of a calibrated light curve was not possible, as the 21" pixel size of TESS means that J0427 is blended with several sources in the images. We constructed custom source and background apertures for each Sector, extracted the source aperture flux and removed background variations and flagged data where the background exceeded 100 electrons s −1 pixel −1 , and then subtracted the mean flux value from the source light curve such that the residual light curve only showed variation in the source aperture. Even after this procedure, there were still many outlier points in the light curve. To remove these, we σ-clipped the data, setting σ = 3.5, removing a further 606 of the 12545 data points. The resulting light curve over all 11 sectors is shown in Figure 1, and the amplitude of the variation was ∼ 1.5 e s −1 . There are still clear systematics in this light curve, particularly at the start and end points of each Sector. The extracted light curve for J0427 from TESS Sectors 1-11. The aperture flux has been background subtracted, median subtracted, and 3.5 σ clipped. The amplitude of the variation is ∼ 1.5 e s −1 . Time is given in Barycentred TESS Julian Date (BTJD), which is BJD-2457000.0. Gaps in the data are due to gaps in data acquisition between sectors.

XMM-Newton
Right: A Lomb-Scargle periodogram of the TESS data. The dashed lines mark the known orbital frequency and its first 8 harmonics.
of 2.4 orbital periods. The European Photon Imaging Camera (EPIC) -pn (Strüder et al. 2001), -MOS1, and -MOS2 (Turner et al. 2001) CCDs were all operated in Full Frame mode with a medium filter inserted. The Optical Monitor (OM; Mason et al. 2001) was operated in fast mode, with a white filter inserted. While both Reflection Grating Spectrographs (den Herder et al. 2001) were operational, these data will not be discussed as no appreciable signal was detected.
The data were reduced using the Science Analysis Software (SAS) v16.1.0. The PN and MOS data were processed using the SAS commands Epproc and Emproc respectively. Unfortunately, the 77.5 ks exposure suffered from periods of severe high background. Figure 2 shows the 0.3-10 keV light curve alongside the high energy background light curve which is used to characterise the intensity of soft proton flaring during observations with the EPIC instruments. Any period when this flaring activity is above 0.5 counts s −1 (∼ 75% of the observation duration) had to be discarded for spectral analysis. All of the X-ray data were considered when looking for X-ray eclipses in the X-ray light curve.

Fermi LAT
For the analyses reported below, we selected Pass 8 (Atwood et al. 2013) data from the Fermi Large Area Telescope (Atwood et al. 2009), collected between 2008-08-04 and 2019-01-08 with reconstructed energies 0.1< E <30 GeV, reconstructed positions lying within 2 • of J0427, and with an event class of 128 1 . The associated source is modeled with a log parabolic shape, and we use the source list and associated diffuse models to compute the probability ("weight", Kerr 2011) that each photon is associated with the FL8Y counterpart or with a background source.
1 In line with the chosen instrument response function of P8R3 SOURCE V2, see https://fermi.gsfc.nasa.gov/ssc/ data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/ IRF_overview.html Background Flux (cnts s 1 ) Figure 2. The extracted X-ray light curve of J0427 (black and grey) alongside the high energy background light curve (red). The grey X-ray data were excluded from spectral analysis due to the high background during these times. Figure 3 shows the 3 bands of ULTRACAM data, the TESS data, the X-ray light curve from the EPIC-pn instrument in the 0.3-10 keV band, the OM light curve, and the Fermi LAT light curve (which is based on an analysis of Fermi data which wll be discussed in Section 3.1), all phased to the orbital period using the ephemeris given in S16 of

PHASED LIGHT CURVE
where T mid is the predicted time of mid eclipse in Barycentric Julian Date and E is the orbit number, with cycle 0 occurring at BJD 2455912.83987. The light curve shows a deep eclipse at phase 0 in each of the data sets, and lasts for 0.08 orbits (∼42 min). Outside of eclipse, the optical and X-ray light curves are dominated by rapid flickering which occurs on a timescale of minutes. This variability is harder to detect in the TESS data, as these data had a cadence of 30 min. Aside from the eclipse and flickering, the ULTRA-CAM and TESS data show a hint of curvature on the orbital timescale which is strongest in the TESS and i s data and is undetectable in the u s data. The TESS data show evidence of a secondary eclipse at phase 0.5, which becomes more obvious after the data have been binned with a bin width of 0.01 in orbital phase. The X-ray eclipse has a similar duration to that of the optical eclipse, and the X-ray light curve is consistent with 0 flux being observed during the eclipse, suggesting the entire X-ray emitting region is being eclipsed.
3.1 γ-ray Eclipse S16 reported tentative evidence for a γ-ray eclipse at the same orbital phase as the eclipses observed in the optical and X-ray bands. To verify this claim, we folded the Fermi photon timestamps with the J0427 orbital period and set the zero phase to fall on BJD 2457527.69139; this epoch is advanced by 0.25 in phase relative to that of S16, so that the eclipses fall at φ = 0. We modelled the eclipse (or excess) using a top hat model such that within an eclipse of width θ, the relative source rate is α. To enforce an average intensity of unity, the rate outside of the eclipse is thus (1− αθ)/(1− θ). With this eclipse model, the Poisson likelihood is where Θ indicates the phases of eclipse andΘ the complement, S is the total expected source counts (S ≈ i w i ), and η Θ represents the fraction of the instrument exposure falling in the eclipse. We can approximate the eclipse shape as a truncated Fourier series, and Kerr (2019) discusses the efficient evaluation of such likelihoods using Fast Fourier Transforms, allowing us to evaluate the likelihood over a wide range of trial orbital frequencies. We perform this search using a 40-term Fourier series to approximate the eclipse profile and search over eclipse width (0.01 < θ < 0.5), position (0 ≤ θ 0 < 1), and amplitude (α > 0), taking the maximum value of the likelihood for each trial frequency. Compared to a uniform signal, we find δ log L = 13.4 for an eclipse of amplitude α = 0.01 and width θ = 0.064 centred at θ 0 = 1.000. Using this shape, but allowing α to vary, we searched 10 5 neighboring orbital frequencies and found no signals with δ log L > 8.8, indicating a lower limit on the chance probability of ∼10 −5 . Moreover, because there is a single degree of freedom (α), δ log L should follow a χ 2 1 distribution, in which case the chance probability of observing δ log L = 13.4 is 2.2 × 10 −7 , indicating the eclipse has > 5σ significance. Moreover, the shape and phase are consistent with eclipses observed at lower frequency.
To further characterize the eclipse, we first performed a simple maximum likelihood analysis to determine the relative source flux at 40 orbital phase bins, following the methods of Kerr (2019, submitted), in which the photon weights are used to approximate a full time-domain likelihood analysis. The measurements (with 1σ errors and 95% confidence upper limits) are shown in blue in Figure 4, clearly show the eclipse. To identify any sharp features, we also apply a Bayesian blocks algorithm (Scargle et al. 2013) using 1000 orbital phase bins and the Poisson likelihood. For a wide range of priors on the number of change points, the algorithm finds only two significant intervals, shown as the red levels in Figure 4, indicating a strong preference for a deep eclipse with sharp edges. To confirm this, we directly fit the unbinned likelihood to find the best-fit values for ingress, egress, and flux within the eclipse, obtaining 0.963 − 1.041 and α = 0.13. The edges are relatively sharply constrained, with steep decreases in L indicating the left edge lies at φ ≥ 0.962 and the right edge at φ ≤ 1.047, while the flux during eclipse is only poorly determined, 0 ≤ α ≤ 0.3. The eclipse is apparently highly symmetric about the compact object superior conjunction, and the egress may be more gradual than the ingress.
To consider a more gradual eclipse, we modelled the eclipse as a Gaussian centred on , and found maximum likelihood parameter values of θ 0 =1.001(6), σ = 0.023(5), and A = −1.09(15). The resulting model is shown as the green trace in Figure 4. However, the improvement in δ log L = 12.42, compared to δ log L = 15.10 for the top hat model, suggests that the more rapid eclipse is the preferred model. Note that the likelihood calculation here uses an exact top hat representation, rather than the truncated Fourier series used in the frequency search reported above, yielding a slightly higher δ log L.

Orbital Period
In order to try to improve the orbital period given by S16, we scaled the i s ULTRACAM data to match the amplitude of the variability in the TESS data, and then computed a Lomb-Scargle periodogram (Lomb 1976;Scargle 1982) of the combined data set. The strongest peak was located at a period of 0.3667 d. We estimated an error on this period by minimising the χ 2 value of a fit to the data using a Super-Smoother algorithm (Friedman 1984), as done in S16. The resulting best fit period was 0.36672±0.00002 d, which is not an improvement on the previous reported period. We also attempted to use Gaussian Process Modelling (which will be further discussed in Section 4.5) with a periodic kernel to determine an accurate period from the TESS data. Such methods have proved effective in the past when dealing with data that contain flickering (Littlefair et al. 2017;Angus et al. 2018). However, when performed solely on the TESS data, this method constrained the orbital period to be 0.366719(8) days, which is consistent with the previous reported error.

Orbital modulation
To better constrain the binary system parameters in J0427, we modelled the i s , g s , and u s ULTRACAM data using the Eclipsing Light Curve (ELC) code (Orosz & Hauschildt 2000). Due to the rapid flickering in the light curve, we median binned each night of data with a bin width of 10 mins, with the error on each binned point taken to be the standard deviation of the points which make up that bin. This binned light curve is shown as the black points in the top three panels of Figure 3. The TESS light curve was excluded from the modelling as the magnitudes could not be calibrated due to  Figure 3. The phased light curve of J0427 in i s , g s , u s bands taken with ULTRACAM (top 3 plots), the TESS data (orange, middle), the XMM-Newton OM data (grey, 3rd from bottom), 0.3-10 keV X-ray data (2nd from bottom, brown), and Fermi LAT data (pink, bottom, upper limits marked). Overlap from the 3 nights of ULTRACAM observations obscures some of the flaring activity around orbital phase 0.6. Orbital phase has been repeated for clarity. The black data points in the ULTRACAM data show the binned light curve which was used for modelling in Section 3.3, while the black points in the XMM-Newton data are the raw data binned with a bin width of 0.025 in orbital phase. The shaded region highlights the eclipse data which were excluded when performing the timing analysis in Section 4.4. Relative Intensity Figure 4. The γ-ray light curve obtained using Fermi-LAT data and following the methods described in Kerr (2019, submitted). The 40 blue points give maximum likelihood estimators for the intensity and its uncertainty, or upper limits, while the red markers give the same estimates for the two intervals determined with the Bayesian Blocks algorithm. For these two intervals, the extent on the x-axis indicates the interval boundaries, while the y-axis values give the intensity and uncertainty. The green line represents the maximum likelihood result from fitting a Gaussian eclipse model directly to the photon phases and weights. Finally, for comparison, the orange points show NuSTAR data, folded on the J0427 orbital period as described in Strader et al. (2016), scaled to give similar intensity to the γ-ray data.
several stars lying within the aperture used for source extraction, while the OM light curve from XMM-Newton was excluded as the data were taken using the white filter, which has a very large bandwidth.
There are many tunable parameters within ELC which are used to described the primary star, accretion disc, and secondary star. The variable parameters were the mass ratio q = M 2 /M 1 (where M 1 is the primary mass and M 2 is the secondary mass), the system inclination i, the outer radius of the accretion disc R out , the temperature of the disc at the inner radius T d , the opening angle of the accretion disc β, the effective temperature of the secondary star T 2 , the irradiation luminosity of the primary which is responsible for heating the secondary star L X , and the apparent radial velocity of the secondary star K 2 . Note that T 2 is the average temperature over the entire surface of the secondary including the irradiated region of the secondary, meaning if irradiation of the secondary is high, then T 2 is likely to be much higher than the night side temperature of the star.
In the following modelling, we assumed that the powerlaw exponent which controls the temperature profile of the disc (that is, where r in is the inner radius of the accretion disc) was ξ = −0.425, larger than what is expected from the typical ξ = −0.75 in "steady-state" accretion discs, but in-line with an irradiated accretion disc (Hayakawa 1981) which is expected in this system due to the strength of the detected X-ray source. To ensure this choice of ξ was not biasing our results, we also ran a Markov chain Monte Carlo (MCMC) as implemented as part of the ELC package (see Tegmark et al. 2004) with ξ free. The only change this made to our final results were slightly higher errors on T d , T 2 , and R out . During this modelling, ξ quickly converged to −0.4 ± 0.1, justifying the above assumption. Table 1 lists all of the relevant physical parameters used by ELC to model J0427, and whether or not the parameter was fixed or allowed to vary.
All 3 optical bands of data were fit simultaneously. The best fit model was found by using a MCMC. We allowed 70 chains to evolve over 8000 steps each, discarded the first 250 steps of each chain as burn-in, and also filtered out models which had a χ 2 which was more than 100 larger than our best fit model. The priors on all parameters were flat with hard edges at minima and maxima values, as given in Table 1.
We also assumed a Gaussian prior on the radial velocity of the secondary star of K 2 = 293 ± 4 km s −1 , in line with the measured value in S16. This K 2 value may be lower than the actual K 2 of the secondary in the system due to heating effects, but is still useful in providing a lower limit on the mass of the companion star.
The final number of samples used for the following analysis was 680,642. The corner plot of the MCMC analysis is shown in Figure 5, and our best fit model is shown alongside each optical band of data in Figure 3.3. The corner plot shows that there are several degeneracies between various parameters (q and i, T 2 and L X , T 2 and T d ) and that several of the parameters are not well constrained using our data. Our final fit had a χ 2 of 165.18 for 185 degrees of freedom, suggesting the error bars on the binned data points were over-estimated. The best constrained parameters are i = 84 ± 3°, T 2 = 5300 ± 600 K, and L X = 1 +0.9 −0.5 × 10 35 erg Table 1. Best fit values from the MCMC analysis of ELC. In the case where a parameter is given, the prior was assumed to be flat with the exception of K 2 , which had a Gaussian prior with a σ equal to the measurement error from S16. If a prior is not given, that means the parameter was fixed to this value.
a R out is expressed as a fraction of the primary star's effective Roche lobe radius, and so must be < 1. b units of log 10 (erg s −1 ) s −1 , where the errors have been scaled such that the χ 2 R of the best fit model was 1. (i) A very rapid reddening of the object as the colour moves from the centre to the top right of the plot which occurs during the eclipse.

Colour information
(ii) A rapid, large variation in the colour which sees J0427 moving from the centre of the plot to the bottom left (the sources of these variations will be discussed in detail in a later section), (iii) A slow, small variation in colour due to the orbital motion of the secondary which sees the colour moving from the centre of the plot towards the upper left between orbital phases 0 and 0.5, and back down towards the bottom right from orbital phase 0.5 to 1.0.
The gradual reddening of the object from orbital phase 0.0-0.5 and the reverse from phase 0.5-1.0 confirms the proposition made at the beginning of this section that the i s light curve has significant curvature over the orbital period while the u s light curve is relatively flat, and is line with the modelling discussed in the previous subsection.

Flux distribution
The ELC model found in Section 3.3 was then subtracted from the i s , g s , and u s ULTRACAM light curves, such that the remaining signal only included noise and short time-scale flickering. We then generated histograms of the residual flux in each band to look for the bimodal distributions which have been seen in several tMSPs during their active states (Shahbaz et al. 2015;Britt et al. 2017). The flux distributions only had a single peak with a high flux tail in each band. The same is true of both the TESS and XMM-Newton OM light curves.
There is a segment of data taken on 2017-10-16 during which the variability of the flux from J0427 decreased for ∼ 15 min. The TESS and XMM-Newton OM data prove that this feature is not an orbital feature, as neither light curves show a drop in the average flux at the same orbital phase that this period of diminished variability was seen in the ULTRACAM light curve. This means that this drop was transient, and does not repeat every orbital period. Figure 8 shows this segment of data alongside a 29 min segment of the light curve of the tMSP PSR J1023+0038 which was also obtained with ULTRACAM when PSR J1023+0038 was in an accreting state on the night of 2019-03-01 using 10 s exposures. Both light curves show the same amplitude of variability. However, there are still differences in the light curves. In particular, the mode switching behaviour of PSR J1023+0038 has a very distinctive step shape, with clear plateaus during the high mode, while the features in J0427 more closely resemble flares rather than the bimodel behaviour seen in PSR J1023+0038.

Photopolarimetry
The photopolarimetric observations were binned to 300 s exposure times for analysis. We did not find any indication of linear polarization above the sky background. However, due to the faintness of J0427 and the moonlit linearly polarized sky, we were only able to place an upper limit of ∼ 3% on the linear polarization of J0427. Circular polarization displayed random excursions around zero percent, consistent with a non-detection with a limit of < 1%.

SHORT TIME-SCALE VARIABILITY
While the orbital modulation dominates the long-term optical variability of J0427 there are also strong short-term flares visible in the data. We begin our discussion on this short-term variability by first considering the TESS data (due to its long cadence), followed by the XMM-Newton OM data, and concluding with the behaviour seen in the highest quality data, the ULTRACAM light curve. The 3σ levels as plotted in the various periodograms in this section were derived as described in Appendix A.

TESS
The periodogram of the TESS data is shown in the right panel of Figure 1 up to a frequency of 0.27 mHz (equivalent to a minimum period of 60 min, which is twice the sampling rate for the TESS data). The periodogram shows strong power at the orbital frequency and its first 8 harmonics. There are a further two potential signals located close to 3 × 10 −4 Hz which are not clearly associated with the orbital period. In order to investigate whether these peaks were real, noise, or related to a beat between the orbital period and the cadence of the data, a function describing the orbital modulation of the light curve was created by fitting the binned TESS light curve shown in Figure 3 using the numpy interpolate feature. This model light curve was then subtracted from the data, and a power spectrum of the residuals taken.  Table 1. There are clear correlations between q and i, q and R out , T d and T 2 , T d and L x , and T 2 and L x . We could not constrain K 2 better than what was done in S16. Dashed vertical lines mark the median value and 3σ confidence intervals.
This power spectrum showed no strong features, suggesting these two peaks seen on the far right of Figure 1 are not intrinsic to the system.

XMM-Newton OM
We next took the optical light curve obtained using the OM on-board XMM-Newton and constructed a power spectrum up to a frequency of 11 mHz (a minimum period of 1.5 min). The power spectrum is shown in Figure 9. There are three frequencies where the detected power surpasses the 3σ level -ν O (the orbital frequency), 3 × ν O , and 0.63 mHz. The detection of the first two frequencies is unsurprising given that the light curve shows eclipses, but the detection of power at 0.63 mHz (equivalent to a period of 26.5 min) is surprising. One possibility is that this is related to the window function of the observations -while the length of a single observation taken with the OM lasts for 20 minutes, there is also 6 minutes of dead time between sequential observations, meaning that the OM data has a sampling rate of 26 minutes. This could be giving rise to the signal at this period. Considering, however, that the 26.5  min period does not show up in our simulated light curves which have the same window function as the actual data, the period is likely not related to this, and is instead intrinsic to the system.

ULTRACAM
After removing the orbital modulation and masking the eclipse from the original light curve, the remaining data showed dramatic variability typically associated with flickering in an accretion disc. We subjected the data to a variety of time-series analysis techniques to look for short term periodicities. We treated each night of data individually to avoid the heavy aliasing which would arise had we combined all of the observations together.

Autocorrelation Function
The autocorrelation function (ACF) for each night of data was calculated after the orbital modulation had been removed. For the data from 2017-10-16 we excluded all data taken after the eclipse so that this feature would not alter    Figure 9. Power spectrum of the data obtained using the OM onboard XMM-Newton (grey). The power spectrum was binned (black points) and fit with a power law (green, dashed line) and broken power law (blue, dashed-dotted line) model. The power spectrum was best fit by the broken power law model, with the break occurring at 1.6 × 10 −3 Hz. The red line represents the 3σ level for a detected signal based on simulated light curves. There are three peaks located above this confidence level -one at the orbital frequency (marked as ν o ), one at three times the orbital frequency (3ν o ), and one at 0.63 mHz (ν QPO ), which corresponds to a period of 26.5 min.
the ACF, and no gaps in the data would be present. The ACFs from each night showed correlations with a time scale of ∼ 2-3 min (see Figure A1 in Appendix A). The ACFs from 2017-10-14 and 2017-11-22 showed no other strong correlations. The same can not be said of the data from 2017-10-16. This ACF has not only a short time-scale correlation of about ∼ 2-3 min, but also a near sinusoidal correlation with a period of 21 minutes which is only coherent for the first 150 min of the observation.

Power Spectrum
Lomb-Scargle periodograms of each night of data showed strong power at varying periods between 10-30 minutes, as shown in Figure 10. However, since the system contains an accretion disc (based on the double peaked optical emission lines reported in S16), we expect the power spectra to be dominated by pink noise. This is correlated noise with a power-law spectrum f α , with α = −1, and is not unusual in systems which are dominated by "flickering" (Lawrence et al. 1987). The presence of pink noise makes the determination of the significance of any peaks in a power spectrum difficult, as the usual metrics such as the False Alarm Probability depend on the noise in the data being uncorrelated (α = 0.0). As such, the 3σ level for whether a period was real or a noise peak for each band and each night was found using the method discussed in Appendix A. The only power spectrum which had a peak detected above the 3σ level was from the data taken on 2017-10-16 in the i s , and the peak was found to lie at ∼21 min. While there is a peak in the g s and u s light curve around the same period, it is not detected at the 3σ level. This may be due to a combination of the lower temporal sampling and lower S/N of the u s and the lower S/N of the g s data when compared to the i s data.

Gaussian Process Modelling
We next attempted to investigate the nature of the flickering using Gaussian Process Modelling (GPM) methods. A Gaussian Process is a collection of random variables, any finite number of which have a joint Gaussian distribution defined by a certain covariance matrix. While traditionally the covariance matrix is calculated from the data, in GPM the matrix is constructed using covariance kernels which attempt to model the correlations between data points. We tested two different kernels when analysing the optical data presented here. The first kernel was the sum of a a Matérn covariance kernel to allow for covariances between data points on a length scale of 1 and a periodic kernel which allowed for variations on a time scale equal to the orbital period. The second kernel had a third component to allow for short-term periodic variations in addition to these two components of the first kernel. A detailed discussion of the kernels used in this analysis is included in Appendix A.
The kernels were implemented using scikit-learn (Pedregosa et al. 2011) and the hyper-parameters were tuned using a Markov Chain Monte Carlo (MCMC) sampler implemented using emcee (Foreman-Mackey et al. 2013). As with our light curve modelling, P or b was fixed to 0.3667200 d. The kernels were initially trained on the combined i s data from 2017-10-14 and 2017-10-16 after converting the data from magnitude to flux. For both of the kernels, we found that length scale over which data was correlated was 2.6±0.3 min, as suggested by the ACF in Section 4.4. While a solution was found using the second kernel with the short time scale periodicity, the marginal log likelihood never exceeded the marginal log likelihood of the first kernel. Even when we limited our analysis to just the data taken on 2017-10-16 prior to the eclipse (the same data used when calculating the middle row of panels of Figure 10), we again find that the simpler first kernel is more likely than the second, suggesting that the periodicity detected using the Lomb-Scargle periodogram is transient in nature.
The best fitting A 1 , , and A 2 are given in Figure 11. Due to the very short time-scale correlation in our kernel of 2.7 minutes, the model quickly loses the ability to accurately predict fluxes outside of the observed data, limiting the use of this technique.

X-RAY SPECTROSCOPY
The extracted X-ray spectra covering 0.3-10 keV from the PN and both MOS1 and MOS2 instruments onboard XMM-Newton are shown in Figure 12. We modelled the X-ray emission using Xspec v. 12.10.1 (Arnaud 1996). We initally began with a simple absorbed power law (const × tbabs × powerlaw, where the const component was included to allow for differences between the PN and MOS instruments and tbabs is the Tuebingen-Boulder ISM absorption model; Wilms et al. 2000). While this model was able to describe the hard X-ray tail of the spectrum (energies > 2 keV) it could not fit the spectrum at soft (0.1-2 keV) energies, with a χ 2 value of 756.656 for 260 degrees of  Figure 11. Results of the MCMC analysis of our data using Gaussian Process modelling. The modelling shows that the data has a short term correlation of 2.6 min ( ) with an amplitude of close to 0.1 mJy (A 1 ). The orbital modulation has an amplitude of 0.26 mJy (A 2 ).
freedom. Since this source is edge on, we expect the amount of absorption due to material within the binary system to be high. As such, we next fit a partially absorbed powerlaw to the data to account for the circumstellar absorption. The model (const × tbabs × pcfabs × powerlaw) fit the data well, with a χ 2 = 272.928 for 259 degrees of freedom. The model is shown alongside the spectra in Figure 12, and the best fit parameters are given in Table 2. We also tried fitting the spectrum with an additional Gaussian emission component to account for any emission at the 6.4 keV Fe kα emission line. The addition of such a component did not increase the goodness-of-fit by a statistically significant amount.
The best fit parameters show that the interstellar absorber accounts for very little absorption in the spectrum, with only an upper limit on the density calculable. Nearly all of the absorption arises due to the high density and covering fraction of the partial absorber, indicative of significant absorption by the accretion disc. This is in line with the high inclination derived from the optical light curve modelling.
We computed the unabsorbed 2-10 keV flux by fitting the spectrum with the Xspec model cflux and using the best-fit parameters found for the above model. The unabsorbed 2-10 keV flux was found to be (1.73 ± 0.08) × 10 −12 erg cm −2 s −1 .

Component Masses
For each model generated by ELC, the software also computes the individual masses using the system inclination, mass ratio, orbital period, and K 2 . The primary and secondary masses for every model computed in our analysis are shown in Figure 13. The masses are not very well constrained, with M 1 = 1.43 +0.33 Residuals Figure 12. The X-ray spectra from 0.3 -10 keV obtained using the EPIC-PN, -MOS1, and -MOS2 instruments onboard XMM-Newton, alongside the best-fit partially absorbed power law model. (1.9 ± 0.8) × 10 −4 a photons keV −2 cm −2 s −1 at 1 keV the 1σ level. We can use the lower bound of K 2 > 283 km s −1 and the measured P orb to put a lower bound on the primary mass of M 1 > 0.86 M by using the mass function formula and letting q → ∞ (for q = M 1 M 2 as defined for ELC) and i → 90°. From our modelling, the primary mass is less than 2.5 M at the 3σ level. However, this upper bound is not entirely reliable, as the observed K 2 may be biased by irradiation of the secondary. Figure 13 also shows why there is a large error on the mass ratio from our modelling -for a very low-mass secondary (< 0.18 M ), the mass ratio of the system can be as large as 10. If we require that the secondary mass > 0.18 M , then the mass ratio is more tightly constrained to q = 3.5 ± 1.0.
We note that our constraints on the individual masses are not as tight as those derived by S16. This is because S16 measured the radial velocity feature of the emission lines and assumed that this velocity was equivalent to the velocity of the primary star. This assumes that the accretion disc extends down to the surface of the primary, which is not entirely obvious. Additionally, any brightness asymmetries in the accretion disc would cause an incorrect measurement of K 1 . Due to these two points, we did not include this value in our modelling, meaning the mass ratio was left uncon- . Component masses of system. The upper limit on the primary mass is close to the theoretical limit on the mass of a neutron star. Contours are at the 2σ and 3σ levels.
strained. A direct measurement of K 1 would become possible if the primary were detected as a radio or X-ray pulsar.

Is J0427 a cataclysmic variable?
The mass of the primary in J0427 has still not been well constrained, and the system could still be a cataclysmic variable (an interacting binary system with a white dwarf primary) with a particularly heavy white dwarf primary. Below we lay out the arguments that suggest the primary is not a white dwarf. First, the optical light curve of the system helps rule out a typical dwarf nova classification. This type of cataclysmic variable exhibits outbursts related to instabilities in their accretion discs which have a typical amplitude of several magnitudes at optical wavelengths. The lack of any observed outbursts in J0427 in our data or in the data presented in S16 help immediately rule out this classification.
The observed parallax of J0427 in the recent Gaia 2nd data release is 0.38 ± 0.07 mas (Brown et al. 2018). This implies a distance to J0427 of 1.8 < d < 5.6 kpc at the 3σ level, estimated using the methods described in Luri et al. (2018) and with a length scale of 1 kpc for the prior. The distance to J0427 combined with the 0.3-10 keV Xray flux of 2.0 +15.9 −0.5 × 10 −12 erg cm −2 s −1 measured here suggests that the 0.3-10 keV X-ray luminosity of the source is 7 × 10 32 < L X < 8 × 10 33 erg s −1 . Such a high X-ray luminosity is incompatible with dwarf novae cataclysmic variables, which typically have X-ray luminosities 10 29−30 erg s −1 . However, cataclysmic variables which have persistently high mass accretion rates (called nova-likes) and cataclysmic variables with magnetic white dwarfs have much higher Xray and optical luminosities than their dwarf nova cousins. In particular, intermediate polars (which are magnetic CVs that harbour accretion discs; IPs) have X-ray luminosities of 10 32−33 erg s −1 (Xu et al. 2016), which is consistent with the lower limit on the X-ray luminosity of J0427. Additionally, an IP at same distance as J0427 would have a similar apparent magnitude at optical wavelengths as J0427. However, we note that the above comparisons are lower limits, since we have assumed the minimum distance possible to J0427, and it is likely that the X-ray flux from J0427 is much higher than the lower limit implies (as hinted at by the ELC modelling). The optical polarimetry does not help rule out a cataclysmic variable classification, as the upper limit of 1% on circular polarisation is higher than the detected circular polarisation in most magnetic CVs (Butters et al. 2009).
The most convincing evidence for ruling out a white dwarf primary is the association of J0427 with a GeV γ-ray source, as no cataclysmic variable has ever been observed to produce such high energy photons in their regular states. The only state in which a system with a white dwarf primary has exhibited γ-ray emission is when a CV undergoes a nova eruption due to thermonuclear burning on its surface (for example V407 Cygni; Abdo et al. 2010), and the γ-rays are only produced during the eruption itself. The confirmation of the persistent γ-ray eclipse in J0427 proves that the γrays are coming from close to or directly from the primary in this system, and are not due to a nova eruption.

Is J0427 a transitional MSP?
If J0427 is not a CV, then it is a low-mass X-ray binary (LMXB) containing either a NS or a black hole primary. The detection of radio or γ-ray pulsations from the primary or the detection of an X-ray burst would confirm the existence of a NS primary, while an accurate mass measurement could be used to distinguish between these two primaries. Regarding the first point, J0427 has yet to be detected as a persistent or periodic radio source, and regarding the Xray flaring, the X-ray light curve shows significant variability (even out side times of high background, as shown in Figure 2). In particular, there is a single event where the count rate increases by a factor of 5. However, these events are not resolved with the 100s bin width which is required to obtain an adequate S/N in the light curve, making the actual duration and shape of the events impossible to determine.
As such, we must rely on the mass constraints given in this paper to determine the nature of the compact object. The maximum mass of the primary star is 2.42 M at the 3σ level, which is close to the theoretical upper limit for the mass of a NS. This means the most likely scenario is that the primary is not a black hole, but a NS. Given this, the question becomes is the system a regular LMXB with a NS primary, or a tMSP in an accreting state?
tMSPs have a combination of optical features and X-ray features which can be used to distinguish them from regular LMXBs -optical/X-ray flares and a bi-modality in the optical/X-ray flux distribution (see Bogdanov et al. 2015 for an example of X-ray bi-modality and flares, Shahbaz et al. 2015 for an example of optical bi-modality, and Kennedy et al. 2018 for examples of optical flares). The origin of these features is still uncertain, and while there are a multitude of models which attempt to explain the features (Linares 2014;Papitto & Torres 2015;Papitto et al. 2019), most of them agree that the features likely arise from interactions deep within the accretion disc, close to the magnetosphere of the NS.
J0427 shows no such bi-modality in its optical flux distributions (after removal of the orbital modulation, the flux distributions in all 3 bands were log-normally distributed), or in the X-ray light curve obtained with XMM-Newton (with the caveat that a majority of the X-ray light curve must be ignored due to high background). Such a discrepancy is not at odds with the classification of J0427 as a tMSP, as Kennedy et al. (2018) showed that in long term optical observations of PSR J1023+0038 taken in its accretion states there were long periods of time when the flux of PSR J1023+0038 did not show a bi-modal distribution.
It may also be that our ability to detect bi-modality is inclination dependent. Since the inner part of the disc is where the bi-modality supposedly originates, and J0427 is an edge-on system, the region where this bi-modality arises may be hidden from us. This allows for the classification of J0427 as a tMSP despite the non-detection of a bi-modal optical and X-ray flux distribution.
Additionally, the derived X-ray flux from the XMM-Newton observations and the measured power law index of Γ = 1.3±0.1 are consistent with a tMSP in an accreting state. The distance limits on J0427 combined with the measured 2-10 keV flux provide limits on the X-ray luminosity of J0427 as 7 × 10 32 < L X < 8 × 10 33 erg s −1 . The measured X-ray luminosity of PSR J1023+0038 during its accreting state is 3 × 10 33 erg s −1 (when in the X-ray high mode). If J0427 were to have a similar luminosity, then it would be located at 4 kpc. The next Gaia data release should have a more precise parallax for J0427 which will then better constrain the X-ray luminosity.
Finally, typical accreting LMXBs do not have γray counterparts, while the 3 confirmed tMSPs (PSR J1023+0038, IGR J18245-2452, and PSR J1227-4853) and 1 of the potential tMSPs (3FGL J1544.6-1125) do. While the definite detection of a γ-ray eclipse in J0427 shown in this paper strongly supports the claim by S16 that J0427 is a tMSP in an accreting state, the conclusive confirmation of J0427 as an accretion-state tMSP requires the detection of millisecond pulsations from the spinning neutron star primary.

Pulsation detection feasibility
While radio pulsations have only been detected from tM-SPs in rotation-powered states, optical and X-ray pulsations have been detected from PSR J1023+0038 in its accreting state (Archibald et al. 2015;Jaodand et al. 2016;Ambrosino et al. 2017;Zampieri et al. 2019). Additionally, while γray pulsations have been detected from a large number of nearby MSPs, it remains an open question as to whether or not γ-ray pulsations from a tMSP are also suppressed during accreting states. Of the three confirmed tMSPs, γ-ray pulsations have only been detected from PSR J1227−4853 during its rotation powered state (Johnson et al. 2015), but the low signal-to-noise ratio and unpredictable variations in the orbital period have prevented the extrapolation of the timing solution back to the pre-transition epoch to check for pulsations in the accreting state.
For J0427 whose pulsation period remains unknown, the detection of γ-ray pulsations would require a multi-dimensional blind search over several timing parameters. Given sufficiently precise knowledge of the orbital parameters (period, phase and projected semi-major axis), such searches are capable of discovering binary MSPs (Pletsch et al. 2012). However, for J0427 the current parameter constraints (from the orbital ephemeris of Equation 1, and the semi-major axis range of 0.7 s (a sin i)/c 1.6s inferred from the parameters in Table 1) still leave a prohibitively large volume to search, and even using e.g. the thousands of computers participating in the Einstein@Home volunteer computing project (Allen et al. 2013;Clark et al. 2017), a full blind search would take more than a year to complete (L. Nieder, private communication). Given that it is still unclear whether or not a tMSP in an accreting state will even emit detectable γ-ray pulsations, we do not believe that such a search would currently be a good use of computing resources, but we note that this may change in the future should the orbital ephemeris be further refined with additional optical observations.

CONCLUSIONS
We have presented further optical, X-ray, and γ-ray observations of J0427. The γ-ray eclipse is now detected with a significance > 5σ. This confirms that the γ-ray emission is associated with the X-ray and optical source, establishing J0427's membership as one of a rare class of accreting binaries. The high time resolution optical data show rapid flickering on a timescale of 2.4 mins, with hints of an underlying 21 min period. Modelling of the optical light curve has placed tight constraints on the inclination of the system to be 84 ± 3°, and we find that there is significant evidence for heating of the secondary star by the primary.
We do not find evidence for bi-modality in the optical or X-ray light curvs. This is still consistent with J0427 belonging to the tMSP class of objects, as tMSPs do not always show bi-modality in their optical and X-ray light curves, and our ability to detect bi-modality may be inclinationdependant, with the bi-modality of high inclination systems such as J0427 difficult to detect due to obscuration of the region associated with this behaviour by the outer parts of the accretion disc. While we have not been able to rule out a white dwarf primary by modelling the optical light curve, it is likely that the primary is too heavy to be a white dwarf, and the now significant detection of a γ-ray eclipse makes any primary other than a NS in a tMSP difficult to explain.
Definitive classification of J0427 as a tMSP requires detection of radio/optical/X-ray pulsations from the primary, or the detection of a state transition. One thing that is beyond debate is that additional optical, X-ray, radio, and γray data should be obtained. In addition to searching for state transitions, further optical photometry can be used to confirm the ∼ 21 min period, and optical spectroscopy would allow the use of Doppler tomography and spectral eclipse mapping to better understand the accretion structures within this system. In order to construct accurate confidence levels for peaks in the power spectrum, it is important to know what the underlying shape of the power spectrum is. Typically power spectra are modelled as power laws (P( f ) ∝ f α where α is the spectral index). A flat spectral index (α ∼ 0) suggests that the light curve is dominated by white noise, while a steep spectral index (α ∼ −2) suggests strong correlated red noise is present in the data. The light curves of systems which contain accretion discs often are often dominated by "flickering", which shows up in power spectra as pink noise (a power-law spectrum with α = −1; Lawrence et al. 1987). For each of the power spectra presented in this paper, we first binned the power spectrum in question such that the distribution of power within each bin followed a log-normal distribution, and then fit both a power law and a broken power law to the binned data.
Once the frequency dependence of the power spectrum was measured, the 3σ level for identifying significant periods present in the power spectrum at each frequency could be estimated. This was done by generating 100,000 light curves which had the same temporal sampling as our data, but were generated using a fake power spectrum with either the best fitting power law or a broken power law from the first step, and then creating the light curve using the algorithm described in Timmer & Koenig (1995) and implemented in Stingray 2 . The power spectrum of each of these fake light curves was taken, and the power at each frequency recorded. The distribution of powers at each frequency were then fit with a cumulative distribution function assuming the noise is Gaussian distributed (equation 53 of VanderPlas 2018), and the 3σ level at each frequency was taken to be the central value of the Gaussian plus 3 standard deviations. The threshold line is not the "single trial" threshold, but shows the 3-sigma level after accounting for the number of independent frequencies in the power spectrum. In this case, this number was assumed to be the length of the frequency range of the power spectrum up to the Nyquist frequency ( f Ny = 0.5 t samp where t samp is the sampling rate) divided by a delta frequency such that the power spectrum was Nyquist sampled (δ f = 1/T obs where T obs is the length of the observations). The results of this error estimation are plotted as the 3σ levels in each of the power spectra in this paper, and code to reproduce these plots is hosted online 3 .

A2 Gaussian Process Modelling of the optical light curve
For a pair of data points, (x i , x j ), we tested two kernels: one of the form k(x i , x j ) = (A 1 ) 2 exp − (x i − x j ) 1 + (A 2 ) 2 exp − 2 ( 2 ) 2 sin 2 π(x i − x j ) P or b , and another of the form k(x i , x j ) = (A 1 ) 2 exp − (x i − x j ) 1 + (A 2 ) 2 exp − 2 ( 2 ) 2 sin 2 π(x i − x j ) P or b + (A 3 ) 2 exp − 2 ( 3 ) 2 sin 2 π(x i − x j ) P . (A2) The first term in each kernel is a Matérn covariance function with ν = 0.5 and allows for covariances between data points on a length scale of 1 (Rasmussen & Williams 2006), the second term in each kernel allows for periodic variations on a time scale equal to the orbital period, and the third term in the second kernel was to allow for short-term periodic variations. 2 is defined such that 2 < 1 allows for rough periodic variations, while 2 > 1 lets the variations be more strictly sinusoidal in their appearance (the same is true of 3 ). We fixed 2 and 3 to 1.0 based on the autocorrelation function calculated previously in Section 4.4. In GPM, A 1 , A 2 , 1 , 2 , and P are often referred to as hyper-parameters, as their physical meanings can, at times, be difficult to understand. This paper has been typeset from a T E X/L A T E X file prepared by the author.