© The Authors 2025

1. Introduction

PSR J1023+0038 (J1023) is the archetypal transitional millisecond pulsar (tMSP). It was classified as a tMSP following the detection of radio pulsations revealing a 1.69-ms spin period neutron star (NS) in a 4.75-hour binary system (Archibald et al. 2009). In 2013, J1023 underwent a sudden multi-wavelength flux increase, accompanied by the disappearance of its radio pulsations and the appearance of optical emission lines. This change marked a transition to a peculiar state characterised by the presence of an accretion disc and an average X-ray luminosity of 10^{33 − 34} erg s⁻¹, known as the X-ray sub-luminous state (Stappers et al. 2014). In this state, J1023 shows three X-ray luminosity modes (Linares 2014; Archibald et al. 2015; Bogdanov et al. 2015; Coti Zelati et al. 2018; Papitto et al. 2019): high (∼7 × 10³³ erg s⁻¹, 70–80% of the time); low (∼10³³ erg s⁻¹, 20–30%); and flaring (≳10³⁴ erg s⁻¹, 1–2%). X-ray, UV, and optical pulsations at the NS spin period are detected only in the high mode and in flaring mode, even though at low significance (Papitto et al. 2019; Miraval Zanon et al. 2022; see also Ambrosino et al. 2017; Jaodand et al. 2021). Radio and millimetre observations reveal bright continuum emission with an anti-correlation between the radio and X-ray bands: higher radio flux is observed during the low mode, with the spectrum shifting from flat (α ∼ 0, where α is the spectral index) to inverted (α > 0; Bogdanov et al. 2018). The radio emission is attributed to synchrotron radiation originating from a Blandford & Königl (1979)-type ‘compact jet’ in the high mode and discrete ejecta in the low mode (Fender et al. 2004, see also Baglio et al. 2023).

To explain the behaviour of J1023, recent scenarios propose the presence of an always-active radio pulsar (Papitto et al. 2019; Veledina et al. 2019; Baglio et al. 2023). In the high mode, pulsed X-ray/UV/optical emission originates from synchrotron-emitting regions in a shock whereby the pulsar wind interacts with the inner accretion flow, ∼100–200 km from the pulsar. In the low mode, the shock may move outwards, reducing the X-ray flux and causing the disappearance of optical-to-X-ray pulsations (Papitto et al. 2019). Alternatively, during low modes, the in-falling matter may penetrate within the light cylinder, pushing J1023 into a propeller regime with lower dissipation, reducing the X-ray luminosity and suppressing pulsations (Veledina et al. 2019). Baglio et al. (2023) detected rapid millimetre flares during the high-to-low mode switch; they attribute them to the ejection of the inner accretion flow as optically thin ejecta, leading to the observed increase in the radio flux during the low mode.

In the high mode, synchrotron emission from the shock contribute ≃5% of the optical flux and ≃90% of the soft X-ray flux (Baglio et al. 2023). The shock region emits polarised radiation at levels of ≈16% in X-rays (in the high mode) and ≈1% in the optical band (on average), with aligned polarisation angles. For low- and flaring-mode data, X-ray polarisation upper limits of 26% and 28%, respectively, were derived. Radio polarisation was not detected, with upper limits of ∼10% in the high mode and ∼19% in the low mode (Baglio et al. 2024).

Despite extensive studies on low and high modes, the flaring mode is less explored, and most research has focused solely on its X-ray evolution. Current scenarios (Papitto et al. 2019; Veledina et al. 2019) suggest that the flaring mode is due to an increase in the thickness of the inner accretion disc intercepting the pulsar wind, though clear observational evidence is lacking. This Letter presents the first multi-band study of the flaring mode of J1023 aimed at providing a comprehensive understanding of the behaviour of tMSPs in this emission mode in the frame of the scenarios described above.

2. Data analysis

Table A.1 reports the observations presented in this work. Data were barycentered at the Solar System using the DE440 ephemeris and the position: RA = 10^h23^m47^s69; Dec = +00°38′40$\stackrel{″}{.}$8 (J2000.0; Deller et al. 2012).

2.1. IXPE

The Imaging X-ray Polarimetry Explorer (IXPE) observed J1023 in two sessions from May 29 to June 14, 2024, totalling ≃675 ks of exposure per detector unit (DU). An in-depth analysis of these data was reported by Baglio et al. (2024). After reducing instrumental background in the data (Di Marco et al. 2023), we selected source photons from a 60″ radius circle, and background photons from an annulus with radii of 150″ and 300″. We extracted the combined background-subtracted light curve from the three DUs using ixpeobssim (Baldini et al. 2022) and FTOOLS.

2.2. Swift

We performed two observations with the Neil Gehrels Swift Observatory (Swift), employing the X-ray Telescope (XRT; Burrows et al. 2005) in photon-counting mode and the Ultra-Violet/Optical Telescope (UVOT; Roming et al. 2005) in event mode with the UVM2 filter. We screened XRT data using standard procedures and collected source photons from a circle with a radius of 47.2″ and background photons from an annulus with radii of 94.4″ and 188.8″ centred on the target. We extracted background-subtracted time series in the 0.5–10 keV range, with 50-s bins. For the UVOT data, we used a circle with a 5″ radius to collect source photons, and a circle away from the source with a 10″ radius for background photons. We extracted background-subtracted time series binned at 30 s using the uvotevtlc tool.

2.3. VLT

J1023 was observed with FORS2 on the Very Large Telescope (VLT) at Cerro Paranal, Chile, on June 4–5, 2024, for 7742 s (≃2.15 hours), covering 45% of its orbital period. The observations, split into two sets (3946 s and 3796 s), produced 112 images with 20 s integration each. We obtained normalised Stokes parameters, Q_opt and U_opt, for linear polarisation and we evaluated the intrinsic polarisation level, P_opt, and angle, θ_opt, following Baglio et al. (2020, 2023). Details are given in Appendix B. The photo-polarimetric data also allowed us to extract the system R-band light curve, since the sum of f^o(Φ) and f^e(Φ) for each image gives the total target flux. We did this calculation for J1023 and two field stars to check for any anomaly and used the light curves of the field stars to perform differential photometry. The resulting light curve of J1023 is shown in Fig. 1.

Fig. 1. Multi-band time series and evolution of optical and radio parameters (including data presented by Baglio et al. 2024). Error bars represent the 1σ confidence levels. In the seventh and eighth panels, horizontal dashed red lines mark the average optical polarisation degree (Popt) and angle (PAopt), with shaded areas indicating one standard deviation about the mean. In the ninth panel, a horizontal dotted line shows a spectral index of 0 (i.e., a flat spectrum). In the final panel, downward arrows denote 3σ upper limits on the radio polarisation degree (Pradio). TDB stands for barycentric dynamical time.

2.4. VLA

Two radio observations of J1023 were conducted in the C band (∼4–8 GHz) with the Karl G. Jansky Very Large Array (VLA). The first observation was conducted on June 4, 2024, simultaneously with IXPE, Swift, and the VLT. The second observation was conducted on September 15, 2024. The observations were calibrated with the Common Astronomy Software Applications VLA pipeline (casa v6.5, CASA Team 2022). The June observations included polarisation calibrators; as a result, following the pipeline, we manually solved for polarisation calibration solutions. We discuss the details of the observations, imaging, and data analysis in Appendix C.

3. Results

3.1. Multiband photometry

J1023 entered the X-ray flaring mode during an observation gap between barycentric modified Julian day (BMJD) 60466.025 and 60466.052. The first IXPE data point post-gap already has a count rate above the typical values in the low and high modes. After an initial flare lasting ≈15 min, the X-ray emission drops before exhibiting more flares. J1023 then resumes alternating between low and high modes after around BMJD 60466.07483.

During the first VLA observation, J1023 underwent a bright (≳ 1 mJy), rapidly evolving flare that lasted for ∼1 h; before and after the flare, the source was significantly dimmer (∼ 240 μJy; see Figs. 1 and D.1). Notably, the radio flux begins rising around BMJD 60466.03, peaks around BMJD 60466.06, and declines to pre-flare levels by BMJD 60466.09. The radio flare is complex, with substructures suggesting repeated events, and the spectral index shifts from optically thick to optically thin on minute timescales. Qualitatively, the radio flares appear to exhibit a small frequency delay, with the evolution at 5 GHz lagging behind that at 7 GHz (see Fig. B.1, middle panel). To quantify the delay, we calculated the cross-correlation function between the sub-bands using the ZDCF algorithm (Alexander et al. 1997), measuring the uncertainty on the most significant lag with PLIKE (Alexander 2013). The ZDCF analysis favoured a non-zero lag of $0.7_{-0.4}^{+0.3}$ min, consistent with our interpretation of a propagating synchrotron emitting plasma (i.e. a jet), though the significance of the lag is marginal (i.e. ∼1.8σ; Fig. D.2).

Alongside the radio flaring, the optical and UV fluxes increase significantly. The first optical data point after the gap is ∼0.6 mag brighter than before. Similarly to the radio band, the optical flare may peak around BMJD 60466.06, but its shorter duration complicates a comparison with the radio flare. The optical light curve shows two small flares early on, peaking at BMJD 60465.997 and BMJD 60466.005, lasting ∼7 and ∼10 minutes, respectively. The first flare has a radio counterpart (peaking at BMJD 60465.997 and lasting ∼10 min) and corresponds to an X-ray low mode (Fig. 1).

3.2. Multi-band polarimetry

Figure 1 shows the time evolution of the optical and radio polarisation properties of J1023. The average optical polarisation fraction is (1.38 ± 0.03)% before the bright flare, decreasing to (0.49 ± 0.04)% between BMJD 60466.051 (the first flare observation) and BMJD 60466.071, before returning to its pre-flare level during the decay. Pre-flare, the Stokes parameters, Q_opt and U_opt, are approximately constant. U_opt remains constant during the flare, but Q_opt evolves, initially decreasing before returning to its pre-flare value by BMJD 60466.071. The average polarisation angle is PA_opt = −1.0 ° ±0.9° before the flare. Similar to the polarisation fraction, the polarisation angle changes during flaring and returns to its pre-flare value during the flare decay.

No significant radio polarisation was detected in our observations. Imaging on 1-min intervals allowed us to place 3σ upper limits of ≲55% on average on the polarisation fraction before and after the radio flare, and of ≲10% during the flare. Imaging on longer time intervals improved the upper limits to < 8.7%, < 2.3%, and < 8.2%, before, during, and after the radio flare. We observed short-timescale variability of the spectral indices, suggesting that the significance of synchrotron self-absorption varies during the sub-flaring. This is reminiscent of the ‘optically thick flares’ observed in X-ray binaries (XRBs; e.g. Fender & Bright 2019), where the rise is optically thick and the decay is optically thin. Given that optically thin synchrotron radiation is more polarised than its optically thick counterpart, we also imaged by combining the times when the radio flux density was rising and decaying (green and yellow regions; Fig. D.1). In both images, J1023 remains undetected with upper limits of < 3.7% and < 3.6% during rise and decay times.

4. Discussion

The observed X-ray flaring has UV, optical, and radio counterparts. Previous studies reported two types of X-ray flares: short (seconds to minutes), isolated flares, and flares happening over extended periods (hours), with intense optical and UV flares coinciding with both types of X-ray flare (Bogdanov et al. 2015; see also Kennedy et al. 2018; Papitto et al. 2019; Miraval Zanon et al. 2022). Flares as long as 10 hours have also been detected (Tendulkar et al. 2014). Jaodand et al. (2021) observed quasi-periodic flares in X-ray, UV, and optical bands, each with an average duration of ≈500 s. This suggests a mechanism that is distinct from high-low mode switching and that continues to operate during flares alongside additional broadband emission. Conversely, radio flares exhibit a more complex connection with X-ray activity, occurring during low-mode episodes, accompanying X-ray flares, or appearing independently (Bogdanov et al. 2018).

In the June VLA dataset, we detected the brightest (≳1.5 mJy at 6 GHz) and longest (∼50 minutes) radio flaring from J1023 since its 2013 state transition. This prominent radio flare is likely unrelated to high-low mode switches, being much longer and brighter than those typically seen during low modes. For comparison, a smaller radio flare around BMJD 60465.997 during an X-ray low mode lasted ∼5 minutes and reached ∼0.46 mJy. The radio flare displays clear substructures, with multiple sub-flares superimposed on an increased baseline radio flux density. A similar substructure is observed in the optical light curve.

The optical polarisation degree decreases during the flare, and the polarisation angle rotates before gradually reverting to its pre-flare value as the flaring activity subsides (Figs. 1 and B.2). This indicates that a polarised component dominates when the flare component is absent or faded, and that an unpolarised component emerges during the flare. Following Baglio et al. (2024) and noting the alignment of optical and X-ray polarisation angles, we propose that the polarised component originates from the shock between the pulsar wind and the inner accretion flow, which is then dominated by the unpolarised flaring emission.

The origin of the observed multi-band flaring behaviour is not well understood. However, the simultaneity of the X-ray, UV, optical, and radio flaring suggests a common underlying mechanism. In tMSPs, X-ray flares likely result from transient increases in inner accretion disc thickness caused by thermal instabilities or magnetic field rearrangements (Veledina et al. 2019). A thicker disk could intercept a larger fraction of the pulsar wind, enhancing energy dissipation and irradiation and thereby increasing the X-ray luminosity. The enhanced irradiation may boost the population of free particles at the light cylinder radius. If these particles reach the jet, the resulting increase in mass-loading would produce a more luminous jet and, hence, stronger optical and radio emission. Strong short-timescale optical variability is a typical feature of jets in XRBs (e.g. BH XRBs GX 339-4, V404 Cygni, MAXI J1535-571, and MAXI J1820+070; Gandhi et al. 2010, 2017; Baglio et al. 2018; Paice et al. 2019; Tetarenko et al. 2021a). In this scenario, the optical flaring from the weakly polarised compact jet dominates over the shock radiation, leading to depolarisation. The jet likely contributes little to the flaring UV emission, which may arise from the same mechanism as X-ray flares (i.e. thickening of the shock region), increasing dissipation and fluxes. Similarly, optical flares might also originate in the disc rather than the jet. In this interpretation, the decreased optical polarisation may result from less ordered magnetic fields in the shock due to inner disc thickening, while the change in polarisation angle could reflect a change in the projection of the average magnetic field direction of the shock on the sky plane.

The radio variability holds an insight into the nature of the jet(s). Similar to what is observed in J1023, compact jets in XRBs are known to flare on minute-to-hour timescales, resulting in light curves that are delayed at lower frequencies (e.g. Tetarenko et al. 2021b). In one interpretation, the radio variability is due to shocks between discrete shells of material moving at different velocities along the jet axis (due to stochastic fluctuations in the mass-loading from the accretion flow, e.g. Malzac 2014). The delays result from propagation times, whereby different frequencies probe distinct regions along the jet axis (higher frequencies are closer to the jet base).

Alternatively, the spectral variability and frequency delays during the radio flaring resemble the optically thick flares typically seen in XRBs and are commonly associated with the launching of discrete jet ejecta (e.g. Fender & Bright 2019). In this scenario, the radio spectrum evolves from self-absorbed to optically thin during the flare rise and decay. Indeed, the behaviour of the radio spectral indices is broadly consistent with this trend, as the maximum (α ∼ 0.5) spectral index occurred as the flux density rose, and the minimum (α ∼ −0.5) occurred when the flux density decayed. This is supported by our integrated images, in which we combined times when the source was rising and decaying to boost the signal-to-noise ratio (yellow and green regions on Fig. D.1). We observe a statistically significant difference in the spectral indices during the rise (α = 0.18 ± 0.04) and decay (α = −0.16 ± 0.04). Smaller flares from J1023 during low-mode episodes have also been interpreted as discrete ejections (Bogdanov et al. 2018; Baglio et al. 2023).

We searched for evidence of jet ejecta in our September follow-up VLA observations ∼3 months after the initial multi-wavelength campaign. While we did not detect any ejecta, with a 3σ upper limit on the flux density of ∼7.5 μJy, this does not rule out that the observed flaring is connected to jet ejections, as the ejecta may be short-lived, similar to those seen in the BH XRB V404 Cygni (Tetarenko et al. 2017; Miller-Jones et al. 2019). Moreover, J1023 was, once again, exhibiting radio flaring (see Fig. C.1), reaching a maximum flux density of ∼1.2 mJy. As a result, slow-moving ejecta (i.e. ≲20 mas/day) would blend into the bright, flaring core emission, making them undetectable. The absence of polarised radio emission (≲2%) also seems to favour a compact jet origin for the flaring. XRB compact jets regularly exhibit polarisations ≲3% (Han & Hjellming 1992; Corbel et al. 2000; Russell et al. 2015), while ejecta can approach the theoretical limit for optically thin synchrotron radiation (∼50%; Brocksopp et al. 2013; Curran et al. 2014). However, some XRBs launch ejecta with polarisation fractions well below that threshold, such as the jets from V404 Cygni (Han & Hjellming 1992; Hughes et al. 2023).

By associating each optically thick-to-thin spectral transition from J1023 with a single ejection, we can estimate the minimum energy of each ejecta using standard minimum energy calculations under equipartition conditions and assuming that the flare spectrum evolution is driven by decreasing optical depth to synchrotron self-absorption (Fender & Bright 2019). We identify at least five individual flares with such transitions (Fig. D.1). We estimate minimum energies of ≃(1.4 − 2.2)×10³⁵ erg, corresponding to a minimum total kinetic feedback of ≈10³⁶ erg within 2 hours, and an equipartition magnetic field of ≃1.4 G per ejection. Assuming each ejection has a duration consistent with the flare rising timescale (e.g. Carotenuto et al. 2024), ≃300 s for J1023, we estimate a jet power of ≃(5 − 8)×10³² erg s⁻¹ (∼1–10% of the X-ray luminosity during flares; Bogdanov et al. 2015).

The radio properties of J1023 closely resemble those of V404 Cygni during its 2015 outburst, including short-timescale (< hour) radio flaring (Tetarenko et al. 2017), rapidly evolving spectral indices (Fender et al. 2023), and very low polarisation levels (0.05–0.3% at 6 GHz; Hughes et al. 2023). VLBI confirmed that V404 Cygni’s flaring was linked to short-lived jet ejecta (Miller-Jones et al. 2019), with multi-wavelength modelling revealing lower jet power than typical XRBs (∼(4 − 3000)×10³² erg s⁻¹; Tetarenko et al. 2017). J1023’s jet power overlaps with the lower range of V404 Cygni’s ejecta. VLBI observations of J1023’s flaring could confirm jet ejecta, with follow-up studies investigating how J1023 and other tMSPs power ejections despite differing geometry and lower accretion power (see Appendix D).

5. Conclusions

Our observations of J1023 have enabled the most comprehensive study of the flaring mode of a tMSP in a sub-luminous state to date. The X-ray flaring has UV and optical counterparts and coincides with the brightest radio flare observed for J1023 after the 2013 state transition. During the flare, the optical polarisation drops, indicating the emergence of an unpolarised component. While X-ray and UV flares may result from the thickening of the inner accretion disc, we argue that the jet could drive the optical and radio flaring. The lack of correlation between the UV and optical light curves supports this interpretation. The larger dissipation may generate more free electrons, contributing to loading the compact jet. The radio properties resemble compact jet flaring seen in XRBs, though the variability and spectral evolution suggest that discrete ejecta may also contribute. The evolution of the radio-band spectral index during flares is consistent with synchrotron self-absorption in jet ejections or shocks in the compact jet. The derived radio polarisation upper limits (< 8.7%, < 2.3%, and < 8.2%, before, during, and after the radio flare) favour a compact jet origin but do not rule out discrete ejecta.

The detection of J1023 in a bright radio state during two brief observations conducted months apart suggests that the source may have entered a previously unobserved phase characterised by frequent radio flaring. Continued monitoring at radio wavelengths could provide further insights into this behaviour. VLBI observations during subsequent flaring may conclusively identify whether the flares result from the launching of discrete ejecta. If confirmed, tMSPs may become crucial laboratories for probing jet launching mechanisms in distinct accretion regimes.

Data Availability

The reduced FORS2/VLT data are available at DOI: 10.5281/zenodo.14840659. Raw data can be obtained from ESO archive upon request (http://archive.eso.org/). This research used data products provided by the IXPE Team (MSFC, SSDC, INAF, and INFN) and distributed with additional software tools by the High-Energy Astrophysics Science Archive Research Center (HEASARC), at NASA Goddard Space Flight Center (GSFC). The data that support the findings of this study are publicly available at their respective online archive repositories (Swift: https://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/w3browse.pl; VLA: https://data.nrao.edu/portal/).

Acknowledgments

MCB is supported by the INAF-Astrofit fellowship. FCZ is supported by a Ramón y Cajal fellowship (RYC2021-030888-I). AKH is supported by NSERC Discovery Grant RGPIN-2021-0400. FC is supported by the Royal Society through the Newton International Fellowship programme (NIF/R1/211296). SC and PD’A are supported by ASI grant I/004/11/5. ADM and FLM are supported by the Italian Space Agency (ASI) via contract ASI-INAF-2022-19-HH.0. AP and DdM are supported by ASI and INAF under agreements ASI-INAF I/037/12/0, ASI-INAF n.2017-14-H.0, INAF Sostegno alla ricerca, and INAF SKA/CTA projects. AP acknowledges support from the Fondazione Cariplo/Cassa Depositi e Prestiti, grant no. 2023-2560. FCZ, SCa, PD’A, AP, DdM, and GI are supported by the INAF-FANS project and MUR PRIN 2020 grant GEMS 2020BRP57Z. DMR and KA are supported by Tamkeen under the NYU Abu Dhabi Research Institute grant CASS. DFT is supported by the grant PID2021-124581OB-I00 funded by MCIU/AEI/10.13039/501100011033 and 2021SGR00426. GI is supported by the AASS Ph.D. program with INAF collaboration. NR is funded by ERC Horizon 2020 grants MAGNESIA (No. 817661) and DeepSpacePulse (No. 101189496). This work is supported by the Spanish Unidad de Excelencia María de Maeztu CEX2020-001058-M and MCIU with EU NextGeneration EU funding (PRTR-C17.I1). Based on observations collected under ESO programme 113.27RE.001.

References

Appendix A: Journal of the observations

Table A.1 reports details on the observations presented in this work.

Appendix B: Polarimetry with VLT/FORS2

The optical polarimetric analysis presented in this work follows that described by Baglio et al. (2020). The VLT/FORS2 polarimeter works with a Wollaston prism that is inserted in the light path and splits the incoming light into two orthogonally polarised beams, and a rotating half-wave plate (HWP) captures images at four angles relative to the telescope axis (Φ_i = 22.5 ° (i − 1) for i = 1, 2, 3, 4). For our dataset, this procedure was repeated 28 times. Images were reduced by subtracting an average bias frame and dividing by a normalised flat field; using daophot (Stetson 1987), aperture photometry was performed with a 6-pixel aperture. Once aperture photometry was performed on all stars in the field using daophotStetson (1987), we extracted the normalised Stokes parameters Q_opt and U_opt for linear polarisation as:

$$\begin{array}{c}\hfill Q=\frac{F({\mathrm{\Phi }}_1)-F({\mathrm{\Phi }}_3)}{2};U=\frac{F({\mathrm{\Phi }}_2)-F({\mathrm{\Phi }}_4)}{2},\end{array}$$(B.1)

where

$$\begin{array}{c}\hfill F({\mathrm{\Phi }}_i)=\frac{f^o({\mathrm{\Phi }}_i)-f^e({\mathrm{\Phi }}_i)}{f^o({\mathrm{\Phi }}_i)+f^e({\mathrm{\Phi }}_i)},\end{array}$$(B.2)

with f^o and f^e being the ordinary and extraordinary beams fluxes, respectively (see Sect. 2.3). The calculated Stokes parameters are shown in Fig. B.1, and are not corrected for instrumental polarisation. This correction can be achieved by observing a non-polarised standard star, as detailed in the FORS2 instrument documentation¹. For FORS2, regular monitoring of unpolarised standard stars ensures accurate measurement of the instrumental polarisation. Over the past decade, these observations have demonstrated that the instrumental polarisation remains consistently low, within the range of 0–0.3% across all bands².

Fig. B.1. Stokes parameters Qopt and Uopt for the optical emission. Data points are colour-coded according to the epoch of their measurements.

We can then calculate the parameter S(Φ), which is the component of the normalised Stokes vector that describes the LP along the direction Φ, for each HWP angle (Φ_i = 22.5 ° (i − 1) for i = 1, 2, 3, 4; Sect. 2.3) as

$$\begin{array}{c}\hfill S(\mathrm{\Phi })=(\frac{f^o(\mathrm{\Phi })/f^e(\mathrm{\Phi })}{f_u^o(\mathrm{\Phi })/f_u^e(\mathrm{\Phi })}-1){(\frac{f^o(\mathrm{\Phi })/f^e(\mathrm{\Phi })}{f_u^o(\mathrm{\Phi })/f_u^e(\mathrm{\Phi })}+1)}^{-1},\end{array}$$(B.3)

where f^e(Φ) and f^o(Φ) are the extraordinary and ordinary fluxes at angle Φ, and f_u^e(Φ) and f_u^o(Φ) are the corresponding fluxes for non-polarised standard stars (see Baglio et al. 2024).

S(Φ) relates to the polarisation level P_opt and angle θ_opt as S(Φ) = P_opt cos2(θ_opt − Φ). P_opt and θ_opt are derived using a two-step process (see Baglio et al. 2020, 2023): initial estimates via Gaussian likelihood maximisation, and refinement of the values using Markov Chain Monte Carlo. The final values are taken as the median of the marginal posterior distributions, with uncertainties from the 16th - 84th percentiles. This method corrects for instrumental and interstellar polarisation, since S(Φ) is normalised using non-polarised stars in the field (in this case, we selected 3 isolated stars in the field to enhance statistics). Given the low line-of-sight absorption and proximity of J1023, interstellar polarisation is likely well-subtracted with this method. A correction of 1.6 °   ±  0.7° was applied on the polarisation angle using the polarised standard star Vela 1-95 (see Baglio et al. 2024).

Figure B.2 shows the trend of P_opt, PA_opt, Q and U as a function of the R-band magnitude of J1023. A Spearman correlation test between P_opt and magnitude yielded a correlation coefficient of ρ = −0.78 and p-value of 3 × 10⁻⁸ (see top panel of Fig. B.2). Using a Bayesian MCMC method with the LINMIX₋ ERR model (Kelly 2007), we derived a median slope of β = −1.0 ± 0.2 and an intrinsic scatter of ${\sigma }_0=0.{07}_{-0.02}^{+0.04}$ from 10000 posterior samples. The middle panel of Fig. B.2 shows that while U_opt remains roughly constant (around 0) during the flare, Q_opt is higher when the flux is lower. These results confirm that an unpolarised component becomes dominant during the flare.

Fig. B.2. Polarisation degree (top), Stokes parameters (middle), and polarisation angle (bottom) as a function of the relative intensity for the optical emission from J1023. The top panel shows the results of the Spearman rank correlation test and of a linear regression algorithm to the data. The shape of the anti-correlation is consistent with ΔPopt / Δmag = -1.0±0.2, with an intrinsic random scatter around the best regression model of ${\sigma }_0=0.{07}_{-0.02}^{+0.04}$ (see text for details).

Appendix C: VLA observations

As in Baglio et al. (2024), the VLA observations used the 3-bit sampler, dividing the 4096 MHz bandwidth into 32 spectral windows with 64 2-MHz channels each. 3C286 served as the flux calibrator, J1024−0052 as the gain calibrator, and J1407+2827 as the leakage calibrator. J1023 was observed over 21 ∼9-min scans, totalling ∼3.15 hr on-source.

We first processed the data using the casa VLA pipeline (v6.5; CASA Team 2022), removing and residual radio frequency interference (RFI) through a combination of manual editing and auto-flagging routings rflag and tfcrop. Currently, the VLA pipeline only calibrates the parallel-hand visibilities (i.e. Stokes I and V). To measure the complete polarisation state, we performed additional calibration of the cross-hand visibilities (i.e. Stokes Q and U) following the standard procedures³. Stokes IQUV images were created using wsclean (v3.4; Offringa et al. 2014), applying a CLEAN-based deconvolution (Högbom 1974). Phase-only self-calibration was performed using quartical (Kenyon et al. 2022). The field of J1023 includes a bright radio galaxy (J102358.2+003826) 3″ away and other, fainter, point sources. We subtracted the visibilities of these background sources, allowing for narrow field-of-view images with minimal aliasing.

Given J1023’s complex evolution, we created a suite of images to probe the range of phenomena; these included 20-second and 1-minute interval images to capture short-timescale variability and integrated images (i.e. by combining independent scans) for better signal-to-noise before, during, and after the flare. For each imaging interval, wsclean produced 16 channelised images and a frequency-averaged Multi-Frequency-Synthesis image. To measure the frequency-dependant time lags during the radio flaring, we also created a Multi-Frequency-Synthesis image for the lower (4–6 GHz) and upper (6–8 GHz) sub-bands. The Stokes I flux density was measured using imfit, and for strong detections (> 20σ), an intra-band spectral index was calculated. We also created linear polarisation images $(P_{\mathrm{radio}}=\sqrt{Q_{\mathrm{radio}}^2+U_{\mathrm{radio}}^2})$. No polarised emission was detected, so we calculated 3σ upper limits following Vaillancourt (2006).

We obtained an additional VLA observation of J1023 to check for evidence of propagation from any ejecta that may have been launched during the June flaring. The second VLA observation (DDT 24B-469) was conducted on September 15, 2024, between 14:51:18 and 17:29:27 UTC (∼2 hr on-source) with the same setup (C band, B configuration) and the same processing and reduction as for the June observation. However, the September observations did not include a leakage calibrator, so we exclusively analysed the Stokes I properties. The total exposure time on J1023 was 128 minutes, and the rms noise of the full observation was ∼2.5 μJy. We obtained a clear detection of J1023 at ∼700 μJy while it was displaying additional flaring activity, as shown in Fig. C.1. The flare shape and duration could be consistent with increased mass-loading of the compact jet or the launch of a discrete ejection, though current data cannot fully distinguish between these scenarios.

Fig. C.1. Top: VLA 6 GHz light-curve of J1023 from the September 2024 observation, with 1-minute time binning. The short-time scale variability from the flaring activity is evident. Bottom: in-band spectral index with 1-minute time bins. The radio spectrum shows transitions between optically thick and thin phases on timescales of minutes.

Appendix D: Jet radiative luminosity

As discussed in Sect. 4, we estimate the jet power to be ≃(5 − 8)×10³² erg s⁻¹, assuming that each flare from the first VLA observation results from a single ejection. Alternatively, if the radio flaring activity is due to increased mass loading in the compact jet, we can tentatively estimate the jet radiative luminosity by assuming a typical compact jet spectrum (Blandford & Königl 1979; Russell et al. 2020). To perform this estimation, we integrate a flat spectrum (consistent with the average flare spectral index) between 10⁹ Hz and the break frequency ν_break = 2.5 × 10¹³ Hz (as estimated by Baglio et al. 2023), anchoring it at the VLA flux measured at 6 GHz. Above the break frequency, we assume a typical optically thin spectrum with a spectral index of α = −0.6 extending up to 10¹⁸ Hz. We note that our results are only weakly dependent on this assumption.

For the peak flux density of ∼1.5 mJy during the flare, we obtain a total jet radiative luminosity of ≃1.5 × 10³⁴ erg s⁻¹. We caution that this estimation is subject to significant uncertainty, as there are currently no measurements of the break frequency during the flaring mode. Notably, the jet luminosity in this mode is comparable to the source’s accretion power in the high mode, GM Ṁ/R_LC ≈ 10³²−10³⁴ erg s⁻¹, estimated considering an accretion rate of Ṁ ≃ 10¹³⁻⁵ × 10¹⁴ g s⁻¹ and a light cylinder radius of R_LC = 80 km (Linares et al. 2022). For the radio flux before the flare (∼200 μJy), we estimate a total jet radiative luminosity of ≃2 × 10³³ erg s⁻¹. This value is again compatible with the accretion power in the high mode, considering the large (orders of magnitude) uncertainty in the mass accretion rate.

While it would be valuable to disentangle the contributions of different emission components through a spectral energy distribution, this is not feasible for the flaring mode. Capturing entire radio and optical flares would require long averaging, which would wash out spectral information on the radio flare substructures and combine X-ray data from different emission modes. This would most likely result in an oversimplified or misleading interpretation, as previous studies indicate that distinct emission components are very likely to be active in different modes.

Expanding on our comparison to V404 Cygni, the accretion luminosity from J1023 is five orders of magnitude lower than V404 Cygni in 2015 which, at the time of flaring, was accreting near the Eddington luminosity (∼ 10³⁹ erg s⁻¹ for a ∼10M_⊙ black hole; Motta et al. 2017), and references therein). In contrast, the jet power in each flare (assuming jet ejecta) is only separated by a factor of ≲3000; while very tentative, these results suggest that tMSPs in the flaring state may be more efficient at powering jet ejections than XRBs in their comparably luminous states (i.e. ‘the very high state’; Done et al. 2007).

Fig. D.1. Top: VLA 6 GHz light-curve of J1023 from the June 2024 observation during the flaring period, with 1-minute time binning. The radio emission shows multiple flares of similar duration, with green and yellow regions indicating flux rise and decay phases, respectively. Second panel: The same light-curve split into two sub-bands centred on 5 and 7 GHz. Third panel: Evolution of the in-band spectral index α on 1-minute time scale. Here the radio flux density follows Sν ∝ να. Rising phases display flat-to-inverted spectra, while decay phases are generally optically thin. Bottom: Time-resolved and integrated upper limit on the linear polarisation fraction of the radio emission.

Fig. D.2. Cross-correlation function calculated with the ZDCF algorithm. We include the full range of time lags and a zoomed-in plot highlighting the peak lag. The solid vertical line and grey region show the best fit time lag and 1σ uncertainties, respectively. The vertical dashed line shows the zero-lag position, offset from the best-fit lag at a 1.8σ significance.

All Tables

All Figures

Fig. 1. Multi-band time series and evolution of optical and radio parameters (including data presented by Baglio et al. 2024). Error bars represent the 1σ confidence levels. In the seventh and eighth panels, horizontal dashed red lines mark the average optical polarisation degree (Popt) and angle (PAopt), with shaded areas indicating one standard deviation about the mean. In the ninth panel, a horizontal dotted line shows a spectral index of 0 (i.e., a flat spectrum). In the final panel, downward arrows denote 3σ upper limits on the radio polarisation degree (Pradio). TDB stands for barycentric dynamical time.

In the text

	Fig. B.1. Stokes parameters Qopt and Uopt for the optical emission. Data points are colour-coded according to the epoch of their measurements.
In the text

Fig. B.2. Polarisation degree (top), Stokes parameters (middle), and polarisation angle (bottom) as a function of the relative intensity for the optical emission from J1023. The top panel shows the results of the Spearman rank correlation test and of a linear regression algorithm to the data. The shape of the anti-correlation is consistent with ΔPopt / Δmag = -1.0±0.2, with an intrinsic random scatter around the best regression model of ${\sigma }_0=0.{07}_{-0.02}^{+0.04}$ (see text for details).

In the text

	Fig. C.1. Top: VLA 6 GHz light-curve of J1023 from the September 2024 observation, with 1-minute time binning. The short-time scale variability from the flaring activity is evident. Bottom: in-band spectral index with 1-minute time bins. The radio spectrum shows transitions between optically thick and thin phases on timescales of minutes.
In the text

Fig. D.1. Top: VLA 6 GHz light-curve of J1023 from the June 2024 observation during the flaring period, with 1-minute time binning. The radio emission shows multiple flares of similar duration, with green and yellow regions indicating flux rise and decay phases, respectively. Second panel: The same light-curve split into two sub-bands centred on 5 and 7 GHz. Third panel: Evolution of the in-band spectral index α on 1-minute time scale. Here the radio flux density follows Sν ∝ να. Rising phases display flat-to-inverted spectra, while decay phases are generally optically thin. Bottom: Time-resolved and integrated upper limit on the linear polarisation fraction of the radio emission.

In the text

	Fig. D.2. Cross-correlation function calculated with the ZDCF algorithm. We include the full range of time lags and a zoomed-in plot highlighting the peak lag. The solid vertical line and grey region show the best fit time lag and 1σ uncertainties, respectively. The vertical dashed line shows the zero-lag position, offset from the best-fit lag at a 1.8σ significance.
In the text