© The Authors 2026

1. Introduction

Millisecond pulsars (MSPs) are old and fast-spinning neutron stars, spun up to periods shorter than 30 ms by the accretion of matter and angular momentum from a companion star (Alpar et al. 1982; Radhakrishnan & Srinivasan 1982; Salvo & Sanna 2021). Most are in low-mass binary systems (< 1 M_⊙), have weak magnetic fields (∼10⁸ − 10⁹ G), and are observed as accreting X-ray pulsars or as rotation-powered radio pulsars. Recent works identified MSPs that can switch between these two regimes, called transitional millisecond pulsars (tMPSs). In particular, these sources experience a transition of state on timescales of a few weeks, passing from a rotation-powered regime, where they act like radio pulsars with winds that prevent accretion onto their surface, to a regime where they accrete matter lost by the companion star, emitting a high amount of X-ray luminosity like X-ray binary systems do (Papitto & De Martino 2021). Only three tMPSs are known so far, PSR J1023+0038 (Archibald et al. 2009), PSR J1227–4853 (De Martino et al. 2010), and PSR J1824-2452 (Papitto et al. 2013).

PSR J1023+0038 was initially classified as a cataclysmic variable (Bond et al. 2002), but between 2000 and 2001 observations showed the existence of an accretion disk, which soon disappeared, leaving a radio pulsar (Thorstensen & Armstrong 2005; Archibald et al. 2009). In 2013, the source underwent a new transition in its emission regime, with the disappearance of the radio pulsations, the reappearance of the accretion disk, and the onset of a strong double-peaked H_α emission observed in the optical spectrum (Halpern et al. 2013; Patruno et al. 2013; Stappers et al. 2014). PSR J1023+0038 is also the only known tMPS that shows detectable UV (Miraval Zanon et al. 2022; Jaodand et al. 2021) and optical pulsations; the latter were first discovered with SiFAP2 at the Telescopio Nazionale Galileo (Ambrosino et al. 2017) and then confirmed with Aqueye+ at the Copernicus telescope in Asiago (Zampieri et al. 2019) and with the panoramic photometer-polarimeter mounted at BTA in Nizhniy Arkhyz (Karpov et al. 2019). The emission mechanism of these optical pulsations is still under debate. Several models have been proposed to explain it (Campana et al. 2019; Papitto et al. 2019; Veledina et al. 2019, and reference therein). In particular, the shock-driven mini-pulsar nebula scenario (Veledina et al. 2019; Papitto et al. 2019) postulates that optical and X-ray pulsations are produced by synchrotron emission caused by shocks that form when the pulsar wind meets the accretion disk. The fact that optical and X-ray pulses show a very short phase lag, between 0 and 0.15, and the fact that both appear only when the source is in the high-mode state, with X-ray pulsations in the high-mode first discovered by Archibald et al. (2015), support this model and a common emission mechanism for the optical and X-ray pulsations (Illiano et al. 2023).

We studied the evolution of the binary system, monitoring the variation in the time of passage at the ascending node, T_asc, and providing a more recent measurement of the phase lag between optical and X-ray pulse profiles. To achieve this, we performed a timing analysis of four Aqueye+ observing runs from 2021 through 2023, as well as of NICER X-ray data from 2023. We derived the timing solution of the source for the 2021 dataset, and we analyzed the evolution of T_asc, incorporating results from previous works (Zampieri et al. 2019; Burtovoi et al. 2020; Illiano et al. 2023; Papitto et al. 2019; Ambrosino et al. 2017; Jaodand et al. 2016; Archibald et al. 2009). Finally, we considered a model for the evolution of the Roche lobe radius of the donor as a function of the mass loss rate, to account for the observed long-term behavior of T_asc.

In Sect. 2 we present the observations and data reduction for the Aqueye+ and NICER observations. In Sect. 3 we describe the correction for the orbital motion and the results of our timing analysis. We discuss our results in Sect. 4, and a summary follows in Sect. 5.

2. Observations and data reduction

We carried out five observing runs of PSR J1023+0038 with the Aqueye+¹ fast photon counter (Barbieri et al. 2009; Naletto et al. 2013; Zampieri et al. 2015) mounted at the 1.8 m Copernicus telescope (Asiago, Italy) in three different years, January 2021, January-February 2022, and January and December 2023 (see Table A.1). The data were reduced with the QUEST software (v. 1.1.5; Zampieri et al. 2015). The arrival time of each photon was referred to the Solar System barycenter using the TEMPO2 package in TDB units (Edwards et al. 2006; Hobbs et al. 2006), using JPL DE405 ephemerides. The position of PSR J1023+0038 was taken from Deller et al. (2012): α = 10h 23m 47.687198s, δ = +00°  38′ 40.84551″ at epoch MJD 54995. We modified the barycentered time series, correcting the photon arrival times for the pulsar orbital motion. The orbital parameters are taken from Archibald et al. (2013): the projected semimajor axis asin(i) = 0.343356 s and the orbital period P_orb = 0.1980963155 d. The reference value for the time of passage at the ascending node is MJD 57449.7258 (Ambrosino et al. 2017). Variations in the time of passage at the ascending node are observed in PSR J1023+0038, the origin of which is not understood but clearly indicates variations in the orbital parameters.

Simultaneous NICER X-ray observations were requested for the Aqueye+ runs of January and December 2023. The corresponding log of the observations is reported in Table A.2. For these observations, we applied the standard data processing procedure², using the script nicerl2 – part of HEASoft³ (v.6.34) software – with version 20240206 of the calibration files. The selected energy range for the analysis is 0.5 − 12 keV. We finally barycentered the data with the script barycorr. Because of the known light-leak issue affecting NICER from May 2023, the observations from December 2023 presented here are only those in which a signal could be detected.

3. Results

We present a timing analysis of the new Aqueye+ optical and NICER X-ray data for PSR J1023+0038, from 2021 through 2023, and analyze the evolution of the passage at the ascending node, T_asc, using observations gathered from 2017 through 2023 and incorporating results from previous works.

3.1. Evolution of the T_asc

The time of passage at the ascending node shows significant variations throughout the epochs. In this work, we report the measure of the T_asc offset for each Aqueye+ run, from January 2021 to December 2023.

For each epoch, we performed an epoch folding search and corrected the barycentered time series for the orbital motion, using the orbital period reported in Jaodand et al. (2016), and varying only the T_asc to find the value that gives the highest signal, i.e., the pulsed profile with the highest χ². We used the T_asc measured by Ambrosino et al. (2017) as the reference value. Then, for each night, we explored a range of 40 values for the T_asc offset (20 in each direction from the expected value), with a time step of 0.5 s. The starting expected value was obtained from the best T_asc offset of the previous run, adding approximately 20 s to it. The final offset was determined by fitting the results of the epoch folding search with a Gaussian. The uncertainty on the best ΔT_asc was calculated as the square root of the sum of the diagonal terms of the covariance matrix of the fit. The final best values of the T_asc are reported in Table 1, along with the final rotational period of the pulsar determined as described in Sect. 3.2. The values of T_asc were calculated with a trial pulsar rotational period P_init = 0.0016879874449 s.

In Fig. 1 we show the T_asc evolution with time using our results and the T_asc offsets obtained from previous works (Jaodand et al. 2016; Ambrosino et al. 2017; Papitto et al. 2019; Burtovoi et al. 2020; Illiano et al. 2023). We then performed a fit with a parabola (Illiano et al. 2023):

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

Fig. 1. Variation in the time of passage at the ascending node, Tasc, concerning the reference value measured in radio reported in Jaodand et al. (2016). The dotted black line represents the fit of the data with Eq. (1).

where T is the time in MJD, and A, B, and C are the free parameters of the fit: A = (2.46 ± 0.17)×10⁴ s, B = −0.87 ± 0.06  s MJD⁻¹, and C = (1.5 ± 0.1)×10⁻⁵  s MJD⁻². The parameter C can be expressed in s/s, C′=(3.52 ± 0.23)×10⁻¹¹. It represents the derivative of the orbital period, and it is higher than the value found in Illiano et al. (2023) by a factor of ∼2. We fitted only the optical data starting from 2017 (Ambrosino et al. 2017), when the variation in T_asc appears to be convincingly consistent with a parabolic trend.

3.2. Folded pulse profiles for the single runs

We report the timing solution for the Aqueye+ data of PSR J1023+0038 for the run of 2021, for which a sufficient number of nights is available to constrain it, and to this purpose, we corrected all the time series for the binary motion using the values of the T_asc (see Table 1) found from the previous analyses, and the orbital period and semimajor axis reported in Jaodand et al. (2016, see our Sect. 2). We first folded the light curves of each night using a guess spin period P_init = 0.0016879874449 s and referring each folded light curve to a reference time t₀ = 59226.0 MJD, i.e., midnight of January 20–21, 2021. The resulting pulsed profiles were fitted with the sum of two harmonically related sinusoids plus a constant (Ambrosino et al. 2017; Zampieri et al. 2019; Burtovoi et al. 2020):

$$\begin{array}{c}\hfill f(x)=K\{1+A_{1}sin(2\pi [x-x_1])+A_{2}sin(4\pi [x-x_2])\},\end{array}$$(2)

where K is a normalization constant, A₁ is the coefficient of the first harmonic, A₂ is the coefficient of the second harmonic, and x₁ and x₂ are the phases of the first and second harmonics, respectively. The measured phases show a drift, Ψ(t), with respect to a uniform rotational spin period, P_init (see Zampieri et al. 2014). To correct the spin period, we fitted this phase drift of the second harmonic with the following first-order polynomial:

$$\begin{array}{c}\hfill \mathrm{\Psi }(t)={\varphi }_0+({\nu }_0-{\nu }_{\mathit{init}})(t-t_0),\end{array}$$(3)

where ϕ₀ and ν₀ are the phase and the spin frequency at the reference time t₀, and ν_init = 1/P_init.

We report in Table 2 the timing solution for the January 2021 run. We show the total pulsed profile of the source in Fig. 2.

Fig. 2. Optical pulse profile of PSR J1023+0038 obtained with Aqueye+ from January 12 to 16, 2021. The red line is the fit with Eq. (2).

The timing analysis for the 2022 Aqueye+ run is reported in Illiano et al. (2023). For the January and December 2023 dataset, the number of nights is not sufficient to perform an analysis similar to that described above; hence, we corrected the data for the binary motion and folded them with the rotational period extrapolated from the accurate timing solution in Burtovoi et al. (2020, values are reported in our Table 1), which differs from the spin period of 2018 by only 4.4 × 10⁻¹³ s. We performed the fit with Eq. (2) and show the final folded profile in Figs. 3 and 4.

Fig. 3. Optical pulse profile of PSR J1023+0038 obtained with Aqueye+ from December 12 and 15, 2023. The red line is the fit with Eq. (2).

Fig. 4. Pulse profiles of PSR J1023+0038 from January 21, 2023, in the optical (Aqueye+, bottom panel, blue dots) and X-ray (NICER, top panel, green dots) energy bands. Both profiles are fitted with the same model (yellow and red lines; see Eq. 2), folded with the same spin period (see Table 1) and with the same t0 = 59965.0 MJD. The dashed navy lines in both panels represent the phase of the second harmonic (x2), used to measure the phase lag, and the light blue shaded bands indicate the corresponding uncertainties.

3.3. Time delay between optical and X-ray data

We measured the delay between the time of arrival of the X-ray and optical pulses using the Aqueye+ and NICER data of January 2023 (Tables A.1 and A.2). We used the simultaneous NICER observation of January 21 to perform an epoch folding search to find the value of the T_asc offset (using only acquisitions when the source was in the high mode for 70 − 80% of the time). We found an offset of 80.25 ± 0.89 s (computed taking the T_asc from Ambrosino et al. 2017 as the reference value). The offset is compatible within the uncertainties with that found from the optical data (Table 1). We corrected both the X-ray and optical datasets for binary motion (considering only the optical observations within the X-ray high-mode window) with the T_asc offset of the X-ray data (80.25 s), the signals of which have a higher root mean square pulse amplitude. We also performed an epoch folding search using all the NICER observations reported in Table A.2 to find the T_asc offset for January and December 2023; we obtained 80.43 ± 0.92  s and 96.86 ± 1.13  s, respectively, which are consistent with the optical data (see Fig. 1).

Then, we folded the light curves of the simultaneous observations using the same t₀ (59965.0 MJD) and the spin period extrapolated from the Burtovoi et al. (2020) long-term timing solution reported in Table 1. We then fitted both profiles with Eq. (2) and measured the related drift of the second harmonic, whose power spectral density is higher than that of the first harmonic. The optical pulses are found to be delayed with respect to the X-ray ones by 112.3 ± 30.7 μs, which corresponds to a phase lag of 0.067 ± 0.018. In Fig. 4 we show the X-ray and optical pulse profiles for January 21, 2023.

4. Discussion

4.1. Evolution of T_asc

In Sect. 3.1 we present the variation in the time of passage at the ascending node from 2021 through 2023. We report our results in Fig. 1, alongside those obtained in previous works. We fitted the ΔT_asc evolution with a second-order polynomial (see Eq. (1)), which is consistent with the observed trend of the data, in particular if we consider the ΔT_asc values after the observations performed by Ambrosino et al. (2017). In the past, the variation in the T_asc was consistent with an irregular evolutionary scenario for the binary system; this scenario is further supported by its initial oscillating behavior after the onset of the accreting phase. Recent results, however, seem to support a different scenario, in which the increase in the T_asc could be driven by a rise in the orbital period and, thus, an increase in the separation between the companion star and the neutron star. In this section we discuss a possible general evolutionary scenario consistent with our recent results and the characteristics of this binary system.

Assuming accretion through Roche lobe overflow (Verbunt 1993; Frank et al. 2002), the exchange of mass leads to a variation in the distance between the two stars, due to the conservation of angular momentum, which translates into variations in the orbital period. Let us consider the general equation of the total angular momentum of the system, J, and its time derivative, $\dot{J}$, assuming a circular orbit:

$$\begin{array}{c}\hfill J=M_1{M}_2\sqrt{\frac{\mathit{Ga}}{M_1+M_2}}{M}_{\odot }^{3/2},\end{array}$$(4)

where G is the gravitational constant, a is the orbital separation, and M₁ and M₂ are the masses of the neutron star and the donor star, respectively. The time derivative of Eq. (4) is

$$\begin{array}{c}\hfill 2\frac{\dot{J}}{J}=\frac{\dot{a}}{a}+2(\frac{\dot{M_1}}{M_1}+\frac{\dot{M_2}}{M_2})-\frac{\dot{M_1}+\dot{M_2}}{M_1+M_2},\end{array}$$(5)

and the total mass and mass ratio of the system are defined as

$$\begin{array}{c}\hfill M=M_1+M_2,q=\frac{M_2}{M_1}.\end{array}$$(6)

We used the equation that approximates the evolution of the Roche lobe radius of the donor star, R_2L, for 0.1 ≤ q ≤ 0.8, which is given by (Frank et al. 2002)

$$\begin{array}{c}\hfill \frac{\dot{R_{2L}}}{R_{2L}}=\frac{\dot{a}}{a}+\frac{\dot{M_2}}{3M_2}-\frac{1}{3}\frac{\dot{M}}{M}.\end{array}$$(7)

The accretion luminosity in X-rays is L_X = 7 × 10³³ erg s⁻¹ in the high-mode state, and L_X = 3 × 10³³ erg s⁻¹ in the low-mode state (Messa et al. 2024). Neglecting emission from the neutron star surface, we estimated the mass accretion rate of the source as

$$\begin{array}{c}\hfill L_{\mathrm{acc}}=\xi \frac{GM{\dot{M}}_{\mathrm{acc}}}{R_{★}},\end{array}$$(8)

where G is the gravitational constant, M the mass of the neutron star, Ṁ_acc the mass accretion rate, R_★ the radius of the neutron star, and ξ the accretion efficiency, which depends on the inner radius, R_in, of the accretion disk (ξ = R_★/R_in). In the case of PSR J1023+0038, in the high-mode state, material channeled along the magnetic field lines is likely not accreted onto the neutron star (as it would otherwise quench the pulsar wind) but rather expelled in a jet (Baglio et al. 2023). In the Papitto et al. (2019) model, the inner radius of the disk, in the high-mode state, is beyond the light cylinder radius, at ≃80; however, when the source is in the low-mode state, the disk recedes up to a distance of ≃1600 km, and most of the X-ray emission is produced by the compact jet along the pulsar rotational axis. Assuming R ≃ 80 km, and considering the X-ray accretion luminosity in the high mode, we obtain

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{acc}}=4.22\times {10}^{14}{\left(\frac{M}{M_{\odot }}\right)}^{-1}\mathrm{g}·{\mathrm{s}}^{-1}.\end{array}$$(9)

For a conservative mass transfer, the value of the mass accretion rate is between Ṁ₂ = −1.1 ×10¹⁹ g s⁻¹ (if q = 0.1) and Ṁ₂ = −2.5 ×10¹⁸ g s⁻¹ (if q = 0.8; see Frank et al. 2002). These values, being at least four orders of magnitude higher than the value obtained from the observed X-ray luminosity, demonstrate the inconsistency of a conservative accretion scenario.

The fact that the pulsar in PSR J1023+0038 is active implies that a powerful pulsar wind hits the donor surface, contributing to slowly evaporate the companion, as observed in other “redback binary systems”. We can then consider a nonconservative case in which the accretion rate onto the neutron star is negligible and the mass loss is dominated by the donor, Ṁ₁ ≃ Ṁ_acc ≪ −Ṁ₂, and Ṁ ≃ Ṁ₂. Then, Eq. (5) assumes the following form:

$$\begin{array}{c}\hfill 2\frac{\dot{J}}{J}=\frac{\dot{a}}{a}+2\left(\frac{\dot{M_2}}{M_1}\right)-\frac{\dot{M_2}}{M_1+M_2}.\end{array}$$(10)

Taking the time derivative of Kepler’s third law, under the same assumptions, we obtain

$$\begin{array}{c}\hfill \frac{\dot{a}}{a}=\frac{1}{3}\frac{\dot{M_2}}{M_2}\frac{M_2}{M}+\frac{2}{3}\frac{\dot{P}}{P},\end{array}$$(11)

where P is the orbital period. Since we do not know exactly how the companion star evolves as it interacts with the pulsar wind, we considered a general case in which the radius of the donor evolves as

$$\begin{array}{c}\hfill \frac{\dot{R_2}}{R_2}=\alpha \frac{\dot{M_2}}{M_2}.\end{array}$$(12)

When α = 1, this relation represents the evolution of a low-mass donor that is losing mass in a quasi-steady state. Assuming that the system is in (or close to) Roche lobe overflow, and hence R₂ ≃ R_2L, we substituted this expression into Eq. (7) and then combined it with Eq. (11). The resulting relation is

$$\begin{array}{c}\hfill \frac{\dot{P}}{P}=\frac{3}{2}\beta \frac{\dot{M_2}}{M_2},\mathrm{where}\beta =(\alpha -\frac{1}{3}).\end{array}$$(13)

As mentioned above, the case α = 1 represents the evolution of a low-mass donor in a quasi-steady state. This is not a plausible scenario since β = −2/3 and Ṗ/P = Ṁ₂/M₂ < 0, which is at variance with the observed increase in the orbital period. The same conclusion holds in general for α > 1/3, while for α < 1/3, Ṗ/P ∝ −Ṁ₂/M₂ > 0, which is consistent with the observations. We note that the case α = 0 (β = 1/3) would represent the evolution of a system with a constant donor and Roche lobe radius. In these conditions, the mass-loss rate of the donor would be −Ṁ₂/M₂ = (1/2)(Ṗ/P) = 8.2 ×10¹⁸ gs⁻¹, higher by four orders of magnitude than the mass accretion rate. As we do not know the physics of the interaction between the pulsar wind and the donor, different values of α would be possible – implying different evolutions for the donor and different values for the mass-loss rate – but they must be on the same order of magnitude as the previous estimate for |α|≈1. At this pace, the donor would be completely consumed in ∼8 M yr. Thus, a scenario with a binary system evolving through nonconservative Roche lobe overflow, similar to the one presented above, appears to be consistent with the observed steady parabolic increase in the orbital period with time.

4.2. Time delay between optical and X-rays

We analyzed NICER and Aqueye+ data from January 2023 to perform new measurements of the time shift between the optical and X-ray pulses of PSR J1023+0034. We managed to measure the shift simultaneously between optical and X-rays only for the January 2023 data (January 21) because for the December 2023 dataset, the duration of the simultaneous X-ray observation (December 15) was too short. The time lag between optical and X-ray pulses is 0.036 ± 0.016, which corresponds to 60.7 ± 27.7 μs. This result is consistent with that of Illiano et al. (2023), who report that, from the second harmonic of the pulsed profiles, the time lag between optical and X-ray pulses lies in the range (0 − 250) μs. This, along with the fact that the optical and X-ray pulses have similar morphologies, confirms that they are produced in the same region and with the same emission mechanism, as proposed by Papitto et al. (2019) and Veledina et al. (2019, see the references therein). They present a possible scenario to explain such a simultaneous emission, the so-called striped wind model. In this scenario, the source is in the disk state but still able to emit a relativistic wind. In a region a few light cylinder radii from the star, the pulsar wind interacts with the accretion disk, producing two shocked symmetric regions; this is where synchrotron radiation is produced and pulsed emission in the X-ray and optical bands originates. They also estimated the synchrotron timescale for optical photons, obtaining a value compatible with the lag observed between optical and X-ray photons.

5. Summary

We have studied the evolution of the PSR J1023+0038 binary system, monitoring the long-term variation in the time of passage at the ascending node, T_asc, and providing new measurements of the phase lag between the optical and X-ray pulse profiles. We performed a timing analysis of three runs of Aqueye+ optical data from 2021 through 2023, as well as NICER X-ray data from 2023. We derived the timing solution of the source for January 2021, and we analyzed the evolution of the T_asc, including also results from previous works. We measured the shift between optical and X-ray pulses using simultaneous Aqueye+ optical data and NICER X-ray data from January 2023. The phase lag measured is 0.036 ± 0.016, in agreement with Illiano et al. (2023) and extending the baseline of these measurements by one year. This result further confirms the common origin and emission mechanism of X-ray and optical pulses. Concerning the evolution of the T_asc, starting from 2017, it follows a parabolic increasing trend, which translates into an increase in the orbital period and orbital separation of the system. A scenario of a binary system evolving through nonconservative Roche lobe overflow in which the donor is interacting with the pulsar wind, losing mass at a rate much higher than the inferred accretion rate, appears to be consistent with the observed steady increase in the orbital period of the system. The large majority of the mass stripped from the donor was likely shocked by the pulsar wind itself and ejected from the system after crossing the Roche lobe. Further optical and X-ray timing observations of PSR J1023+0038 are needed to continue to monitor the evolution of the parameters of the binary system and to determine the long-term rotational history of the pulsar.

Acknowledgments

SC is supported by the INAF Doctorate program. SC and LZ acknowledge financial support from the INAF Research Grant “Uncovering the optical beat of the fastest magnetized neutron stars (FANS)”. Based on observations collected at the Copernico 1.82m telescope (Asiago Mount Ekar, Italy), INAF – Osservatorio Astronomico di Padova. This work has made use of data products and software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. We would like to express our thanks to Alessandro Papitto, Giulia Illiano, and Arianna Miraval Zanon for the very interesting discussions on the content of this article, in particular Alessandro Papitto for his support in training S.C. on pulsar timing techniques.

References

Appendix A: Observations

All Tables

All Figures

	Fig. 1. Variation in the time of passage at the ascending node, Tasc, concerning the reference value measured in radio reported in Jaodand et al. (2016). The dotted black line represents the fit of the data with Eq. (1).
In the text

	Fig. 2. Optical pulse profile of PSR J1023+0038 obtained with Aqueye+ from January 12 to 16, 2021. The red line is the fit with Eq. (2).
In the text

	Fig. 3. Optical pulse profile of PSR J1023+0038 obtained with Aqueye+ from December 12 and 15, 2023. The red line is the fit with Eq. (2).
In the text

Fig. 4. Pulse profiles of PSR J1023+0038 from January 21, 2023, in the optical (Aqueye+, bottom panel, blue dots) and X-ray (NICER, top panel, green dots) energy bands. Both profiles are fitted with the same model (yellow and red lines; see Eq. 2), folded with the same spin period (see Table 1) and with the same t0 = 59965.0 MJD. The dashed navy lines in both panels represent the phase of the second harmonic (x2), used to measure the phase lag, and the light blue shaded bands indicate the corresponding uncertainties.

In the text