© The Authors 2023

1. Introduction

Transitional millisecond pulsars (tMSPs) are rapidly rotating (P ≲ 10 ms), weakly magnetized (≈10⁸ − 10⁹ G) neutron stars (NSs) that have been observed to swing between distinct states within a few days that are likely powered by different physical mechanisms. These sources are part of binary systems with low-mass (M ≲ 1 M_⊙) companion stars. At high-mass accretion rates, an accretion-powered pulsar that is able to channel the inflowing matter toward the NS magnetic poles is observed. When the accretion rate decreases, the pulsar wind sweeps the matter transferred through Roche-lobe overflow away, and a rotation-powered radio pulsar is seen. Variability in the mass inflow rate is the driver of these state changes.

Three confirmed tMSPs are known to date: PSR J1023+0038 (Archibald et al. 2009), XSS J1227−4859 (de Martino et al. 2010; Bassa et al. 2014), and IGR J1824−2452 (Papitto et al. 2013). All of them have also been observed in an intermediate state, called subluminous disk state. The physical mechanism that powers this state still remains to be fully understood. In this state, tMSPs show the presence of an accretion disk and γ-ray emission that is up to ten times more intense than the emission that is observed during the rotation-powered state (Papitto & de Martino 2022, and references therein). This is at variance with low-mass X-ray binaries (LMXBs), which generally do not emit a detectable γ-ray flux. Another peculiarity of tMSPs is the X-ray luminosity (L_X ∼ 10³³ − 10³⁴ erg s⁻¹), which is lower than what is usually observed in the outburst phase of accreting millisecond X-ray pulsars (AMXPs) (L_X ∼ 10³⁶ erg s⁻¹), but higher than the luminosity in the rotation-powered state (L_X <  10³² erg s⁻¹). The X-ray emission is also variable over timescales of a few tens of seconds: different intensity modes (high, low, and flaring modes) have been observed in the X-ray light curves (see, e.g., de Martino et al. 2013; Linares 2014; Bogdanov et al. 2015; Archibald et al. 2015; Coti Zelati et al. 2018). J1023+0038 has been in the subluminous disk state from June 2013 (Patruno et al. 2014; Stappers et al. 2014) until the time of writing (October 2022) and is in high mode for ∼80% of the time (Bogdanov et al. 2015; Archibald et al. 2015; Jaodand et al. 2016). Transitions from high to low mode, and vice versa, occur unpredictably on a timescale of ∼10 s. The duration of these modes varies from a few tens of seconds to a few hours (Papitto & de Martino 2022).

J1023+0038 is the first millisecond pulsar that was ever observed as an optical pulsar (Ambrosino et al. 2017; Zampieri et al. 2019; Karpov et al. 2019; Burtovoi et al. 2020). Optical pulsations have recently also been observed from the AMXP SAX J1808.4−3658 (Ambrosino et al. 2021), but J1023+0038 so far remains the only tMSP with detectable pulsed emission in the optical band. Optical and X-ray pulsations from J1023+0038 were detected simultaneously in the X-ray high modes, but they disappeared when the source transited in the low modes. This suggests a common emission mechanism (Papitto et al. 2019). Recently, pulsations have been observed in the UV band. The UV emission, like the optical emission, undergoes transitions between high and low intensity modes that occur simultaneously with those observed in the X-ray band (Jaodand et al. 2021; Miraval Zanon et al. 2022).

Previous works attempted to determine the physical origin of optical pulsations (see, e.g., Campana et al. 2019; Papitto et al. 2019; Veledina et al. 2019). Individually, the standard rotation- and accretion-powered mechanisms hardly explain the observed optical pulsed luminosity (Ambrosino et al. 2017; Papitto et al. 2019). X-ray pulsations were first interpreted as resulting from the channeling of matter along the magnetic field lines with the subsequent formation of accretion columns at the NS poles (Archibald et al. 2015). However, the optical luminosity of accretion columns is expected to be much lower than the observed optical pulsed luminosity (L_opt ∼ 10³¹ erg s⁻¹; Ambrosino et al. 2017). Even assuming that optical pulses are generated by cyclotron emission from electrons falling into the accretion columns, the expected luminosity would be about 40 times lower than observed (Ambrosino et al. 2017; Papitto et al. 2019). Similar energetic arguments can also exclude the reprocessing of accretion-powered X-ray emission at the surface of the companion star and/or in the outermost regions of the disk, as occurs in some X-ray binaries with a strongly magnetized and slowly rotating accreting pulsar (Ambrosino et al. 2017).

On the other hand, optical emission driven by the rotation of the NS magnetic field would require an efficiency in converting the spin-down power to the pulsed optical emission up to 10⁴ times higher than the values (∼5 × 10⁻⁶ − 2 × 10⁻⁹) measured for the five isolated rotation-powered pulsars from which optical pulses were detected (Cocke et al. 1969; Mignani 2011; Ambrosino et al. 2017). Moreover, the fraction of spin-down power converted into X-ray pulses would be much higher than that of almost all rotation-powered pulsars (Papitto et al. 2019). Consequently, the rotation-powered mechanism by itself cannot be the common origin of X-ray, UV, and optical pulsations from J1023+0038.

These implications led to the proposal of the shock-driven mini pulsar nebula scenario. Optical and X-ray pulsations originate from synchrotron emission in a shock that forms beyond the light cylinder radius, where the striped pulsar wind meets matter from the inner accretion disk (Papitto et al. 2019; Veledina et al. 2019). In this region, electrons are accelerated to relativistic speeds and emit synchrotron radiation by interacting with the magnetic field in the shock region. This configuration permits a higher fraction of the spin-down energy to be converted into X-ray pulses compared to the previously discussed cases. Pulsar wind nebulae indeed radiate up to a few percent of the pulsar spin-down power (Kargaltsev & Pavlov 2008, 2010; Vink et al. 2011; Torres et al. 2014).

Analyses of simultaneous observations performed in May 2017 with XMM-Newton and the fast optical photometer SiFAP2, mounted at the INAF Telescopio Nazionale Galileo (TNG), found that optical pulses lag the X-ray pulses by ∼200 μs (Papitto et al. 2019). The proposed model interprets this time lag in terms of the different timescales that synchrotron X-ray and optical photons take to be emitted. However, the above measurement was affected by the absolute timing accuracy of SiFAP2 (∼60 μs; Papitto et al. 2019) and by that of XMM-Newton/EPIC.

In a Calibration Technical Note of May 2022¹, the uncertainties on arrival times acquired through XMM-Newton/EPIC were reviewed and reached a value of 100 μs for the timing mode, which is more than twice that calculated by Martin-Carrillo et al. (2012) of 48 μs and considered by Papitto et al. (2019) to estimate the significance of the optical/X-ray pulse phase lag they measured. The absolute timing accuracy of these instruments makes the systematic error associated with time lag estimated in Papitto et al. (2019) compatible with the measure itself. This highlights the importance of presenting an in-depth study of the relation between optical and X-ray pulsations that also analyzes simultaneous observations acquired only with NICER and the fast optical photometer Aqueye+, mounted at the Copernicus Telescope in Asiago. The NICER absolute timing accuracy is estimated to be < 300 ns², while that of Aqueye+ is < 0.5 ns (Zampieri et al. 2015). We report here a detailed timing analysis performed on optical/X-ray simultaneous or quasi-simultaneous observations to elucidate the physical mechanisms that cause the observed pulsations. Our data were acquired with the XMM-Newton and NICER X-ray satellites and with the fast optical photometers SiFAP2 and Aqueye+ over a time ranging from May 2017 to the beginning of February 2022.

Section 2 is dedicated to the description of the observations and the data processing techniques. In Sect. 3, we perform the phase analysis of simultaneous or quasi-simultaneous observations in the optical and the X-ray bands. By modeling the pulse profiles as the sum of two harmonic components, we study the time lags between the pulsations in the two different observational bands. We discuss our results in Sect. 4 and constrain the synergistic model proposed to explain the emission mechanisms of optical and X-ray pulsation from J1023+0038. Last, in Sect. 5, we summarize our main results and outline future prospects.

2. Observations

Table 1 lists the observations analyzed in this paper. They were selected with the aim of studying (quasi-)simultaneous observations of J1023+0038 in the optical and X-ray bands in a time interval of about five years. In the following, we detail the analysis of the different data sets.

2.1. X-ray observations

2.1.1. NICER

We present the analysis of NICER observations of J1023+0038 that were simultaneously or quasi-simultaneously performed with optical observations from January 2019 to February 2022. The events were reduced and processed using HEASoft version 6.28 and NICERDAS version 7a. We corrected the photon arrival times to the Solar System Barycenter (SSB) using the JPL ephemerides DE405. We adopted the source coordinates RA (J2000) = 10:23:47.687198(2) and Dec (J2000) = +00:38:40.84551(4) (Deller et al. 2012). We use these coordinates throughout the rest of the work. We estimated the background contributions to our data with the tool nibackgen3C50 (Remillard et al. 2022).

Generally, it was sufficient to analyze the individual observations listed in Table 1 to derive an orbital solution (Table 2), except in the case of the June 2019 data set, for which the close-by observation of 2019 June 13 (Obs. ID: 2531010401) had to be used as well.

2.1.2. XMM-Newton/EPIC

The XMM-Newton observations were performed on 2017 May 23 and 24 (presented in Papitto et al. 2019), and on 2018 December 11 and 15 (Table 1). The data were reduced using the Science Analysis Software (SAS) v.16.1.0. In each observation, the EPIC-pn was operating with a time resolution of 29.5 μs (timing mode) and a thin optical blocking filter. The photon arrival times observed by XMM-Newton were reported to the SSB, using the JPL ephemerides DE405 and the barycen tool from HEASoft. We defined source and background regions with coordinates RAWX = 27–47 and RAWX = 3–5, respectively, and retained good events characterized by a single or a double pattern.

2.2. Optical observations

2.2.1. TNG/SiFAP2

We analyzed the optical observations reported in Table 1 acquired with the SiFAP2 fast optical photometer (Meddi et al. 2012; Ambrosino et al. 2016, 2017) mounted at the TNG. The arrival times of each photon were referred to the SSB through the TEMPO2 package (Hobbs et al. 2006), using the JPL ephemerides DE405.

In the January 2019 data set, the statistics was too poor to obtain precise measurements of the epoch of passage at the ascending node, T_asc, and of the spin period, P (Table 2). In this part of the analysis, we merged that observation with SiFAP2 observations carried out on 2019 February 1, with a total exposure of ∼10 ks.

The SiFAP2 quartz clock is characterized by drifts with respect to the actual time measured by two global positioning system (GPS) pulse-per-second (PPS) signals that are used to mark the beginning and end of each observation. Following Ambrosino et al. (2017), the arrival times recorded by SiFAP2, t_SiFAP2, were corrected assuming the linear relation t_arr = t_SiFAP2 × (Δt_GPS/Δt_SiFAP2), where Δt_GPS and Δt_SiFAP2 are the total elapsed time measured by the GPS and SiFAP2 clocks, respectively. In January 2019, we had Δt_SiFAP2 − Δt_GPS = −1.164 ms, with the value of Δt_GPS reported in Table 1. During the observations of June 2019, we had Δt_SiFAP2 − Δt_GPS = −2.323, −1.156 ,and − 1.511 ms. Finally, for January 2020 observations, we had Δt_SiFAP2 − Δt_GPS = +0.580, +0.751, +0.822, and +0.663 ms.

2.2.2. Copernicus/Aqueye+

Aqueye+ is an ultra-fast optical single photon counter mounted at the Asiago 1.8-m Copernicus Telescope with the capability of time-tagging the detected photons with subnanosecond time accuracy (Zampieri et al. 2015). The chosen Aqueye+ observations (Table 1) were reduced with the QUEST software (v. 1.1.5, see Zampieri et al. 2015). The arrival times of each photon were referred to the SSB through the TEMPO2 package, using the JPL ephemerides DE405.

3. Data analysis

For each data set, we first corrected the photon arrival times for the pulsar orbital motion in the binary system. We set the orbital period and the projected semimajor axis equal to the values found in the timing solution of Jaodand et al. (2016) (see Table 2), and we performed a search on the epoch of passage at the ascending node, T_asc. We used a grid of T_asc values spaced by 0.37 s (Caliandro et al. 2012) around the estimate extrapolated from Jaodand et al. (2016). We carried out an epoch-folding search on each time series by sampling each spin period, P, with 16 phase bins. The final best T_asc was determined by fitting the peak of the χ² distribution with a Gaussian function. We considered the half width at half maximum (HWHM) of the Gaussian as the uncertainty to associate with T_asc. To improve the spin period estimate obtained from the T_asc search, we corrected the photon arrival times with the best-fitting values of the orbital parameters and performed an epoch-folding search with 16 phase bins. The best P value was then estimated by modeling the peak of the χ² distribution with a Gaussian function. The associated error was calculated using Eq. (6a) from Leahy (1987). Table 2 summarizes the best-fitting values of T_asc and P found for each data set.

We verified that the timing results were compatible between simultaneous X-ray and optical observations. This allowed us to use the values of T_asc and P found from the X-ray timing (which are more accurate than those obtained from the analysis of optical observations due to the higher root mean square (rms) pulse amplitude of the signal, i.e., the pulse amplitude divided by the square root of 2) to correct the photon arrival times for the pulsar orbital motion and to perform the phase analysis of simultaneous optical and X-ray observations. Since the folded pulse profiles are double-peaked (see, e.g., Archibald et al. 2015; Ambrosino et al. 2017; Papitto et al. 2019), we modeled them using a decomposition function with two harmonic terms,

$$\begin{array}{c}\hfill F(\varphi )=K\{1+{\displaystyle \sum _{i=1}^2}{r}_{i}sin[2\pi i(\varphi -{\varphi }_i)]\},\end{array}$$(1)

where K is the average count rate, and the free parameters, r_i and ϕ_i (i = 1, 2), are the fractional amplitude and the phase of the two harmonics, respectively. Uncertainties of our best-fitting values were estimated from the parameter range required to increase the χ² from the fit of a quantity Δχ²(α = 68%) = 1.0 (Lampton et al. 1976; Avni 1976; Yaqoob 1998).

3.1. Evolution of the time of passage at the ascending node

Figure 1 shows the difference ΔT_asc between our values of T_asc in NICER observations (Table 2) and that computed using the radio timing solution (Archibald et al. 2013; Jaodand et al. 2016) as a function of the number of orbital cycles since T_ref = 57897.027668 MJD. A similar approach was adopted in previous works (see, e.g., Jaodand et al. 2016; Papitto et al. 2019; Burtovoi et al. 2020). NICER observations in addition to those chosen in this paper simultaneously with optical data (Table 1) were analyzed to study the T_asc long-term evolution (Table 3). We selected observations with an exposure > 10 ks in order to have good statistics for the timing analysis and at least two measurements of T_asc per year.

Fig. 1. Long-term evolution of Tasc as a function of the number of orbital cycles since Tref = 57897.027668 MJD. Red points are the values found by Jaodand et al. (2016), which deviate from the almost linearly increasing trend found in this work. Blue points are from Papitto et al. (2019), and lighter blue points are from Burtovoi et al. (2020). Triangles, associated in the legend with the different NICER observations, are from this work. The thick dashed line indicates our increasing roughly linear trend. The dotted gray line indicates the sinusoidal trend that is similar to the trend found in previous work (e.g., Jaodand et al. 2016; Papitto et al. 2019).

We found an increasing trend of ΔT_asc with time, as already inferred since May 2017 by Burtovoi et al. (2020), who emphasized that the steady increase may indicate a systematic underestimation of the orbital period of the system. We modeled the second part of the data in Fig. 1, that is, from T_ref = 57897.027668 MJD, with the following expression:

$$\begin{array}{c}\hfill \mathrm{\Delta }{T}_{\mathrm{asc}}(N_{\mathrm{orb}})=A+B{N}_{\mathrm{orb}}+\frac{1}{2}{P}_{\mathrm{orb}}C{N}_{\mathrm{orb}}^2.\end{array}$$(2)

P_orb is the orbital period, which was fixed at the value given in Table 2, while A, B, and C are free parameters. The integer number of orbital cycles since T_ref is N_orb = int[(T − T_ref)/P_orb]. We obtained A = (− 27.3 ± 0.8) s, B = (5.0 ± 0.4) × 10⁻³ s, and C = (1.7 ± 0.5) × 10⁻¹¹ s s⁻¹. In general, redbacks, that is, millisecond pulsars in a close orbit with a low-mass companion star such as the source under consideration, show unpredictable variations of the orbital phase (see, e.g., Jaodand et al. 2016, Fig. 2). Therefore, we caution that the B term might be related to the correction we should apply to the orbital period to have constant residuals over time and C as the orbital period derivative. Since the orbital period difference would not significantly affect the results of the analysis presented here, we retained the value reported by Jaodand et al. (2016). Future studies will allow us to firmly establish the binary evolution.

3.2. Nonselection of intensity modes

Except for the May 2017 data set, we did not distinguish among low, high, and flaring modes, differently from what was done in other works (see, e.g., Bogdanov et al. 2015; Papitto et al. 2019). As we show in the following, the nonselection of modes does not change the phase of our signal. The reason for this choice lies in the fact that different modes are typically identified in the X-ray band: just a fraction of the optical data considered in this work is strictly simultaneous to X-ray observations, differently from the case of the 2017 observations analyzed in Papitto et al. (2019, hereafter P19). It is not possible to precisely distinguish high and low modes in the optical band because the corresponding variations in intensity are fainter than those observed in the X-rays and are preferentially observed in the red part of the visible spectrum (see, e.g., Shahbaz et al. 2015, 2018; Kennedy et al. 2018; Papitto et al. 2018). To verify whether the pulse phase changes depend on either selecting or not selecting the different modes, we analyzed the optical observation acquired in May 2017 (Table 1) in three different ways: by selecting only the high modes; by using the whole observation, therefore without distinguishing between low, high, and flaring modes; and by visually selecting and removing flares, hence keeping both low and high modes. We first selected high modes only in simultaneous intervals with the X-ray observation, adopting the definition of X-ray modes from Bogdanov et al. (2015). Because our values of T_asc and P (Table 2) are compatible with those used by P19, equal to P_ref = 1.687987446019 ms and T_asc = 57896.82926(1) MJD, we used the latter values for a better comparison. The first and second harmonics of the optical pulsation lag the X-ray harmonics by δτ₁ = (223 ± 31) μs and δτ₂ = (195 ± 7) μs (Table 5), respectively. Our results are compatible with those of P19. We note that the phases of the optical pulse profile obtained by selecting only the high modes are ϕ₁ = 0.4726 ± 0.0053 and ϕ₂ = 0.5500 ± 0.0023 (Table 4). Second, analyzing the entire optical observations of May 2017 without selecting the different intensity modes, we obtained ϕ₁ = 0.4714 ± 0.0070 and ϕ₂ = 0.5615 ± 0.0034. These differ by 0.2σ and 5σ, respectively, from the phases estimated when only high-mode intervals were selected. When we removed the flaring intervals (i.e., analyzing both high and low modes), the results are compatible with those obtained by selecting only high modes. In this third case, we indeed obtained ϕ₁ = 0.4666 ± 0.0064 and ϕ₂ = 0.5501 ± 0.0029, both compatible within 3σ with the results from the data set in which only the high modes were analyzed. This is expected because optical (and X-ray) pulses are not detected during the low modes, whereas optical pulsations appear during flares, although the pulse amplitude is six times smaller than in the high modes (Papitto et al. 2019). For this reason, we paid special attention to removing flaring intervals in the remainder of this work, but we did not distinguish between high and low modes. Optical flares are indeed easy to identify (see, e.g., Bogdanov et al. 2015; Papitto et al. 2018, 2019). As an example, Fig. 2 shows the visual selection of flares in the Aqueye+ observation of January 2020. We conclude that, with these caveats, the decision not to distinguish high and low modes does not produce different results in the phase analysis.

Fig. 2. Aqueye+ light curve observed in January 2020, binned every 10 s. The count rate is normalized at the maximum count rate. The vertical colored region indicates the visually identified flaring interval.

3.3. Phase analysis results

We summarize our main results in Tables 4 and 5. We filtered X-ray and optical observations to analyze simultaneous intervals. The cases in which this was not possible are marked in the tables with an asterisk. The reason was either the absence of exact simultaneity or the short duration of these intervals, which made the statistics insufficient to detect the pulsed profile. Table 4 lists the best values for the phases and amplitudes, ϕ_i and r_i, with i = 1, 2, obtained by modeling the pulse profiles with a function consisting of two harmonics (Eq. (1)). T_ref denotes the reference epoch (the same for simultaneous optical/X-ray observations) against which we calculated the phases. Fractional amplitudes were converted to background-subtracted rms amplitudes, that is, $r^{\prime }=(1/\sqrt{2})r{c}_{\mathrm{tot}}/(c_{\mathrm{tot}}-c_{\mathrm{bkg}}),$ where r is the pulse amplitude, c_tot is the total count rate, and c_bkg is the count rate associated with the sky background. $R=(r_1^2+r_2^2{)}^{1/2}$ is the total rms amplitude. We associated the statistical errors computed through the least-squares method, setting Δχ²(α = 68%) = 1.0 (Sect. 3). The absolute timing accuracies for XMM-Newton and SiFAP2 are ∼100 μs and ∼60 μs, respectively, while systematic errors are negligible for NICER and Aqueye+. Table 5 shows the corresponding lags between optical and X-ray pulsations.

Figure 3 focuses on the results from the second harmonic, whose power spectral densities are more than three times higher than those of the first harmonic. Colored error bars represent 1σ statistical uncertainties, and the black error bars indicate the total error. We note the influence of the absolute timing accuracy of XMM-Newton (∼100 μs) and SiFAP2 (∼60 μs). The time lags of optical and X-ray pulsations are always within the range (0 − 250) μs, that is, a phase lag of δϕ ∈ (0 − 0.15), even in observations acquired over five years.

Fig. 3. Time lags relative to the second harmonic term (Table 5). The y-axis is not temporally equispaced. The purple point indicates the values for simultaneous observations between XMM-Newton and SiFAP2, blue points show values between XMM-Newton and Aqueye+, light blue points indicate values between NICER and SiFAP2, and green points show the values between NICER and Aqueye+. The dashed line indicates a zero time lag. Colored error bars represent 1σ statistical uncertainties, while the black error bars indicate the total error. The influence of the absolute timing accuracy of XMM-Newton (∼100 μs) and SiFAP2 (∼60 μs) is visible.

3.3.1. NICER and Aqueye+ simultaneous observations

In this section, we discuss the results of (quasi-)simultaneous observations between Aqueye+ and NICER. They are most valuable for our analysis because the absolute timing uncertainties are negligible compared with the effects we aim to measure. They were also taken with a different set of instruments than in P19, thus with disconnected systematics and of a much smaller magnitude. We have simultaneous observations between these two instruments in February 2019, January 2020, and January–February 2022 (Table 1).

From January 29 to February 2, 2022 an observational campaign was carried out with Aqueye+ and NICER. However, optical data were affected by the bad weather conditions, and the X-ray statistics was often low due to short NICER exposure. The top and middle panels of Fig. 4 show the X-ray and optical rms amplitudes in February 2019, in January 2020, and over the five-day observational campaign in 2022. When optical pulsations were not detected, we estimated upper limits on the pulse amplitude by computing the Fourier power spectral density and measuring the power of the first and second harmonics of the spin frequency. We then converted these power pairs into rms amplitudes at 3σ confidence level according to the procedure described by Vaughan et al. (1994) to take into account that the probability distribution of total power in a frequency bin of a Fourier spectrum containing both signal and noise is more complicated than a χ² distribution. The bottom panel shows the phase of the two harmonic components of the corresponding optical and X-ray pulse profiles. In the main analysis of this work (see Tables 4 and 5, and Fig. 3), we report the results from the first day of the 2022 observational campaign (January 29), when we had the longest X-ray observation that was suitable to provide statistically acceptable results and the weather conditions were good enough to detect optical pulsations. The phase of the second harmonic term of optical pulses lags that in the X-rays by δτ₂ = (108 ± 36) μs.

Fig. 4. Background-subtracted rms amplitudes of X-ray (top panel) and optical (middle panel) pulse during NICER/Aqueye+ simultaneous observations in February 2019, January 2020, and January–February 2022. The arrows represent the upper limits converted into rms amplitudes at 3σ confidence level computed when the optical pulsations were not detected. In the lower panel, we show the phases of the optical pulse (red dots and empty yellow squares show the first and second harmonic, respectively), and the phases of the X-ray pulse (blue dots and light blue empty squares show the first and second harmonic, respectively) during these simultaneous observations.

Special attention was paid to filtering out flare intervals present in the optical data of February 2019 and January 2020. We obtained time lags from the second harmonic component of δτ₂ = (254 ± 47) μs and (119 ± 31) μs, respectively (Table 5).

Figure 5 shows the pulsed profiles of the NICER/Aqueye+ simultaneous observations discussed here.

Fig. 5. Pulse profiles obtained by folding with 16 phase bins of simultaneous observations with NICER (blue points, left y-axis scale) and Aqueye+ (red points, right y-axis scale) in February 2019 (top panel), January 2020 (middle panel), and January 2022 (bottom panel). Solid lines indicate the best-fitting two-harmonic functions. The zero of the pulse phase was shifted arbitrarily in each panel by the same amount for both instruments.

Among the three sets of NICER/Aqueye+ (quasi-) simultaneous observations discussed in this section, the February 2019 data feature the longest interval of simultaneity between NICER and Aqueye+ (Table 4), combined with the highest significance of the pulse profiles. Therefore, they can be considered most valuable for estimating the optical/X-ray time lag. This data set is also of particular interest because it shows the largest time lag with respect to the second harmonic term (Fig. 3). To strengthen the time lag measure obtained from this data set, we performed a cross-correlation analysis of the two pulse profiles as a function of the optical profile offset against that in X-rays. The cross-correlation function (Derrick & Thomas 2004) was corrected for the counting statistics, and the data were wrapped around themselves to avoid leaving unmatched points. For instance, when we shifted the optical pulse profile of one phase bin to the right, the number of counts in the last bin of the optical profile was matched with the number of counts of the first bin of the X-ray profile. We folded the Aqueye+ and NICER pulse profiles of February 2019 with a 64 phase bins. Figure 6 shows that the maximum value of the correlation coefficient is at a lag of 9 phase bins. When 64 phase bins are used, each of them has a length equal to P/64 ∼ 26 μs, where P is the spin period from the X-ray timing in μs (Table 2). Therefore, this phase lag corresponds to a time lag of δτ ∼ (237 ± 13) μs, where the associated error is half the bin width. This result is significantly different from zero at more than 5σ and is fully compatible with the analysis result for the second harmonic phase. This conclusion also shows that the second component prevails in the computation of the correlation between the optical and X-ray pulse profiles.

Fig. 6. Cross-correlation of the pulse profiles from NICER and Aqueye+ observations of February 2019 as the lag varies. The top and middle panels show NICER and Aqueye+ pulse profiles with 64 phase bins, respectively. Solid lines indicate the best-fitting two-harmonic function. The bottom panel shows the cross-correlation coefficient as a function of the lag with which the optical pulse profile is shifted. We used phase lags ranging from −32 to +32, including the zero lag (dashed vertical line).

3.3.2. Other observations

We report the results of the remaining data sets, including observations carried out with XMM-Newton and/or SiFAP2, which must be taken with caution as they are affected by larger absolute timing uncertainties than those of NICER and Aqueye+. Simultaneous or quasi-simultaneous observations performed between XMM-Newton and Aqueye+ in December 2018 and between NICER and SiFAP2 in January 2019 and in January 2020 return time lags that are compatible with the absence of a phase shift considering the systematic errors. The observational campaign of June 2019 performed with SiFAP2 and NICER provides δτ₂ = (124 ± 34) μs, considering the statistical errors alone. When the SiFAP2 absolute timing accuracy is added in quadrature, we obtain δτ₂ = (130 ± 69) μs. An improvement in the absolute timing accuracy of XMM-Newton and SiFAP2 would lead to statistically stronger results. Although the systematic errors affecting these estimates often make the lags compatible with zero, these observations are important to confirm that the time lag relative to the second harmonic term between optical and X-ray pulses is always in the ∼(0 − 250) μs range, that is, is not randomly distributed.

4. Discussion

This paper presented a detailed timing analysis of optical/X-ray (quasi-)simultaneous observations of the tMSP PSR J1023+0038, focusing on the study of the time lags between the pulses in the optical and X-ray bands. We folded the data at the spin periods found in the X-ray timing, compatible with the values estimated in the simultaneous optical observations (Table 2), and we modeled the pulse profiles with two harmonic terms.

The optical pulses have total rms pulsed amplitudes of ∼0.1 − 0.8%, while the X-ray total rms pulsed amplitudes are in the range 4.3 − 10.8% (Table 4), in agreement with what was found in previous works (Archibald et al. 2015; Bogdanov et al. 2015; Ambrosino et al. 2017; Papitto et al. 2019; Zampieri et al. 2019; Miraval Zanon et al. 2022).

Moreover, the rms pulsed fraction is variable over time (e.g., top and middle panels of Fig. 4), as was also found in Miraval Zanon et al. (2022).

Although the estimated time lags (Table 5) are not consistent with being modeled by a single value, we report the weighted averages obtained from the studies of the first and second harmonic terms, respectively: ${\overline{\delta \tau }}_1=(175\pm 22)\mathrm{\mu }\mathrm{s}$ and ${\overline{\delta \tau }}_2=(162\pm 6)\mathrm{\mu }\mathrm{s}$. These results are compatible with each other within 1σ.

By focusing on the second harmonic of the pulse profiles (because their power spectral densities are higher than those of the first harmonic), we found that the time lag between optical and X-ray pulsations always lies in the limited range of (0 − 250) μs, also taking into account observations acquired over about five years (Fig. 3). We note the large total errors due to the absolute timing accuracy of XMM-Newton (∼100 μs) and SiFAP2 (∼60 μs), highlighting the value of the results from NICER/Aqueye+ simultaneous observations.

A timescale of ∼250 μs is virtually equal to the light cylinder radius light-crossing time (for a NS spin frequency of ν ≃ 592.42 Hz, R_LC ≃ 80 km), indicating that the optical pulsed luminosity might be produced by the reprocessing of accretion-powered X-ray emission. However, the observed pulsed luminosity in the visible band is too high to be explained in this way (Ambrosino et al. 2017; Papitto et al. 2019). The properties of optical and X-ray pulsations, such as their simultaneous detection during the X-ray high modes and their disappearance in the low modes (Papitto et al. 2019), the similar pulse shape, and the very limited range of the estimated time lags might indicate that they are related, either by originating in the same region or in a region that lies very near, or that they are connected to the same emission process.

X-ray pulsations observed in the high modes of J1023+0038 were at first interpreted as due to the accretion of matter channeled along the magnetic field lines of the pulsar (Archibald et al. 2015). However, this model fails to explain the observed optical pulsed luminosity of ≈10³¹ erg s⁻¹. When we assume that optical pulsations result from cyclotron emission by electrons that are in-falling in the accretion columns on the NS hotspots, the luminosity produced in this way is about 40 times lower than the observed one (Ambrosino et al. 2017; Papitto et al. 2019).

On the other hand, indicating the efficiency in converting the spin-down power into the optical pulsed luminosity as ${\eta }_{\mathrm{opt}}=L_{\mathrm{opt}}/{\dot{E}}_{\mathrm{sd}}$, where ${\dot{E}}_{\mathrm{sd}}$ is the pulsar spin-down power, the measured values for the five rotation-powered pulsars from which optical pulses were detected (Cocke et al. 1969; Mignani 2011) are in the range from η_opt ∼ 5 × 10⁻⁶ down to ∼2 × 10⁻⁹ (see red points in Fig. 3 from Ambrosino et al. 2017). These NSs are all young (10³ − 10⁵ years), isolated, and possess high-magnetic fields (> 10¹² G). A higher value, equal to η_opt ∼ 2 × 10⁻⁵, is found for J1023+0038. In addition, the fraction of the spin-down power converted into X-ray pulses would be much higher than that of almost all rotation-powered pulsars (Papitto et al. 2019). The rotation-powered mechanism can hardly be the only common source of X-ray and optical pulsations from J1023+0038.

A model was thus proposed in which, despite the presence of an accretion disk, a rotation-powered pulsar is active in the system. In the shock-driven mini pulsar nebula scenario (Papitto et al. 2019; Veledina et al. 2019), optical and X-ray pulses are produced via synchrotron emission from a shock that is formed where the pulsar wind meets the accretion disk, within a few light cylinder radii from the pulsar. Inside the light cylinder, the magnetic field is described in terms of a split-monopole (Bogovalov 1999). The two monopoles with opposite sign give rise to a current sheet on the equatorial plane, which expands as an Archimedean spiral until it reaches the shock in two opposite spots where electrons are accelerated to relativistic speeds. Indicating the shock distance from the pulsar as ∼k R_LC, for low values k ≃ 1 − 2, just beyond the light cylinder radius, the magnetic field is still so intense that synchrotron emission provides the dominant cooling mechanism for shock-accelerated electrons. This model is compatible with the presence of two harmonics. For quite large inclination angles, the emission from the spot on the observer’s side is obscured by the disk. In contrast, the emission from the farthest spot, which is modulated sinusoidally during its rotation, is more easily observed (see Fig. 14 of Papitto et al. 2019). The higher intensity of the second compared to the first harmonic may derive from an asymmetry of the system that causes them to be obscured in a different way.

We used a parametric value of the time lag between optical and X-ray pulsations to estimate some physical quantities of the system within this model. Based on our most accurate results, that is, those from NICER/Aqueye+ simultaneous observations (Sect. 3.3.1), we adopted δτ_p ∼ 150 μs. Synchrotron photons are emitted by the shock-accelerated electrons on a timescale (Eq. (7) of Papitto et al. 2019)

$$\begin{array}{cc}\hfill t_{\mathrm{sync}}& \simeq \frac{\gamma m_{\mathrm{e}}{c}^2}{P_{\mathrm{sync}}}\simeq \frac{9}{4}\frac{m_{\mathrm{e}}^3{c}^5}{e^4{B}_{\mathrm{s}}^2\gamma }\hfill \\ \hfill & \simeq 2.2(\frac{ϵ}{10\mathrm{keV}}{)}^{-1/2}(\frac{B_{\mathrm{s}}}{4.5\times {10}^5\mathrm{G}}{)}^{-3/2}\mathrm{\mu }s.\hfill \end{array}$$(3)

The electron energy is γm_ec² and P_sync ≃ 4/3 σ_Tcγ²U_B is the average synchrotron power per relativistic electron in a source with an isotropic pitch-angle distribution. Here σ_T is the Thomson scattering cross section, U_B is the magnetic energy density, and we used the typical energy of synchrotron photons written as ϵ = ℏω_sync = 3ℏeB_sγ²/(2m_ec). The magnetic field intensity in the post-shock region is given by (Arons & Tavani 1993)

$$\begin{array}{c}\hfill B_{\mathrm{s}}=3(\frac{\sigma }{1+\sigma }{)}^{1/2}(\frac{\dot{E}}{c{f}_{\mathrm{p}}{r}^2}{)}^{1/2}\simeq 4.5\times {10}^5{k}^{-1}{f}_{\mathrm{p}}^{-1/2}\mathrm{G},\end{array}$$(4)

where $\dot{E}=4.43\times {10}^{34}\mathrm{erg}{\mathrm{s}}^{-1}$ is the total spin-down power (Archibald et al. 2013) and r = k R_LC is the shock distance. The magnetization parameter σ (Kennel & Coroniti 1984; Bogdanov et al. 2011) is ≫1 close to the light cylinder where the wind is released (Arons et al. 2002), and f_p represents a geometric factor indicating the sky fraction in which the pulsar wind is emitted (f_p = 1 for an isotropic pulsar wind; Bogdanov et al. 2011). We assumed 1 ≲ k ≲ 2, that is, the shock region where optical and X-ray pulses are produced is just beyond the light cylinder (for k <  1 the matter of the disk would enter the light cylinder, preventing the formation of a relativistic wind from the pulsar).

From Eq. (3), we found that it takes

$$\begin{array}{c}\hfill t_{\mathrm{sync}}(ϵ=5\mathrm{keV})\simeq 3(\frac{B_{\mathrm{s}}}{4.5\times {10}^5\mathrm{G}}{)}^{-3/2}\mathrm{\mu }s\end{array}$$(5)

to emit X-ray photons. We assumed the average energy of NICER photons to be ϵ ∼ 5 keV³. On the other hand, the synchrotron timescale for optical photons, with ϵ ∼ 1 eV, is ∼70 times longer, such that the time lag between optical and X-ray pulsations is

$$\begin{array}{c}\hfill \delta \tau \simeq 220(\frac{B_{\mathrm{s}}}{4.5\times {10}^5\mathrm{G}}{)}^{-3/2}\mathrm{\mu }s.\end{array}$$(6)

Similarly, we can predict the time lags with UV and near-infrared (NIR) pulsations from J1023+0038. Assuming an energy of ϵ ∼ 5 eV for UV photons, we would have a time lag of ∼122 μs between optical and UV pulsations, and of ∼95 μs between UV and X-ray pulses. For a NIR observation in the K band, the photon energy is ϵ ∼ 0.6 eV: NIR pulsations would lag optical pulsations by ∼64 μs, UV pulsations by ∼186 μs, and X-ray pulsations by ∼281 μs.

When we use our parametric value of the time lag, δτ_p ∼ 150 μs, and invert Eq. (6), the intensity of the magnetic field behind the shock is

$$\begin{array}{c}\hfill B_{\mathrm{s}}\simeq 5.8\times {10}^5\delta {\tau }_{\mathrm{p},150}^{-2/3}\mathrm{G},\end{array}$$(7)

where δτ_p, 150 is in unit of 150 μs. For a surface magnetic field of B_surf ≃ 9.6 × 10⁷ G (Deller et al. 2012), we computed the dipolar decrease at the light cylinder radius as B_LC ≃ B_surf (R_LC/R_NS)⁻³, assuming a NS radius of about R_NS ≃ 10 km. Since the magnetic field in the post-shock region is three times the intensity before the shock (Arons & Tavani 1993), we found a value similar to that in Eq. (7).

Combining Eqs. (4) and (7), we obtained

$$\begin{array}{c}\hfill f_{\mathrm{p}}\simeq (\frac{B_{\mathrm{s}}}{4.5\times {10}^5\mathrm{G}}{)}^{-2}\frac{1}{k^2}\lesssim 0.6\delta {\tau }_{\mathrm{p},150}^{4/3}.\end{array}$$(8)

A geometric factor < 1 indicates that the pulsar wind is concentrated on the equatorial plane, as also suggested by several models of the Crab pulsar (see, e.g., Kirk 2006).

Variations in disk truncation radius, that is, in the k parameter, are to be expected as they were found in several magnetohydrodynamic simulations (see, e.g., Parfrey & Tchekhovskoy 2017, Fig. 2, panel d). They clearly show that this is a highly turbulent system. This variability has been proposed to explain the disappearance of optical and X-ray pulsations during the low-intensity modes, when the inner disk boundary is pushed farther out by the pulsar wind (Campana et al. 2019; Papitto et al. 2019). The disk truncation radius may vary also during the high-mode intervals, with smaller changes involved to avoid transitions between different intensity modes while causing variations in pulse amplitude (Miraval Zanon et al. 2022) and phase. Combining Eqs. (4) and (6), we found that the expected time lag between optical and X-ray pulsations scales with the k parameter as $\delta \tau \simeq 220k^{3/2}{f}_{\mathrm{p}}^{3/4}\mathrm{\mu }\mathrm{s}$. If k oscillates between 1 and 2 (i.e., the truncation radius of the disk oscillates between ∼R_LC and ∼2R_LC), the time lag would vary by a factor of ∼2.8. This, as well as shock height changes, may explain some of our measured variations in the time lags (Table 5). Another reason may be that some data sets do not have intervals of exact simultaneity between optical and X-ray observations. There may also be some dependence on source parameters, such as luminosity and orbital phase, which it will be necessary to investigate in future works.

J1023+0038 and the AMXP SAX J1808.4−3658 are the only two millisecond pulsars with detectable optical pulsations so far (Ambrosino et al. 2017, 2021; Papitto et al. 2019). Interestingly, optical pulsations observed from SAX J1808.4−3658 during an accretion outburst in August 2019 seem to be almost in antiphase with those in the X-rays (Ambrosino et al. 2021), as opposed to what was observed in J1023+0038. Therefore, optical and X-ray pulsations from SAX J1808.4−3658 can hardly be explained by the same physical mechanism. An important step to confirm the origin of optical and X-ray pulsations from J1023+0038 and test the proposed model may come from the observation of other transitional millisecond pulsars in the subluminous disk state, such as XSS J12270−4859 (see, e.g., de Martino et al. 2010, 2013; Bassa et al. 2014), and candidates such as 3FGL J1544.6−1125 (see, e.g., Bogdanov & Halpern 2015; Britt et al. 2017), CXOU J110926.4−650224 (Coti Zelati et al. 2019, 2021), PSR J0337+1715 (Strader et al. 2016), and XMM J083850.4−282759 (Rea et al. 2017).

5. Conclusions

We presented a detailed timing analysis of (quasi-)simultaneous observations in the X-rays and in the optical band of the transitional millisecond pulsar PSR J1023+0038. The analyzed data cover the time interval from May 2017 to January 2022, when the system was in a subluminous disk state. They were acquired with XMM-Newton and NICER X-ray satellites, and with the fast optical photometers SiFAP2 and Aqueye+. Our main results are summarized below.

1. Although the estimated time lags between optical and X-ray pulsations are not consistent with being modeled by a single value, the weighted averages obtained from the studies of the first and second harmonic terms are ${\overline{\delta \tau }}_1=(175\pm 22)\mathrm{\mu }\mathrm{s}$ and ${\overline{\delta \tau }}_2=(162\pm 6)\mathrm{\mu }\mathrm{s}$, respectively. When we focus on the second harmonic because of its higher power spectral densities, the time lag lies in a limited range of values, ∼(0 − 250) μs. This is maintained over the years, supporting the hypothesis that both pulsations originate from the same region and that their emission mechanisms are intimately linked.

2. From NICER/Aqueye+ simultaneous observations, we found that the second harmonic of the optical pulse lags that of the X-ray pulse by δτ₂ = (254 ± 47) μs in February 2019, (119 ± 48) μs in January 2020, and (108 ± 36) μs in January 2022. These results are not affected by the systematic errors, while the previous measurement by Papitto et al. (2019) was affected by the absolute timing uncertainty of SiFAP2 and XMM-Newton.

3. Our results find a convincing interpretation in the shock-driven mini pulsar nebula scenario (Papitto et al. 2019; Veledina et al. 2019), which suggests an origin of optical and X-ray pulses based on synchrotron radiation emitted from a shock formed where the striped pulsar wind meets the accretion disk, within a few light cylinder radii from the pulsar. The time lag is interpreted in terms of the different timescale that synchrotron X-ray and optical photons take to be emitted.

Within the proposed model, variations of the estimated time lags of optical pulses relative to the X-ray pulses may be due to the variability in disk truncation radius and/or to changes in the shock height. Another reason may be that some data sets lack exactly simultaneous intervals between optical and X-ray observations. There may also be some dependence on source parameters, such as luminosity and orbital phase, that must be investigated in future works.

Acknowledgments

This work is partly based on observations made with the Italian Telescopio Nazionale Galileo (TNG), operated on the island of La Palma by the Fundación Galileo Galilei of the Istituto Nazionale di Astrofisica (INAF), which is installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. Part of this paper is based on observations obtained with the Copernicus Telescope (Asiago, Italy) of the INAF-Osservatorio Astronomico di Padova. Some of the scientific results are obtained from observations with XMM-Newton, which is an European Space Agency (ESA) science mission with instruments and contributions directly funded by ESA member states and NASA, and NICER, which is a NASA mission. We are grateful to the NICER directors for scheduling ToO observations in Director’s Discretionary Time, as well as to all the teams involved for their effort in scheduling the simultaneous ToO observations of PSR J1023+0038. This research also made use of data and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC). G.I. and A.P. wish to thank A. Sanna for useful discussions. G.I. is supported by the AASS Ph.D. joint research program between the University of Rome “Sapienza” and the University of Rome “Tor Vergata”, with the collaboration of the National Institute of Astrophysics (INAF). F. A., G. I., A.P., L.S., D.d.M., and L.Z. acknowledge financial support from the Italian Space Agency (ASI) and National Institute for Astrophysics (INAF) under agreements ASI-INAF I/037/12/0 and ASI-INAF n.2017-14-H.0, from INAF ‘Sostegno alla ricerca scientifica main streams dell’INAF’, Presidential Decree 43/2018. F. A., G. I., A.P., and L.S. also acknowledge funding from the Italian Ministry of University and Research (MUR), PRIN 2020 (prot. 2020BRP57Z). A.P. and D.d.M. acknowledge financial contribution from INAF ‘SKA/CTA projects’, Presidential Decree 70/2016. A.M.Z. is supported by PRIN-MIUR 2017 UnIAM (Unifying Isolated and Accreting Magnetars, PI S. Mereghetti). F.C.Z. is supported by a Juan de la Cierva fellowship. The work of F.C.Z and D.F.T. is also partially supported by the program Unidad de Excelencia María de Maeztu (CEX2020-001058-M). D.F.T. work has been supported by the grants PID2021-124581OB-I00 and acknowledges as well USTC and the Chinese Academy of Sciences Presidential Fellowship Initiative 2021VMA0001. M.T. acknowledges funding from the ERC under the EU’s Horizon 2020 research and innovation programme (consolidator grant agreement no. 101002352).

References

All Tables

All Figures

Fig. 1. Long-term evolution of Tasc as a function of the number of orbital cycles since Tref = 57897.027668 MJD. Red points are the values found by Jaodand et al. (2016), which deviate from the almost linearly increasing trend found in this work. Blue points are from Papitto et al. (2019), and lighter blue points are from Burtovoi et al. (2020). Triangles, associated in the legend with the different NICER observations, are from this work. The thick dashed line indicates our increasing roughly linear trend. The dotted gray line indicates the sinusoidal trend that is similar to the trend found in previous work (e.g., Jaodand et al. 2016; Papitto et al. 2019).

In the text

	Fig. 2. Aqueye+ light curve observed in January 2020, binned every 10 s. The count rate is normalized at the maximum count rate. The vertical colored region indicates the visually identified flaring interval.
In the text

Fig. 3. Time lags relative to the second harmonic term (Table 5). The y-axis is not temporally equispaced. The purple point indicates the values for simultaneous observations between XMM-Newton and SiFAP2, blue points show values between XMM-Newton and Aqueye+, light blue points indicate values between NICER and SiFAP2, and green points show the values between NICER and Aqueye+. The dashed line indicates a zero time lag. Colored error bars represent 1σ statistical uncertainties, while the black error bars indicate the total error. The influence of the absolute timing accuracy of XMM-Newton (∼100 μs) and SiFAP2 (∼60 μs) is visible.

In the text

Fig. 4. Background-subtracted rms amplitudes of X-ray (top panel) and optical (middle panel) pulse during NICER/Aqueye+ simultaneous observations in February 2019, January 2020, and January–February 2022. The arrows represent the upper limits converted into rms amplitudes at 3σ confidence level computed when the optical pulsations were not detected. In the lower panel, we show the phases of the optical pulse (red dots and empty yellow squares show the first and second harmonic, respectively), and the phases of the X-ray pulse (blue dots and light blue empty squares show the first and second harmonic, respectively) during these simultaneous observations.

In the text

Fig. 5. Pulse profiles obtained by folding with 16 phase bins of simultaneous observations with NICER (blue points, left y-axis scale) and Aqueye+ (red points, right y-axis scale) in February 2019 (top panel), January 2020 (middle panel), and January 2022 (bottom panel). Solid lines indicate the best-fitting two-harmonic functions. The zero of the pulse phase was shifted arbitrarily in each panel by the same amount for both instruments.

In the text

Fig. 6. Cross-correlation of the pulse profiles from NICER and Aqueye+ observations of February 2019 as the lag varies. The top and middle panels show NICER and Aqueye+ pulse profiles with 64 phase bins, respectively. Solid lines indicate the best-fitting two-harmonic function. The bottom panel shows the cross-correlation coefficient as a function of the lag with which the optical pulse profile is shifted. We used phase lags ranging from −32 to +32, including the zero lag (dashed vertical line).

In the text