© The Authors 2026

1. Introduction

Millisecond pulsars (MSPs) constitute a subclass of rapidly rotating neutron stars (NSs) that demonstrate coherent millisecond pulsations in one or more energy ranges from the radio to gamma-rays. The first representative of this class, radio pulsar PSR B1937+21, was discovered in 1982 (Backer et al. 1982) and its nature was immediately explained on the basis of the recycling scenario (Alpar et al. 1982; Radhakrishnan & Srinivasan 1982). According to this view, first introduced theoretically by Bisnovatyi-Kogan & Komberg (1974), the NS is spun-up by accretion of matter and angular momentum from a low-mass companion star in a close binary system. During the accretion stage, the NS is expected to exhibit X-ray pulsations due to magnetic channeling of the accreting matter and release of its kinetic energy in the magnetic pole regions. Thus, the discovery of a recycled accreting X-ray MSP, SAX J1808.4-3658 (Wijnands & van der Klis 1998), in 1998 was a strong confirmation of the presented picture (see also detection of millisecond pulsations from 4U 1728-34 during X-ray bursts by Strohmayer et al. 1996).

However, the most compelling evidence came with the discovery of transitional MSPs. The first one, PSR J1023+0038 (hereafter J1023), was initially classified in early 2000s as a cataclysmic variable due to its radio emission, optical variability and prominent double-peaked emission lines indicating presence of an accretion disk (Bond et al. 2002; Szkody et al. 2003). Since early 2002, the source has shown no further sign of accretion or emission lines (e.g., Thorstensen & Armstrong 2005; Homer et al. 2006). Moreover, in June 2007 a bright MSP was found in the radio (Archibald et al. 2009), undoubtedly indicating that the system performed the transition from the low-mass X-ray binary (LMXB) to the radio pulsar state. For six years, J1023 behaved like a typical redback pulsar (Archibald et al. 2010; Bogdanov et al. 2011; Archibald et al. 2010), but in June 2013 it transitioned back to the subluminous disk state (Patruno et al. 2014; Stappers et al. 2014) where it has remained since then. At the moment, only two other transitional MSPs have been conclusively identified, IGR J18245-2452 (Papitto et al. 2013) and XSS J12270-4859 (Bassa et al. 2014; Roy et al. 2015). A few candidates with a similar morphology await confirmation (Bogdanov & Halpern 2015; Strader et al. 2016; Li et al. 2020; Coti Zelati et al. 2019; Miller et al. 2020; Strader et al. 2021; Kyer et al. 2025; Zyuzin et al. 2025).

J1023 is characterized by extremely rich variety of observational properties revealed through a numerous multi-wavelength campaigns in recent years. We refer to Table 2 from Jaodand et al. (2016) and Table 6.1 from Papitto & de Martino (2022) for a summary of its main parameters. In the LMXB state, J1023 exhibits sporadic switches between three different X-ray modes (Bogdanov et al. 2015). About 70 percent of the time, it remains in the high mode with the X-ray luminosity in the 0.3–10 keV band L_X ≃ 3 × 10³³ erg s⁻¹, with rapid switches to the low mode (L_X ≃ 5 × 10³² erg s⁻¹) and occasional flares (L_X ≃ 10³⁴ erg s⁻¹), also observed in the optical band (Papitto et al. 2019). These luminosities are considerably lower than the values typical for accreting MSPs during outbursts (∼10³⁶ − 10³⁷ erg s⁻¹, Campana & Di Salvo 2018). Moreover, for the NS mass of 1.7 M_⊙ (Deller et al. 2012) and typical radius of 10 km, the corresponding accretion rates $\dot{\mathfrak{M}}\simeq 2\times {10}^{12}-4\times {10}^{13}\text{g}{\text{s}}^{-1}$ are insufficient to overcome the centrifugal barrier imposed by a rapidly rotating magnetosphere of the NS (Archibald et al. 2015; Papitto et al. 2019). Coherent X-ray pulsations are detected only during the high mode (Archibald et al. 2015; Papitto et al. 2019), although Bogdanov et al. (2015) reported a marginal (with a single-trial significance of ∼2σ) signal in the flare mode.

The isotropic optical luminosity in the 320–900 nm band is L_opt ∼ 10³³ erg s⁻¹ (Papitto et al. 2019). In the optical and infrared, J1023 displays sporadic activity on seconds to hours time scales (Shahbaz et al. 2015, 2018; Hakala & Kajava 2018; Kennedy et al. 2018; Messa et al. 2024). The minimal stochastic variability time scale of the optical emission is as short as fractions of seconds (Shibanov et al. 2017; Beskin et al. 2018), which is close to the characteristic time scales of variability due to matter fragmentation in a propeller regime demonstrated by MHD simulations (Romanova et al. 2018). The presence of some kind of an outflow is suggested by highly variable, flat-spectrum continuum radio emission anti-correlated with X-rays (Deller et al. 2015; Bogdanov et al. 2015, 2018).

Surprisingly, coherent optical (Ambrosino et al. 2017; Zampieri et al. 2019; Karpov et al. 2019) and near-UV (Jaodand et al. 2021; Miraval Zanon et al. 2022) pulsations at the NS spin frequency were detected. Simultaneous X-ray and optical observations (Papitto et al. 2019) indicated that optical pulsations are present during the high X-ray mode, having a double-peaked profile similar to that seen in X-rays. Analogous behavior is observed in the UV (Miraval Zanon et al. 2022). Pulsations in all three energy ranges are likely produced in regions located within a few kilometers from each other and by the same mechanism. Papitto et al. (2019) also detected optical pulsations with reduced amplitude in the flare mode.

A number of models have been proposed to explain such a wealth of observational phenomena seen in J1023. The main questions include whether the ejector mechanism continues to operate in the presence of the accretion disk, how strong the propeller mechanism is, whether a significant fraction of matter ends up reaching the NS surface, what the origins of the X-ray and optical emission are, and what process governs the switches between high, low, and flaring X-ray modes (see e.g. Coti Zelati et al. 2014; Li et al. 2014; Stappers et al. 2014; Takata et al. 2014; Papitto & Torres 2015; Campana et al. 2016). The discovery of the optical pulsations essentially ruled out accretion onto the NS surface as the mechanism behind the generation of the pulsing emission in the optical and X-rays. Papitto et al. (2019) and Veledina et al. (2019) suggested instead that such emission originates in the shock powered by collision between the pulsar wind and the disk just beyond the light cylinder during the high mode. A number of observational works generally support this scenario (Campana et al. 2019; Burtovoi et al. 2020; Illiano et al. 2023; Baglio et al. 2023, 2025a,b, also see however Jadoliya et al. 2026). We discuss this in more detail in Sect. 4.1.

In 2017–2020, we carried out an observational campaign with the 6-m BTA telescope, focusing on the coherent optical millisecond pulsations of J1023. The system was observed over seven nights, providing about 10 hours of total exposure time. Here, we present an analysis of the data obtained over two consecutive nights on November 14–15, 2017, featuring a puzzling behavior of the optical pulse profile never seen before in pulsars with smooth, isotropic-like light curves. The rest of the data, along with the general properties of observed millisecond pulsations and subsecond flaring activity will be discussed in subsequent papers.

Observations and data reduction are described in Sect. 2, while the results are presented in Sect. 3. In Sect. 4, we discuss possible scenarios capable of explaining the observed behavior of the system. We draw our conclusions in Sect. 5. In Appendix A, we briefly consider Swift/XRT X-ray observations of J1023 two weeks after our data set.

2. Observations and data reduction

We observed J1023 with the Multi-channel panoramic photo-polarimeter (MPPP, Plokhotnichenko et al. 2021) mounted on the 6-m BTA telescope of the Special Astrophysical Observatory of the Russian Academy of Sciences (SAO RAS). The observation log is presented in Table 1. The device operated in dual-channel regime allowing to collect photons simultaneously with the two MCP-based panoramic photon counters, the “red” one with the GaAs photocathode sensitive at the 5670 Å effective wavelength, and the “blue” one with the multi-alkali photocathode sensitive at the 4320 Å effective wavelength. Photons in a 10″ × 10″ diaphragm around the object were detected and registered with the effective temporal resolution of ∼1 μs.

For to technical reasons, no comparison star was observed, precluding the standard flux correction procedure. To roughly estimate the total background contribution, which is essential for calculation of the pulsed fraction of the source, we considered 10 second-long frames, taking the mean flux in the source-free area of each frame as the sky contribution to the source flux, similarly to what is usually done in X-ray astronomy. We note that this approach is possible since we are only interested in the total background contribution to the pulse profile, rather than its short-timescale fluctuations. As a result, we estimated the number of the background photons as 3461883 ± 1860 (0.47 of the total) and 945797 ± 973 (0.34 of the total) in the red and blue channels for the first night, and 6074874 ± 2465 (0.38 of the total) and 1398271 ± 1182 (0.24 of the total) in the red and blue channels for the second night.

We followed the standard procedure to search for periodic signal in our data. The photon times of arrival (ToAs) were converted to the barycenter of the Solar system using DE405 ephemeris. To account for orbital motion of the pulsar in the binary system, the ToAs are usually converted to the epoch of passage at the ascending node, T_asc,

$$\begin{array}{c}\hfill t=t_0-a_{0}sin2\pi \phi ,\end{array}$$(1)

where t₀ and t are the ToAs before and after correction, φ = 2π(t − T_asc)/P_orb is the orbital phase of the NS¹, a₀ is the projection of the semi-major axis of the orbit to the line of sight², and P_orb is the orbital period.

The orbital phase of J1023 (or, equivalently, the time of passage at the ascending node) is known to exhibit rather peculiar behavior. After the transition in 2014 from the radio pulsar state to the subluminous disk state, the time of passage at the ascending node shifted by about 25 seconds (back in time) and then for a few years varied around this value (e.g. Jaodand et al. 2016; Zampieri et al. 2019; Papitto et al. 2019). Sometime in 2019 these variations turned into an increasing trend, as was shown by Illiano et al. (2023) using NICER observations and recently confirmed by Conforti et al. (2026) on a larger data set (see their Fig. 1), which could be associated with underestimation of the orbital period. In any case, fine tuning of the orbital phase is required. We calculated the ephemeris orbital phase based on the solution by Jaodand et al. (2016) and considered a set of corrections in the interval from −0.0025 to 0.0025 with the resolution of 10⁻⁵. This interval covers all stochastic variations of the orbital phase ever observed in J1023. For each correction we calculated the Z_K²-statistic (Buccheri et al. 1983) as

$$\begin{array}{c}\hfill Z_K^2=\frac{2}{N}{\displaystyle \sum _{k=1}^K}[{({\displaystyle \sum _{i=1}^N}cos2\pi k{\varphi }_i)}^2+{({\displaystyle \sum _{i=1}^N}sin2\pi k{\varphi }_i)}^2],\end{array}$$(2)

taking into account the first two harmonics (K = 2) of the pulsar frequency since the optical pulse profile is known to have a smooth double-peaked structure. Here, N is the total number of collected photons and ϕ_i is the rotation phase of the NS corresponding to the ToA of the i^th photon, which depends on the assumed pulsar frequency and orbital phase (through the ToAs correction given by Eq. (1)). The orbital period of the system and the projection of the semi-major axis of the orbit to the line of sight were fixed at the values reported in Jaodand et al. (2016). The pulsar frequency was fixed at f = 592.421467517559 Hz, which is the Jaodand et al. (2016) frequency corrected for the epoch difference using the frequency derivative $\dot{f}=-3.0413\times {10}^{-15}$ Hz s⁻¹. The result is shown in Fig. 1. A strong periodic signal (Z₂² ≃ 250) is found at the orbital phase correction of 0.00129(9), translating to −22 s shift of the time of passage at the ascending node with respect to 54905.97140075 MJD (Jaodand et al. 2016). The associated uncertainty is the half width at half maximum (HWHM) of a Gaussian fitted to the peak. The corrected value at our epoch (58072 MJD) is T_asc = 54905.97115(2) MJD. The observations were taken in 2017, so this value is compatible with the behavior observed back then. The significance of the peak is very high. For example, for a 5σ confidence level with 500 independent trials, the detection threshold for Z₂² (which for the Poisson noise is distributed as χ² with 4 degrees of freedom) is 48. We also checked that the ephemeris pulsar frequency is indeed the optimal frequency in terms of the signal strength and no correction is needed in this regard.

Fig. 1. Z22-test for both nights joined. The orbital phase correction is shown with respect to the ephemeris value obtained from the timing solution by Jaodand et al. (2016). The maximum at 0.00129 corresponds to the shift of the epoch of passage at the ascending node by −22 seconds with respect to 54905.97140075 MJD. The aliases around the main peak emerge due to the presence of gaps in each data set and between them.

After the specified reduction and optimization of orbital parameters, we constructed binned pulse profiles of the source at the rotation period. The numbers of phase bins were estimated using Eq. (1) from Becker & Trümper (1999), which compromises between the pulse profile resolution and signal to noise ratio.

Preliminary search for periodic signal in the data with a 100-second-wide sliding window performed by Karpov et al. (2019) revealed one short interval of peculiar activity with strong, single-peaked pulsations. Having that in mind, in addition to folding the entire data set, we also analyzed individual chunks of the data on time scales from tens to hundreds of seconds.

To characterize the strength of the pulsing component, we calculated the background-corrected pulsed fractions in the obtained pulse profiles. Different quantities are used in the literature as a representation of the pulsed fraction (see e.g. Appendix C of Hare et al. 2021, for discussion). We define p_area as the ratio of the area between the pulse profile and the lowest phase bin level to the full area under the curve,

$$\begin{array}{c}\hfill p_{\mathrm{area}}=\frac{C-c}{C-b},\end{array}$$(3)

where C is the total amount of counts, c is the unpulsed level equal to the number of counts in the lowest phase bin times the number of bins, and b is the estimated total amount of background counts. For comparison with previous studies, we also provide the values of the root mean square (rms) pulsed fraction, p_rms,

$$\begin{array}{c}\hfill p_{\mathrm{rms}}={\left(\frac{Z_2^2}{N}\right)}^{1/2},\end{array}$$(4)

which essentially characterizes the rms amplitude of the harmonic representation of the pulse profile (see also e.g. Papitto et al. 2019, for a similar definition). We note that p_rms is systematically lower than p_area (by a factor of $\sqrt{2}$ for a sinusoidal signal; Hare et al. 2021).

3. Results

For most of our data, we observed a typical double-peaked pulse profile very similar to those reported in a number of previous studies (Ambrosino et al. 2017; Zampieri et al. 2019; Papitto et al. 2019; Burtovoi et al. 2020; Illiano et al. 2023). The averaged background-corrected pulsed fractions for the red and blue channels combined (all the detected photons considered together) are p_area = 0.65 ± 0.04 percent and p_rms = 0.43 ± 0.03 percent. The two peaks are approximately half of a phase apart, with the first one about two times higher than the second, although the heights of the peaks show strong variability on shorter timescales. In contrast to previous studies, the system was observed simultaneously in two optical bands. The “normal” pulse profiles in the red and blue channels are shown on the top panels of Fig. 2.

Fig. 2. Top: Typical double-peaked pulse profiles obtained by phase folding together the Nov 14 and Nov 15 light curves and excluding the single-peaked interval. The red and blue channels are shown on the left and right, respectively. Bottom: Averaged single-peaked pulse profiles of the 220-second-long interval. Corresponding pulsed fractions parea and prms are shown in each panel.

However, toward the end of the second night an outstanding 220-second-long event took place, starting at 58072.11364 MJD. In the bottom panels of Fig. 2, the pulse profiles of this interval are shown. The first thing that catches an eye is a completely different pulse shape with a broad, smooth, slightly skewed single peak. It is important to note that the maximum of this peak does not match either of the two in the normal state. The averaged background-corrected pulsed fractions for the sum of the red and blue channels are p_area = 3.3 ± 0.2 percent and p_rms = 1.9 ± 0.1 percent (i.e., about four to five times higher than that of the double-peaked pulse profile). Inspecting the light curve of the source, we found that the single-peaked interval is encompassed by flaring activity (Fig. 3). While flares are consistently observed in J1023 (e.g. Shibanov et al. 2017; Beskin et al. 2018; Shahbaz et al. 2018; Papitto et al. 2019), we note that neither the two nights under consideration nor any other data we obtained in 2017–2020 contain similar single-peaked and/or high-amplitude pulsation intervals.

Fig. 3. Segment of the light curve in the red (top) and blue (bottom) channels containing the interval with single-peaked pulsations (shown with a hatched stripe). Flaring activity around the interval is readily seen. The time resolution is 1 s. Zero time corresponds to 15 Nov 2017 02:37:29 UT.

Thanks to much stronger signal, we are able to investigate temporal variability of the pulse profile during this event at timescales as short as about 10 seconds. Figure 4 provides an example of such variability. Prominent pulsations seen in the first 10 seconds of the event (panel a of Fig. 4, p_area = 5.2 ± 0.9 percent, with a detection significance of 3.8σ) almost cease in noise in the next 10 seconds (panel b, p_area = 2.4 ± 1.1 percent, 1.7σ) and rise up again to p_area = 5.3 ± 0.9 percent (panels c, 2.8σ) and p_area = 4.7 ± 0.9 percent (panel d, 3.6σ). Similar behavior is observed throughout the whole interval. Toward the end of the event, the pulse profile first brings back its double-peaked shape (Fig. 5), and only after that, with a few tens of seconds delay, the pulsed fraction drops to its normal level of ≲1 percent, adding one more piece to the puzzle.

Fig. 4. Pulse profile variability during the first four 10-second-long segments (a, b, c, and d) of the single-peaked interval. Profiles observed within the time ranges indicated in the legends are plotted in black. Time is counted from the onset of the single peak. Pulsed fractions and detection significance levels are shown in each panel. Blue and orange filled step plots show the averaged pulse profile of the full single-peaked interval and the averaged over both nights typical double-peaked pulse profile, respectively.

Fig. 5. Same as Fig. 4, but for the first 40 seconds after the single-peaked interval. The pulsed fraction remains high, but the peak is shifted in phase and the pulse profile starts to resemble the normal double-peaked one, although the lower peak is lost in noise at such a short integration time.

The onset of the event takes place very quickly and the pulsed fraction is the highest during the first seconds of the interval. As we can see from Fig. 6, the single-peaked profile becomes prominent at the integration time as short as 5 seconds from the onset of the event and in this interval the pulsed fraction p_area exceeds 5 percent.

Fig. 6. Top to bottom: Pulse profiles (red and blue channels combined) at different integration times (indicated in the legends) from the onset of the single-peaked interval. Corresponding background-corrected pulsed fractions and significance levels are given in each panel.

Simultaneous observations in the red and blue channels allow us to roughly analyze the behavior of the pulsar in two optical bands. The averaged background-corrected pulsed fractions of the double-peaked pulse profile are p_area = 0.69 ± 0.05 percent and p_rms = 0.44 ± 0.04 percent in the red channel, and p_area = 0.60 ± 0.07 percent and p_rms = 0.39 ± 0.05 percent in the blue channel. For the single-peaked interval, the averaged values are p_area = 3.1 ± 0.3 percent and p_rms = 2.1 ± 0.2 percent in the red channel, and p_area = 3.1 ± 0.3 percent and p_rms = 1.6 ± 0.2 percent in the blue channel, respectively.

In the absence of simultaneous spectroscopic observations, we are not able to reliably translate the observed pulsed fractions to the fluxes of the pulsed component. However, we can roughly analyze the relative evolution of the energy distribution in the pulsed emission during the observed transformation to the single-peaked state (as we show below), assessing the long-term (years) and short-term (minutes) variability of the spectrum shape. There are a number of spectra published in the literature taken after the transition of the system to the subluminous disk state. Here, we considered three of them, which we refer to as spectrum 1 (Bogdanov et al. 2015), spectrum 2 (Hernandez Santisteban 2016), and spectrum 3 (Messa et al. 2024, with reduced spectra available³ online). We assumed a simple power-law energy distribution of the pulsed component, F_ν ∼ ν^−α, and first examined whether its slope varies significantly from one spectrum to another. For this purpose, for each spectrum, we calculate the flux densities at the effective frequencies of the red (ν_eff = 5.36 × 10¹⁴ Hz) and blue (ν_eff = 7.01 × 10¹⁴ Hz) channels. The resulting values are F_eff ≃ 5.1 × 10⁻²⁷ erg cm⁻² s⁻¹ Hz⁻¹ (red) and F_eff ≃ 3.8 × 10⁻²⁷ erg cm⁻² s⁻¹ Hz⁻¹ (blue) for spectrum 1, F_eff ≃ 8.0 × 10⁻²⁷ erg cm⁻² s⁻¹ Hz⁻¹ (red) and F_eff ≃ 6.1 × 10⁻²⁷ erg cm⁻² s⁻¹ Hz⁻¹ (blue) for spectrum 2, and F_eff ≃ 5.0 × 10⁻²⁷ erg cm⁻² s⁻¹ Hz⁻¹ (red) and F_eff ≃ 4.0 × 10⁻²⁷ erg cm⁻² s⁻¹ Hz⁻¹ (blue) for spectrum 3. Then we multiply these values by the corresponding pulsed fractions, p_area, measured in the double-peaked and single-peaked states to obtain the flux densities of the pulsed emission, and calculate the corresponding spectral indices α. For the double-peaked state, we find α = 1.64 ± 0.52 for spectrum 1, α = 1.52 ± 0.52 for spectrum 2, and α = 1.38 ± 0.52 for spectrum 3. These values are in a fairly good agreement with each other. Also, inspection of spectra obtained by Messa et al. (2024) in high-temporal-resolution spectroscopic observations shows that the shape of the continuum does not change significantly on timescales of a few minutes. For the single-peaked state, which lasted only ∼3.5 minutes, α = 1.12 ± 0.50 for spectrum 1, α = 1.00 ± 0.49 for spectrum 2, and α = 0.86 ± 0.50 for spectrum 3. So, the results of our analysis of the relative evolution of the spectrum during the transformation to the single-peaked state hint at hardening of the energy distribution in the pulsed component, although the uncertainties of the estimated values of α are large. Simultaneous spectral and timing observations of similar events are needed in the future to confirm this result.

4. Discussion

In this paper, we report the discovery of a transformation of the optical pulse profile of J1023 from its normal, double-peaked shape to the single-peaked and back. This 220-second episode was accompanied by flaring activity observed in the light curve. In the absence of simultaneous X-ray observations, we were not able to identify the X-ray mode of the system during the event. At the same time, coordinated X-ray and optical observations show that prominent optical flares usually occur during X-ray flaring mode (Papitto et al. 2019). Interestingly, in this paper, the authors report the detection of double-peaked pulsations during a strong 3-ks flare with the rms amplitude of 0.16 ± 0.02 percent, which is more than five times lower than during the high mode. They note that such a low value cannot be explained only by addition of an unpulsed flux to the high mode level, which means that the pulsed flux was lower during the flare. In contrast, in our event we observed a ∼5 times increase of the pulsed fraction, which hints at a different nature of this particular flare. Other flares observed in our data are much shorter than the one reported in Papitto et al. (2019). They probably also contain pulsations with reduced amplitudes, but we were not able to distinguish them due to insufficient count statistics.

4.1. Double-peaked pulsations

It was suggested by Papitto et al. (2019) and Veledina et al. (2019) that generation of the non-thermal optical and X-ray emission during the high mode is powered by collision of the pulsar wind with the accretion disk beyond the light cylinder resulting in formation of a shock, acceleration of electrons, and, consecutively, synchrotron emission in a wide range of frequencies due to motion of accelerated electrons in magnetic fields. Veledina et al. (2019) further linked the flaring mode to rapid inflation of the inner parts of the accretion disk, which naturally explains the reduction of the pulsed fraction. This scenario appears promising, and important constraints were obtained in a number of observational works (Campana et al. 2019; Burtovoi et al. 2020; Illiano et al. 2023; Baglio et al. 2023, 2025a,b, see, however, Jadoliya et al. 2026). Still, many details remain to be understood. In particular, it is unclear what mechanism exactly forces the resulting emission to pulsate. In the Papitto et al. (2019) interpretation, two synchrotron-emitting spots travel along the inner disk surface as the NS rotates, and while the emission from one half of the shock (that is closer to an observer) is absorbed by outer parts of the disk, the emission from the other half is sinusoidally modulated with the period of P_s/2. On the other hand, Veledina et al. (2019) noted that due to transparency of the emitting region the observer basically sees the same synchrotron source all the time, regardless of the rotation phase. They propose instead that the observed modulation is caused by the anisotropic angular distribution of the radiation emerging as a consequence of the presence of the regular large-scale magnetic field.

The observed properties of the pulsed emission allowed to estimate physical conditions inside the volume (e.g. Papitto et al. 2019; Veledina et al. 2019). Here, we briefly reproduce these estimates since they are important for discussion of the transformation of the pulse profile to the single-peaked state. First of all, the injected energy should be radiated sufficiently fast. The characteristic cooling time scale of accelerated electrons (e.g. Papitto et al. 2019) is expressed as

$$\begin{array}{c}\hfill t_{\mathrm{sync}}=\frac{\gamma m_{\mathrm{e}}{c}^2}{W_{\mathrm{sync}}}=\frac{3m_{\mathrm{e}}^3{c}^5}{2e^4{B}^2\gamma },\end{array}$$(5)

where (Ginzburg & Syrovatskii 1965)

$$\begin{array}{c}\hfill W_{\mathrm{sync}}=\frac{2e^4{B}^2}{3m_{\mathrm{e}}^2{c}^3}{\gamma }^2\end{array}$$(6)

is the synchrotron emission power. Here, c is the speed of light, e and m_e are the electron charge and mass, B is the magnetic field strength, and γ is the Lorentz factor. As long as t_sync is less than roughly a quarter of the pulsar period,

$$\begin{array}{c}\hfill t_{\mathrm{sync}}\lesssim \frac{P_{\mathrm{s}}}{4},\end{array}$$(7)

double-peaked pulsations are observed (provided that the size of the emitting region is small compared to cP_s). On the other hand, the peak frequency of radiation from an electron with the Lorentz factor γ is (Ginzburg & Syrovatskii 1965)

$$\begin{array}{c}\hfill \nu (\gamma )=0.29\frac{3eB{\gamma }^2}{4\pi m_{\mathrm{e}}c}.\end{array}$$(8)

These two conditions on the cooling time scale (Eq. (7)) and the emission frequency (Eq. (8)) lead to the following expressions for γ and B,

$$\begin{array}{c}\hfill \{\begin{array}{cc}& \gamma \lesssim {\gamma }_{\max }\simeq 3.3\frac{e^{2/3}{P}_{\mathrm{s}}^{1/3}{\nu }^{2/3}}{m_{\mathrm{e}}^{1/3}c},\hfill \\ \hfill & B\gtrsim B_{\min }=B({\gamma }_{\max })\simeq 14\frac{m_{\mathrm{e}}c\nu }{e}\frac{1}{{\gamma }_{\max }^2}>1.4\frac{m_{\mathrm{e}}^{5/3}{c}^3}{e^{7/3}{P}_{\mathrm{s}}^{2/3}{\nu }^{1/3}}.\hfill \end{array}\end{array}$$(9)

Solving this system for the period of J1023, P_s = 1.69 ms, and an optical frequency of ν = 6 × 10¹⁴ Hz yields γ_max ≃ 60 and B_min ≃ 10⁵ G. The latter is consistent with the magnetic field strength close to the light cylinder, B = B_s(R_ns/R_lc)³, where B_s ≃ 10⁸ G is the surface magnetic field of the NS, R_ns ≃ 10 km is its radius, and R_lc = cP_s/2π ≃ 80 km is the radius of the light cylinder of J1023.

The single optical electron luminosity is L ∼ 0.1 erg/s (as follows from Eq. (6) with γ = γ_max and B = B_min) and the total optical pulsed component luminosity L_p ∼ 10³¹ erg/s (e.g. Papitto et al. 2019); thus, the number of emitting particles N_e ∼ 10³². Since the size of the emitting region, l, has to be smaller than c P_s/4 ∼ 10⁷ cm for pulsations to be observed, the number density is greater than 10¹¹ cm⁻³.

Pulsed emission is observed from optical up to ∼45 keV X-rays (Papitto et al. 2019), which translates to Lorentz factors of emitting electrons from ∼60 to ∼8000, assuming the magnetic field strength of ∼10⁵ G. This energy is deposited to the inner parts of the accretion flow by the pulsar wind. The spin-down luminosity of J1023 in the radio pulsar state was 4.3 × 10³⁴ erg s⁻¹ (Deller et al. 2012). Since, most of the time, the disk is located outside the light cylinder, we can expect that this value would not change significantly with the transition of the system to the subluminous disk state. Jaodand et al. (2016) used X-ray observations obtained with XMM-Newton in 2013–2015 to find that the spin-down rate of J1023 increased by 26.8(4) percent after the transition to the subluminous disk state. Moreover, basing on optical observations of the system obtained with Aqueye+ in 2018–2020, Burtovoi et al. (2020) reported that the spin-down rate is only ∼5 percent faster than that in the radio pulsar state. These small changes can be attributed to quite rare episodes of interaction between the accretion flow and closed magnetosphere, when the disk makes its way beneath the light cylinder. Veledina et al. (2019) showed that even in the case of an isotropic wind, the observed high-mode X-ray luminosity of ∼3 × 10³³ erg s⁻¹ can be powered by the pulsar rotation energy losses under the reasonable assumption that the half-thickness of the inner edge of the disk is about 0.1 of its radius. The pulsed X-ray emission is ≲10³² erg s⁻¹, while the pulsed optical emission is an order of magnitude lower. At a reasonable radius of the inner edge of the disk of 2R_lc (Papitto et al. 2019), the emission region of a linear size of 10⁷ cm constitutes about 0.1 of the full circumference, so the spin-down power is enough also to explain the pulsed luminosity coming from a relatively compact volume.

4.2. Transparency conditions

Some information, in principle, can be obtained from the spectral index α of the synchrotron emission. If the emission is isotropic, then the electron spectral emissivity per unit volume J(ν), on the one hand, is

$$\begin{array}{c}\hfill J(\nu )=\frac{4\pi D^2{F}_{\nu }}{V},\end{array}$$(10)

where V ≈ l³ is the volume of the emission region, D is the distance, and F_ν is the pulsed flux density. On the other hand, assuming the electron energy distribution of the form n(E)dE = κE^−pdE, we can write J(ν) in cgs units (Rybicki & Lightman 1979) as

$$\begin{array}{c}\hfill J(\nu )=\frac{\sqrt{3\pi }{e}^{3}B\kappa }{2c^2{m}_{\mathrm{e}}(p+1)}{\left(\frac{2\pi \nu m_{\mathrm{e}}^3{c}^5}{3eB}\right)}^{-(p-1)/2}a(p),\end{array}$$(11)

where a(p) is a combination of gamma functions on the order of unity. Then, κ can be expressed from Eqs. (10) and (11), while the spectral index p is related to α as p = 2α + 1.

If the density of the medium exceeds some threshold, it becomes opaque to the synchrotron radiation resulting in a strong intensity decrease. This effect may be related to the observed strong variability of the amplitudes of the peaks, with occasional episodes of complete suppression of pulsations. The condition of transparency is expressed as χ_νl ≪ 1, where χ_ν is the absorption coefficient dependent on the radiation frequency (Rybicki & Lightman 1979), expressed as

$$\begin{array}{c}\hfill {\chi }_{\nu }=\frac{\sqrt{3\pi }{e}^{3}B\kappa }{16\pi m_{\mathrm{e}}}{\left(\frac{3eB}{2\pi m_{\mathrm{e}}^3{c}^5}\right)}^{p/2}{\nu }^{-(p+4)/2}b(p),\end{array}$$(12)

where b(p) is another combination of gamma functions on the order of unity. Substituting κ yields

$$\begin{array}{c}\hfill {\chi }_{\nu }={\left(\frac{3Be}{128c{m}_{\mathrm{e}}^3{\nu }^5}\right)}^{1/2}\frac{4\pi D^2{F}_{\nu }}{V}(p+1)\frac{b(p)}{a(p)},\end{array}$$(13)

and for p ≈ 2.4, b(p = 2.4)/a(p = 2.4)≈1, we have

$$\begin{array}{c}\hfill {\chi }_{\nu }l\approx 0.02\times {\left(\frac{B}{{10}^5\text{G}}\right)}^{1/2}{\nu }_{\mathrm{red}}^{-5/2}{D}_{1.37}^2{F}_5{l}_7^{-2}{p}_{2.4},\end{array}$$(14)

where D_1.37 = D/(1.37 kpc) (Deller et al. 2012), ν_red = ν/(5.3 × 10¹⁴ Hz) is the effective frequency of the red detector, F₅ = F_ν/(5 μJy) is the pulsed flux density in the red channel estimated from the optical spectrum of J1023 (Hernandez Santisteban 2016) and the pulsed fraction p_area. Let us rewrite this condition in terms of electron number density. Spectral emissivity per unit volume given by Eq. (11) can also be written as

$$\begin{array}{c}\hfill J(\nu )\approx W_{\nu }\frac{N_{\mathrm{e}}}{l^3},\end{array}$$(15)

where in the frequency range of interest, W_ν ≈ W_sync/ν. Then,

$$\begin{array}{c}\hfill {\chi }_{\nu }l\approx 0.01\times {\left(\frac{B}{{10}^5\text{G}}\right)}^{5/2}\left(\frac{N_{\mathrm{e}}}{{10}^{32}}\right){\left(\frac{\gamma }{60}\right)}^2{\nu }_{\mathrm{red}}^{-7/2}{l}_7^{-2}.\end{array}$$(16)

While the obtained value is rather small, it strongly depends on the characteristics of the emitting volume (in particular, its size l). Thus, in summary, the upper limit on l is set by the pulsar period, while the lower limit follows from the transparency condition, roughly constraining l in the range of 10 − 100 km (for specified B, γ, and N_e).

4.3. Transformation scenarios

Visually inspecting the behavior of the double-peaked and the single-peaked pulse profiles with respect to the rotation phase (see Fig. 7), we note that (1) the phase of the minimum of the single-peaked profile is approximately the same (roughly in the range of 0.55–0.65) as the phase of the minimum between D1 and D2 peaks of the double-peaked profile; and (2) the descending edge of the single-peaked profile is similar to the descending edge of D1 peak of the double-peaked profile (both cover the range of phases of about 0.30–0.55). The shape of the peaks is the result of complex interplay between different factors, such as the time of propagation of the injected energy inside the emitting region, the synchrotron cooling time t_sync, the geometric configuration, and so on. Changes in the physical conditions may affect the shape of D2 peak, stretching its leading and trailing edges and resulting in the observed quick merging of D1 and D2 peaks into S1. From Eqs. (5)–(9), it follows that slight changes in γ and B strongly affect the cooling time scale. We should note, however, that both pulses are generated by roughly the same volume when it is hit by diametrically opposed “lobes” of the pulsar wind, and why the shape of D2 peak is affected by these changes while the shape of D1 is not remains unclear.

Fig. 7. Top: Pulse profile (red and blue channels combined) integrated over both nights with exclusion of the single-peaked interval. “D1” and “D2” labels denote the higher and the lower peaks of the double-peaked profile. Bottom: Average pulse profile of the 220-second-long single-peaked interval. The peak is referred as “S1”. Two phase cycles are plotted for clarity.

Let us consider the changes of the parameters of the emitting plasma within a toy model, assuming that the transformation to the single-peaked state is only governed by the broadening of D2 peak and not by a complete reconfiguration of the structure of the pulsar wind-disk shock. We suppose the total optical pulsed luminosity of ∼10³¹ erg s⁻¹ is apportioned between D1 and D2 as 6 × 10³⁰ erg s⁻¹ and 4 × 10³⁰ erg s⁻¹ in the double-peaked state (roughly according to the relative areas under the peaks), while the luminosity of ∼4 × 10³¹ erg s⁻¹ (since the pulsed fraction is about 4 times higher) is equally apportioned between D1 and D2 (2 × 10³¹ erg s⁻¹ each) in the single-peaked state. These values are somewhat arbitrary due to strong variability of the peaks, but here we are interested in general behavior of γ, B, and n with respect to the observed properties of pulsations. Additionally, we suppose that the trailing edges of the peaks are governed by the synchrotron cooling timescale, t_sync. Then, using Eqs. 6 and 9, we can write γ_max = 786 t_sync^1/3, B_min = 787 t_sync^−2/3 G and L_p = 6.1 × 10⁻⁴N_et_sync^−2/3 erg s⁻¹, where t_sync is given in seconds. In the double-peaked state, t_sync ≃ P_s/4 = 0.42 ms for both D1 and D2 peaks, while in the single-peaked state it stays the same for D1 and roughly doubles for D2. The corresponding parameters of the emitting plasma for the peaks in these two states are listed in Table 2. Our toy model thus suggests that if in the double-peaked state the values of γ and B are close to γ_max and B_min, then moderate variations of these parameters by about 30 percent may lead to broadening of D2 peak. At the same time, these variations does not significantly change the luminosity L_p, which mainly depends on the number of the emitting particles, N_e. The increase by a factor of 6 − 7 in the number densities may be caused, for instance, by the accumulation of matter in the inner parts of the disk due to temporal increase in the accretion rate, while moderate differences in the values of N_e for D1 and D2 peaks may be related to slightly different sizes of the shocks generated by the opposed lobes of the pulsar wind. Veledina et al. (2019) associate flares with a rapid increase of the disk thickness at the inner edge. In this case a larger fraction of the pulsar wind is intercepted by the disk, which may lead, among other things, to changes of physical conditions inside the emitting volume. However, the pulsed fraction is expected to drop in such events (as was observed during a strong 3-ks flare in Papitto et al. 2019). The mechanism for a five-time increase of the pulsed fraction during the single-peaked interval still remains to be understood.

Transformation of the pulse profile in principle could be related to global reconfiguration of the magnetosphere of the NS. Such events have been discussed in the context of magnetar outbursts (e.g., Eichler 2002). A complex pulse profile behavior was observed in magnetar flares and afterglows (Mazets et al. 1999; Feroci et al. 2001; Palmer et al. 2005). Recently, peculiar mode changes and pulse profile phase shifts of the γ-ray pulsar PSR J2021+4026 were associated with reconfigurations of the dipole and quadrupole components of the magnetic field (Razzano et al. 2023; Fiori et al. 2024). However, it is difficult to conceive that something similar could happen in a very old NS with a much lower magnetic field on the order of 10⁸ G. Some kind of interaction between the charged particles produced in the pulsar magnetosphere and the portion of the infalling plasma that managed to reach the closed magnetosphere and penetrate into it probably could change the morphology of the pulsar wind without the reconfiguration of the magnetosphere itself. However, this scenario remains rather unclear and requires further investigation.

5. Summary

We report the discovery of the rapid transformation of the optical pulse profile of the transitional MSP J1023 from its normally observed double-peaked shape to the single-peaked and back again, which we were able to precisely trace in our high temporal resolution observations. The transformation occurred on a timescale of a few seconds or less. At the start of the single-peaked interval, the pulsed fraction exceeded 5 percent, which is almost an order of magnitude larger than the average pulsed fraction in the normal state with the double-peaked pulse profile. Toward the end of this peculiar 220-second-long event, the pulse profile started to take its normal shape and it was only afterward (with a few tens of seconds delay) that the pulsed fraction dropped back to the level of ≲1 percent. Strong flaring activity was observed during this event. Unfortunately, the lack of simultaneous X-ray observations makes it very difficult to understand the processes going on in the system during this change in state.

The shock near the light cylinder of the NS, where the pulsar wind hits the infalling matter, is usually considered as the source of the pulsed emission. The interaction of diametrically opposed lobes of the pulsar wind with the inner radius of the accretion disk gives rise to two emission peaks per one rotation of the NS. At the same time, the details of this scenario remain rather unclear. An aspect that is even more puzzling is the temporary transformation of the pulse profile presented in this paper. It might be related to changes of the physical conditions inside the emitting volume resulting in merging of the two peaks into one. More exotic scenarios include global reconfiguration of the NS magnetosphere or pattern change of the pulsar wind.

Our study reveals one more enigmatic feature of J1023. We stress the importance of further monitoring of the system, especially with regard to the need to carry out coordinated X-ray and optical timing and spectral observations that will make it possible to perform more detailed analyses of similar events in the future.

Data availability

Our optical data are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/709/A170

Acknowledgments

We dedicate this publication to the memory of our dear colleague, Yuri Shibanov, who passed away on September 14, 2025. His contributions to pulsar research in various wavelength ranges, both theoretical and observational, and in particular to the study of PSR J1023+0038, are exceptionally important in modern astrophysics. He also possessed remarkable human qualities – openness, kindness, willingness to share his knowledge and to support his students and colleagues. We are grateful to the anonymous referee for important remarks that helped to significantly improve the quality of the paper. SVK was supported by the EU and by the Czech Ministry of Education, Youth and Sports (project CZ.02.01.01/00/22_008/0004632 – FORTE). The work of AST was supported by the Russian Science Foundation Grant No. 24-12-00320. The work of GMB and VLP was carried out within the framework of the state assignment of the Special Astrophysical Observatory of the Russian Academy of Sciences, approved by the Ministry of Science and Higher Education of the Russian Federation. DAZ was supported by the baseline project FFUG-2024-0002 of the Ioffe Institute. This study is based on the data obtained at the unique scientific facility the Big Telescope Alt-azimuthal of SAO RAS equipped with the MANIA complex. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. DAZ thanks Pirinem School of Theoretical Physics for hospitality.

References

Appendix A: Swift/XRT observations

Properties of X-ray emission of J1023 is one of the best probes of its current state. In the absence of simultaneous X-ray observations we are not able to compare the behavior of the system in the optical and X-ray bands for our data set and, in particular, during the puzzling episode of the structural change of the pulse profile. However, it is worth checking the state of the system near our epoch.

Swift/XRT observed J1023 in the Photon Counting mode on November 28 (ObsID 00033012183), two weeks after our observations. The count rate during the 1.1 ks exposure was 0.131 ± 0.013 s⁻¹ in the 0.3–10.0 keV energy range. Inspection of the light curve (Fig. A.1) hints that after first ∼400 seconds the system probably switched from the high to the low mode.

Fig. A.1. Swift/XRT light curve in the 0.3–10 keV range. Zero time corresponds to 28 Nov 2017 01:16:06 UT, two weeks after our observations. The time bins are 100 seconds wide. Apparently, a switch from the high to the low mode occurred at around 400 seconds from the start.

We extracted the spectrum in the 0.3–10 keV range using the Swift/XRT data products generator (Evans et al. 2009), and fitted it with an absorbed power law model using the X-Ray Spectral Fitting Package (XSPEC) v.12.11.1 (Arnaud 1996). The hydrogen column density along the line of sight cannot be reliably constrained from this observation due to the low count statistic. We fixed it at the value of N_H = 3.1 × 10²⁰ cm⁻² reported by Bogdanov et al. (2015) for the high mode. We obtained a good fit with the spectral photon index Γ = 1.8 ± 0.2 and the unabsorbed flux F_X = (4.6 ± 0.8)×10⁻¹² erg cm⁻² s⁻¹ in the 0.3–10 keV range, which at the parallax distance of 1.37 kpc (Deller et al. 2012) translates to the luminosity L_X = (1.0 ± 0.2)×10³³ erg s⁻¹. We conclude that the system was in its normal subluminous disk state during the Swift/XRT observation performed two weeks after our optical observation.

All Tables

All Figures

Fig. 1. Z22-test for both nights joined. The orbital phase correction is shown with respect to the ephemeris value obtained from the timing solution by Jaodand et al. (2016). The maximum at 0.00129 corresponds to the shift of the epoch of passage at the ascending node by −22 seconds with respect to 54905.97140075 MJD. The aliases around the main peak emerge due to the presence of gaps in each data set and between them.

In the text

	Fig. 2. Top: Typical double-peaked pulse profiles obtained by phase folding together the Nov 14 and Nov 15 light curves and excluding the single-peaked interval. The red and blue channels are shown on the left and right, respectively. Bottom: Averaged single-peaked pulse profiles of the 220-second-long interval. Corresponding pulsed fractions parea and prms are shown in each panel.
In the text

	Fig. 3. Segment of the light curve in the red (top) and blue (bottom) channels containing the interval with single-peaked pulsations (shown with a hatched stripe). Flaring activity around the interval is readily seen. The time resolution is 1 s. Zero time corresponds to 15 Nov 2017 02:37:29 UT.
In the text

Fig. 4. Pulse profile variability during the first four 10-second-long segments (a, b, c, and d) of the single-peaked interval. Profiles observed within the time ranges indicated in the legends are plotted in black. Time is counted from the onset of the single peak. Pulsed fractions and detection significance levels are shown in each panel. Blue and orange filled step plots show the averaged pulse profile of the full single-peaked interval and the averaged over both nights typical double-peaked pulse profile, respectively.

In the text

	Fig. 5. Same as Fig. 4, but for the first 40 seconds after the single-peaked interval. The pulsed fraction remains high, but the peak is shifted in phase and the pulse profile starts to resemble the normal double-peaked one, although the lower peak is lost in noise at such a short integration time.
In the text

	Fig. 6. Top to bottom: Pulse profiles (red and blue channels combined) at different integration times (indicated in the legends) from the onset of the single-peaked interval. Corresponding background-corrected pulsed fractions and significance levels are given in each panel.
In the text

	Fig. 7. Top: Pulse profile (red and blue channels combined) integrated over both nights with exclusion of the single-peaked interval. “D1” and “D2” labels denote the higher and the lower peaks of the double-peaked profile. Bottom: Average pulse profile of the 220-second-long single-peaked interval. The peak is referred as “S1”. Two phase cycles are plotted for clarity.
In the text

	Fig. A.1. Swift/XRT light curve in the 0.3–10 keV range. Zero time corresponds to 28 Nov 2017 01:16:06 UT, two weeks after our observations. The time bins are 100 seconds wide. Apparently, a switch from the high to the low mode occurred at around 400 seconds from the start.
In the text