© The Authors 2023

1. Introduction

The evolutionary path that leads a neutron star (NS) in a low-mass X-ray binary (XRB) system to become a millisecond radio pulsar involves different stages, the most enigmatic of which is represented by transitional millisecond pulsars (tMSPs). These sources alternate between a radio pulsar state and an active low-luminosity X-ray disc state on timescales from weeks to months (for reviews, see Campana & Di Salvo 2018; Papitto & de Martin 2022). The archetypal tMSP is PSR J1023+0038 (hereafter referred to as J1023), which was discovered in 2007 as a 1.69 ms radio pulsar orbiting a low-mass (≃0.2 M_⊙) companion star with a period of 4.75 h (Archibald et al. 2009). In 2013, J1023 showed a sudden brightening in emission at X-ray and gamma-ray frequencies (by a factor of 5–10) as well as at ultraviolet (UV) and optical frequencies (by 1–2 mag), which coincided with the disappearance of the pulsed radio signal (Stappers et al. 2014; Patruno et al. 2014). Double-peaked optical emission lines were soon detected, indicating the formation of an accretion disc (Coti Zelati et al. 2014). J1023 has remained in this active state ever since, at an X-ray luminosity of L_X ≃ 7 × 10³³ erg s⁻¹ in the 0.3–80 keV energy range (Coti Zelati et al. 2018) for a distance of 1.37 kpc (Deller et al. 2012).

During the sub-luminous X-ray state, the X-ray emission of J1023 is observed to switch between two distinct intensity modes, dubbed high and low, with sporadic flares in between (e.g. Bogdanov et al. 2015). The high mode occurs 70–80% of the time, while the low mode occurs 20–30% of the time and typically lasts from a few tens of seconds to minutes. The drop and rise timescales involved in the mode switching are of the order of ≈10 s. X-ray, UV, and optical pulsations occur simultaneously in the X-ray high-intensity mode and disappear in the low mode (Archibald et al. 2015; Papitto et al. 2019; Miraval Zanon et al. 2022; Illiano et al. 2023). The UV, optical, and near-infrared (NIR) emission is dominated by the accretion disc and irradiated companion star and also shows flaring and flickering activities (Shahbaz et al. 2015, 2018; Kennedy et al. 2018; Papitto et al. 2018; Hakala & Kajava 2018). Optical polarimetric observations reveal a linear polarisation (LP) at a level of ∼1% (Baglio et al. 2016; Hakala & Kajava 2018) of uncertain origin. Bright, variable radio continuum emission is also detected (Deller et al. 2015). In particular, the low X-ray mode episodes are accompanied by an increase in the radio flux by a factor of ≃3 on average, with fairly equal rise and decay times and an overall duration that matches low X-ray mode intervals (Bogdanov et al. 2018; see also our Fig. 1). These re-brightenings are associated with a temporal evolution of the radio spectrum from inverted to steep. Currently, J1023 spins down at a rate that is ≃4% faster than in the radio pulsar state (Burtovoi et al. 2020). This overall behaviour is unique in the landscape of XRBs, and no comprehensive explanation for it has been found.

Fig. 1. Temporal evolution of the X-ray, UV, optical, NIR, and radio emissions of J1023 on the night of 2021 June 3–4. The light curves are shown in decreasing order of energy band from top to bottom, with the XMM-Newton/EPIC one as the reference for all others (NuSTAR, HST/STIS, XMM-Newton/OM, NTT/SOFI, and VLA). The XMM-Newton X-ray light curve is shown over the time interval in which the three EPIC instruments collected data simultaneously. The OM residual count rates were obtained after removing the orbital modulation from the time series (see the main text for details). For displaying purposes, only VLA detections are shown. The dark yellow segment in the bottom panel indicates the time range covered by the FAST observation. The orange shaded areas mark the time intervals of the low X-ray mode detected by XMM-Newton/EPIC covered by observations with other instruments.

This paper presents the results of the largest multi-wavelength campaign ever performed on J1023. Our aim is to ultimately understand the root cause of the mode changes in the active sub-luminous X-ray state and its implication for our understanding of the accretion-ejection coupling in XRB systems that harbour rapidly spinning NSs. The observations of our campaign were conducted on two different nights in June 2021. The first night (Night 1) involved XMM-Newton, the Nuclear Spectroscopic Telescope Array (NuSTAR), the Hubble Space Telescope (HST), the Son OF ISAAC (SOFI) instrument mounted on the New Technology Telescope (NTT), the Karl G. Jansky Very Large Array (VLA), and the Five-hundred-meter Aperture Spherical Telescope (FAST); the second night (Night 2) involved the Neutron Star Interior Composition Explorer Mission (NICER), the Neil Gehrels Swift Observatory (Swift), the Very Large Telescope (VLT), the Rapid Eye Mount (REM) telescope, and the Atacama Large millimetre/sub-millimetre Array (ALMA).

The paper is structured as follows: in Sect. 2 we present the observations and data analysis. The results are then reported in Sect. 3. In Sect. 4 we propose a scenario for the mode switching, present the results of the modelling of the spectral energy distributions (SEDs) in the two modes, and discuss the results in relation to our scenario. Conclusions follow in Sect. 5.

2. Observations and data analysis

Table 1 lists all the observations of J1023 presented in this paper. In the following, we describe these observations and the procedures adopted to process and analyse these data.

2.1. XMM-Newton (Night 1)

XMM-Newton (Jansen et al. 2001; Schartel et al. 2022) observed J1023 on 2021 June 3–4 (Obs. ID 0864010101) using the European Photon Imaging Cameras (EPIC) and the Optical/UV Monitor (OM) telescope (Mason et al. 2001). The EPIC-pn (Strüder et al. 2001) was set in fast timing mode (time resolution of 29.52 μs), the two EPIC-MOS (Metal Oxide Semi-conductor; Turner et al. 2001) in full frame mode (time resolution of 2.6 s) and the OM in the fast window mode (time resolution of 0.5 s) with the B filter (effective wavelength 4330 Å; full width at half maximum 926 Å). The observation data files were processed and analysed using the Science Analysis System (SAS; v. 20.0). The X-ray data were analysed in the same way as described by Miraval Zanon et al. (2022). No periods of high background activity were detected. The background-subtracted light curve was extracted over the time interval during which all three EPIC instruments collected data simultaneously and was binned with a time resolution of 10 s. The optical background-subtracted light curve was extracted using the OMFCHAIN pipeline with default parameters. The light curve was binned at a time resolution of 120 s to ensure a high S/N and detect any changes in the optical brightness at least during prolonged periods of low X-ray mode. The net count rates were then converted to magnitudes in the Vega system using the formula provided on the XMM-Newton online threads¹. To align the X-ray and optical light curves with those extracted from other observatories, the arrival times of photons were converted from the terrestrial time standard to the Coordinated Universal Time (UTC) standard without applying a barycentric correction².

2.2. NuSTAR (Night 1)

NuSTAR (Harrison et al. 2013) observed J1023 for a total of ∼48.6 ks and a net exposure of ∼22 ks (Obs. ID 30601005002), almost fully overlapping with the XMM-Newton observation. Data were processed and analysed using the NuSTAR Data Analysis Software (NUSTARDAS, v.2.1.2) with the instrumental calibration files stored in CALDB v20220525 and the same prescriptions reported by Miraval Zanon et al. (2022). J1023 was detected up to energies of ∼50 keV in both focal plane modules (FPMA and FPMB). Light curves were extracted separately for the two modules, combined to increase the S/N, and binned at a time resolution of 50 s.

2.3. HST (Night 1)

The Space Telescope Imaging Spectrograph (STIS; Woodgate et al. 1998) on board HST observed J1023 for 4750 s using the Near-UV Multi-Anode Micro-channel Array (NUV-MAMA) detector in TIME-TAG mode (time resolution of 125 μs); the data were collected with the G230L grating equipped with a 52″ × 0.2 slit with a spectral resolution of ∼500 over the nominal range (first order). The data analysis procedure that we adopted is the same as that described by Miraval Zanon et al. (2022). In summary, we used the stis_photons package³ to correct the position of the slit channels and to assign the wavelengths to each time of arrival. We selected times of arrival belonging to channels 993–1005 of the slit and in the 165–310 nm wavelength interval to isolate the source signal.

2.4. NTT/SOFI (Night 1)

We collected NIR high time-resolution imaging data using the SOFI (Moorwood et al. 1998) instrument mounted on the European Southern Observatory (ESO) 3.6 m NTT at La Silla Observatory (Chile). Observations were carried out using the J filter from 2021 June 3 at 23:49:09 UTC to June 4 at 02:42:38 UTC (Table 1). In order to achieve 1 *s time resolution, observations were performed in Fast Photometry mode, which restricts the area on which the detector is read. Specifically, the instrument was set to read a 400 × 400 pixels window (1.9′×1.9′). During the observation, a series of 100 frames each with an exposure time of 1.0 s were stacked together to form a ‘data cube’. There was a gap of ∼3 s between each cube. Photometric data were extracted through aperture photometry with an aperture radius of 1.4″ using a custom pipeline based on the DAOPHOT package (Stetson 1987). In order to minimise systematic errors, we performed photometry on J1023 and a reference star (at position RA = 10^h23^m50^s54, Dec = +00°38′164, J2000.0, and with a magnitude of J = 15.44  ±  0.05 mag in the Vega system from the Two Micron All-Sky Survey (2MASS) point source catalogue), assumed to have constant brightness. We then re-binned the light curve at a time resolution of 10 s to reduce the scatter and investigate the presence of variability at the same time resolution as the light curves in the X-ray and UV bands. For a more conservative approach, we evaluated the errors on the magnitudes of the re-binned light curves as the standard deviation of the magnitudes measured in each bin.

2.5. VLA (Night 1)

The VLA observed J1023 for 4 h, starting on 2021 June 3 at 22:30 UTC. The array was in its C-configuration and data were taken at a central frequency of 10 GHz, with a bandwidth of 4 GHz. J1331+3030 was used as flux and band-pass calibrator, while J1024−0052 was used as phase calibrator. The separation between the target and the phase calibrator is 0.3°. To calibrate and analyse the data we followed the prescription of Deller et al. (2015) and Bogdanov et al. (2018). Here we report the main steps. After initial calibration using the VLA pipeline within the Common Astronomy Software Application (CASA) package (version 5.1.1), manual flagging of bad data was performed. Unfortunately, the first half of the observation turned out to be unusable, possibly due to the bad weather conditions that occurred during the observation, and was completely flagged out.

Subsequently, we self-calibrated the data using a solution interval of 75 s for the phase calibration and 45 min for the amplitude+phase calibration. Upon inspection of the field of view (FOV), three main sources were identified: the target J1023, a faint source to the north, namely J102348.2+004017, and the brightest source in the field, located to the south-east, namely J102358.2+003826. J102358.2+003826 is a well-known galaxy at redshift z = 0.449 (Ahn et al. 2014) with a resolved extended emission whose secondary lobes heavily affect the overall root mean square (rms) noise level (Deller et al. 2015). To remove the problematic contribution of J102358.2+003826, we first cleaned the image using a mask with the tclean task. We set the nterms parameter equal to 2 to take into account the spectral properties of the galaxy over the wide bandwidth of the observation (Deller et al. 2015). This step provided the model visibilities of J102358.2+003826. Next, we used the uvsub task to subtract these model visibilities from the self-calibrated visibilities and deconvolved the residuals. The resulting image therefore does not have any contribution from J102358.2+003826.

As pointed out by Deller et al. (2015), the gain solutions from the self-calibration obtained before the subtraction of J102358.2+003826 were dominated by the pointing errors of this source. This source was located beyond the half-power point of the antenna, resulting in inaccurate solutions. Therefore, we had to invert the amplitude and phase self-calibration after subtracting this source from the field. We used the CASA task invgain (Hales 2016) to invert the calibration table and applied it to the self-calibrated visibilities. The final image obtained had a beamwidth size of 3.6″ × 2.1″, and an rms noise level of ≃4 μJy beam⁻¹.

To extract the light curve, we imaged the field using a 60 s time interval and a mask around the target with the CASA tclean task, setting niters=10 and cycleniter=10. We selected these values after evaluating random images where the target was either detected or not detected, and ensuring that this choice did not result in false positive detections. Then, we used the imfit task to obtain the peak flux density of the target and the check source, J102348.2+004017, for each 60 s image to investigate possible systematic effects. To measure the rms noise level around J1023 and check source, we placed four boxes around each source and used the imstat task to average the rms noise level within each box. To build the light curve, we considered a source as detected whenever the flux density value obtained with imfit was ≥3σ; otherwise, we took three times the rms noise level as an upper limit. To take into account possible systematics affecting the FOV, we investigated possible correlations between the light curves of the target and the check source. No evidence of such correlations was found. Finally, we further flagged bad 60 s images and removed intervals with odd results, such as sudden increases in the flux density of J102348.2+004017.

2.6. FAST (Night 1)

FAST (Jiang et al. 2019, 2020; Qian et al. 2020) observed J1023 starting on 2021 June 4 at 08:47:45 UTC. The total integration time was 3 h (the first two minutes and the last two minutes were dedicated to the injection of an electronic noise signal, referred to as a CAL). The observation was planned so as to cover a range of orbital phases from approximately 0.5 to 1.1, and hence to encompass the epoch of passage of the pulsar at the inferior conjunction of the orbit (corresponding to a phase of 0.75). This strategy should allow us to minimise the impact of eclipses of the radio signal caused by the gas out-flowing from the companion star due to irradiation from the pulsar wind. The observation sampled a total of 21 low mode episodes according to the simultaneous XMM-Newton observation (see Fig. 1). The central observing frequency was 1.25 GHz, with a range from 1.05 GHz to 1.45 GHz, including a 20 MHz band edge on each side. The average system temperature was 25 K. The recorded FAST data stream for pulsar observations is a time series of total power per frequency channel, stored in psrfits format (Hotan et al. 2004) from a ROACH-2-based back-end⁴, which produces 8-bit sampled data over 4k frequency channels at a cadence of 49.152 μs.

We searched for pulsations with either a dispersion signature or instrumental saturation in all data collected during the observation. We performed four different types of data processing.

Data folding based on the ephemeris derived from the X-ray data. We analysed the data using the ephemeris derived from the strictly simultaneous X-ray data collected by XMM-Newton (Miraval Zanon et al. 2022). We separately folded the FAST data from the entire observation and those obtained by stacking the data acquired during all the periods of high X-ray mode and all the periods of low X-ray mode. However, no radio pulses were detected with S/N > 6 in either case.

Dedicated search for persistent periodic radio pulses. For the sake of robustness, we created de-dispersed time series for each pseudo-pointing over a dispersion measure (DM) range between 3 pc cm⁻³ and 100 pc cm⁻³. This range covers all uncertainties of the semi-blind search, as J1023 was discovered at a DM of 14.3 pc cm⁻³ (Archibald et al. 2009). We set the step size between subsequent trial DMs equal to the intra-channel smearing over the entire band. This ensures that any trial DM deviating from the source DM by ΔDM does not cause extra smearing beyond the intra-channel smearing⁵. We performed a ‘jerk’ acceleration search (Andersen & Ransom 2018; see also Wang et al. 2021) for periodic signals in the de-reddened power spectrum of each de-dispersed, barycentred time series. We restricted the search to various frequency ranges around the nominal spin frequency of J1023, and allowed the power of the highest harmonic component of the signal to drift in the Fourier domain up to a maximum number of 200 frequency bins and 50 frequency derivative bins. The powers of the first eight harmonics were summed (in powers-of-two harmonics). All identified candidate periodicities underwent a sifting procedure to reject less significant, duplicated (i.e. detected at different DMs and accelerations), and/or harmonically-related candidates. Using the prepfold tool, the data were folded at the nominal period, period derivative and DM of the sifted candidates. Corresponding diagnostic plots including the integrated pulse profile, the signal intensity as a function of time and frequency, and the DM versus χ² values were generated and saved for visual inspection. The candidates were judged based on persistent and broadband emission, as well as distinct peaking of the signal’s significance at specific DM values. No signal with DM, period, and period derivative close to the expected values was detected.

Search for single pulses. We used the dedicated search scheme in (ii) to de-disperse the time series. We then applied 14 box-car width-match filter grids distributed in logarithmic space from 0.1 to 30 ms. A zero-DM matched filter was used to mitigate radio frequency interference (RFI) in the blind search. All the resulting candidate plots were visually inspected, but most candidates were determined to be RFI. No pulsed radio emission with a dispersive signature was detected with S/N > 5.

Search for baseline saturation. We acknowledge that FAST would become saturated if the radio flux density reached levels as high as 1 kJy to 1 MJy. Therefore, we also searched for saturation signals in the recorded dataset. We looked for epochs where at least half of the channels satisfied one of the following conditions: (1) the channel was saturated (255 value in 8-bit channels), (2) the channel had a zero value, or (3) the rms deviation of the bandpass was < 2. We did not find any instances of saturation lasting > 0.5 s during the observational campaign.

We then calibrated the noise level of the baseline with noise CAL injection, and derived a 5-σ flux density upper limit for persistent radio pulsations of 2.0 ± 0.5 μJy at 1.25 GHz considering the entire integration time of 2.8 h, assuming a pulse duty cycle between 0.05–0.3 (which is typical for millisecond pulsars; see the ATNF pulsar catalogue; Manchester et al. 2005) and disregarding the possible scintillation effect. Within the uncertainty of the flux measurements and stability threshold, no flux variations were seen during the observation, consistent with an apparent lack of scintillation. The 5-σ flux upper limits on the pulsed flux density of J1023 during the high and low modes can be evaluated by scaling the above-mentioned limit for the specific integration time spans during the high and low-mode periods in the X-rays: we obtain upper limits of 2.2 ± 0.6 μJy in the high mode (2.3 h) and 4.7 ± 1.0 μJy in the low mode (0.5 h).

2.7. NICER/XTI observations (Night 2)

The X-ray Timing Instrument (XTI) mounted on NICER (Gendreau et al. 2012) observed J1023 between 2021 June 26 at 21:56:04 UTC and June 27 at 01:33:51 UTC, for an exposure time of ∼6.1 ks. The raw, unfiltered data were processed and analysed using tools within the NICERDAS software package and the latest version of the calibration files at the time of the analysis (CALDB v. 20221001)⁶. We extracted the background-subtracted light curve in the 0.5–10 keV energy range with a time bin of 10 s using the NICERL3-LC pipeline. This pipeline utilises the ‘Space Weather’ background model, that is to say, it estimates the background level on the basis of the local spacecraft environment (see Fig. 2, first panel).

Fig. 2. Zoom-in of the time series collected on Night 1 over the time span of maximum overlap among the different telescopes. Arrows indicate 3σ upper limits.

2.8. Swift/XRT observations (Nights 1 and 2)

The Swift (Gehrels et al. 2004) X-ray Telescope (XRT; Burrows et al. 2005) observed J1023 in the photon-counting mode (time resolution of 2.5 s) on both Night 1 and Night 2. The first observation started on 2021 June 3 at 18:29:46 UTC and consisted of 11 snapshots, with a total elapsed time of ≃51.3 ks and an effective exposure time of ∼17.9 ks. The second observation started on 2021 June 26 at 22:19:35 UTC and consisted of three snapshots, resulting in a total elapsed time of ≃7.3 ks and an effective exposure time of ∼2.4 ks. The raw, level-1 data were processed and screened using standard criteria. Source photons were collected from a circular aperture of radius 47.2″ centred on J1023 and background photons were collected from an annular region with inner and outer radii of 94.4″ and 188.8″, respectively, also centred on J1023. We extracted a background-subtracted time series with a time bin of 100 s, suitable to sample adequately the X-ray mode switching (see Figs. 1 and 3).

Fig. 3. Temporal evolution of the X-ray, UV, optical, and millimetre emissions and of the optical LP of J1023 on the night of 2021 June 26–27. The light curves are shown in decreasing order of energy band from top to bottom (NICER, Swift/XRT, Swift/UVOT, REM, VLT/FORS2 polarisation, and ALMA). In the bottom panel, the ALMA time series are plotted using either bins of length 4–5 min (cyan) or bins of variable length corresponding to the duration of the X-ray mode episodes identified by the strictly simultaneous NICER and Swift data (blue). In some panels, the sizes of the markers are bigger than the uncertainties and/or the time bins. Arrows indicate 3σ upper limits. The orange shaded areas mark the time intervals of low X-ray mode detected by NICER and Swift/XRT covered by observations with other instruments. A plot of the X-ray and millimetre light curves around the epoch of the millimetre flares is shown in Fig. 6.

2.9. Swift/UVOT observations (Nights 1 and 2)

The Swift Ultra-Violet/Optical Telescope (UVOT; Roming et al. 2005) observed J1023 simultaneously with the XRT. The observations were performed in both image and event modes using the UVM2 filter (central wavelength 2260 Å; full width at half maximum 527 Å). Data were processed and analysed in the standard way. For the source photons, we used a circular extraction region of radius of 5″ centred on the source position. For the background photons, we adopted a circular extraction region of radius 20″ away from source. For the data in image mode, photometry was carried out using the uvotsource task. For the data in event mode (time resolution of 11 ms), we used the coordinator tool to convert raw coordinates to detector and sky coordinates, uvotscreen to filter out hot charge-coupled device (CCD) pixels and obtain a cleaned event list, and uvotevtlc to extract a background-subtracted light curve with a time bin of 120 s. This time bin provided a sufficiently high S/N to detect any changes in the UV brightness at least during the longest periods of low X-ray mode.

2.10. REM photometric observations (Night 2)

We observed J1023 on Night 2 with the 60-cm REM telescope (Zerbi et al. 2001; Covino et al. 2004) located at La Silla Observatory (Chile). We obtained 48 strictly simultaneous exposures with the ROSS2 imager, each with an integration time of 60 s, with Sloan Digital Sky Survey (SDSS) griz optical filters (Table 1). We reduced the images following standard procedures (bias subtraction and division by a normalised flat-field), and performed aperture photometry with the PHOT tool in IRAF using an aperture size of 6 pixels. We performed flux calibration using stars in the field from the American Association of Variable Star Observers (AAVSO) Photometric All-Sky Survey (APASS) catalogue⁷ (Henden 2019). We obtained simultaneous K-band observations with the REM infrared (REMIR) camera. A total of 300 15 s integration images were acquired, which we combined 5 by 5 to correct for the background contribution. Aperture photometry was performed with IRAF as described above. The images were flux-calibrated against a selection of stars in the field from the 2MASS point source catalogue⁸.

2.11. VLT/FORS2 polarimetric observations (Night 2)

We observed J1023 with the FOcal Reducer/low dispersion Spectrograph 2 (FORS2; Appenzeller et al. 1998) mounted on the 8.2-m VLT at Paranal Observatory (Chile) in polarimetric mode using the optical filter R_SPECIAL + 76 (R; central wavelength 655 nm, full width at half maximum 165 nm). The observations were performed on Night 2 under thin cirrus clouds conditions, with seeing varying from 0.4″ to 0.9″. A total of 156 short (10 s) exposures were taken, resulting in an observing time of ∼2.41 h, overheads included (about half the orbital period).

A Wollaston prism was inserted in the optical light path to split the incident radiation into two orthogonally polarised light beams called ordinary (o) and extraordinary (e) beams. A Wollaston mask was introduced to avoid overlap of the beams on the CCD. In addition, a rotating half-wave plate (HWP) was inserted, allowing images to be taken at four different angles (Φ_i) of the plate with respect to the telescope axis: Φ_i = 22.5° (i − 1), i = 1, 2, 3, 4. This step is essential to obtain a polarisation measurement, since the set of images at the four different angles must be combined to evaluate the level of LP (as described below). Therefore, a total of 39 polarisation epochs were acquired with the programme. All images were processed by subtracting an average bias frame and dividing by a normalised flat field. Aperture photometry was performed using the daophot tool (Stetson 1987), using an aperture of 6 pixels. The normalised Stokes parameters Q and U for the LP were evaluated following the methods described by Baglio et al. (2020) and references therein. These values were then corrected for instrumental polarisation by observing a non-polarised standard star⁹ (WD 1620-391) using the same setup on the same night. A polarisation level of < 0.2% was measured, which is consistent with the typical instrumental polarisation level expected for the FORS2 instrument as reported on the ESO website¹⁰ (< 0.3% in all optical bands). The Stokes parameters evaluated for the non-polarised standard star were then subtracted from the Stokes parameters of J1023 at each epoch.

We then used an algorithm to estimate the level and angle of LP. This algorithm is based on the evaluation of the quantity S(Φ) for each of the HWP angles (e.g. di Serego Alighieri 1998 and references therein; see also Covino et al. 1999; Baglio et al. 2020). S(Φ) is the component of the Stokes vector that describes the LP along the Φ direction (e.g. S(0°) = Q). It is defined as

(1)

where Φ is the angle of the HWP, f^o(Φ) and f^e(Φ) are the ordinary and extraordinary fluxes for a given angle Φ, and and are the same quantities evaluated for the non-polarised standard star. S(Φ) is related to the polarisation level P and angle θ by S(Φ) = P cos2(θ − Φ). The parameters P and θ were estimated by maximising the Gaussian likelihood function and using this as the starting point for a Markov chain Monte Carlo (MCMC) procedure (Hogg & Foreman-Mackey 2018; see Baglio et al. 2020 for further details on the algorithm). The best fit parameters were then calculated as the median of the marginal posterior distributions of these parameters, while the uncertainties were estimated as the 16–84th percentiles of the distributions. For non-detections, the 99.97% percentile of the posterior distribution of the parameter P was used to calculate an upper limit.

After fitting the S(Φ) curve, we obtained an LP value that is already corrected for instrumental effects. However, the polarisation angle still needs to be corrected using the Stokes parameters of a polarised standard star. For this purpose, we observed the star Hiltner 652 and measured its Stokes parameters to determine its polarisation angle. This was then compared to the tabulated value, and a difference of 1.71° ± 0.30° was applied to the polarisation angle of J1023. However, the resulting LP will still include an interstellar component. Evaluating the interstellar contribution is not straightforward. In principle, we could use a group of reference stars in the same field as the target in Eq. (1), under the assumption that these stars are non-polarised. This way, the LP of the target would be corrected for both instrumental and interstellar contributions. However, it is important to ensure that the chosen reference stars all have the same level of polarisation, which would be a good indication of their intrinsic non-polarised nature (see e.g. Baglio et al. 2020). Unfortunately, not only is the FOV of J1023 sparse in the optical, but the very short exposure times (10 s) exclude all faint (R > 19) stars from this study. We are therefore left with only one possible comparison star in the field, which makes it difficult to verify whether the chosen star is indeed intrinsically non-polarised. As a result, we used the fluxes of the non-polarised standard star in Eq. (1). The field star will be used as a comparison for the analysis.

We can evaluate the maximum interstellar contribution to the LP of a source, referred to as P_int, through considerations of its absorption coefficient in the V band, A_V: P_int(%) < 3A_V (Serkowski et al. 1975). For J1023, the best-fitting absorption column density, N_H, according to XMM-Newton observations is cm⁻² (Campana et al. 2019). Using the relation between N_H and A_V in our Galaxy from Foight et al. (2016), we calculate A_V = 0.17 ± 0.01 mag and therefore P_int < 0.52%. This is very low and comparable to the instrumental polarisation measured for the instrument configuration. It is also comparable to the upper limit of the interstellar polarisation measured for J1023 (Baglio et al. 2016).

2.12. ALMA observations (Night 2)

The ALMA observations were conducted using the 12 m array and 38 antennas in the C43-7 configuration (baselines 22–2492 m), with a FOV of 59.72″ and a maximum recoverable scale of 11.5″. The total integration time on source was ∼2 h and the phase centre was RA = 10^h23^m47^s68; Dec = +00°38′4100 (ICRS). The weather conditions were good and stable, with an average precipitable water vapour of ∼1.2 mm and an average system temperature of ≃65 K. The band-pass and flux calibrators used were J1058+0133 and J1256−0547, while J1028−0236 was used as the phase calibrator. The continuum image was obtained by averaging the four base-bands of width 1.875 GHz, each centred at 90.5003 GHz, 92.5003 GHz, 102.5003 GHz, and 104.5003 GHz, covering a total bandwidth of 7.5 GHz. No spectral lines were detected within the observed frequency range above the 3σ noise level per spectral channel of ∼0.48 mJy beam⁻¹.

Data reduction was carried out following standard procedures using the ALMA pipeline within CASA (v. 6.1.1-15)¹¹. Continuum images were produced using tclean interactively and applying a manually selected mask. These images were corrected for primary beam attenuation and were produced using both natural weighting and a Briggs robust=0 weighting (Briggs 1995). The former maximises the sensitivity whereas the latter provides higher angular resolution. The mean frequency of the continuum image is 97.51 GHz, which corresponds to a wavelength of 3.1 mm. Table 2 lists the synthesised beam and the rms noise of each image.

3. Results

3.1. X-ray and UV emissions

The X-ray light curve extracted from XMM-Newton data on the first night reveals flaring activity over a time span of ∼1.5 h at the beginning of the observation. A total of 288 switches between high and low modes were observed during the observation (see Fig. 1). Based on our thresholds for the EPIC net count rates in the different modes (see Miraval Zanon et al. 2022), we estimated that J1023 spent ≃65, 17 and 1% of the time in the high, low and flaring modes, respectively, and the rest of the time switching between modes. The good time intervals corresponding to the high and low modes were then used to extract the background-subtracted spectra, response matrices, and ancillary files separately for each mode. The NuSTAR light curve is characterised by a lower S/N than the XMM-Newton/EPIC light curve. Nevertheless, it is possible to recognise the same flaring episodes and mode switchings detected using XMM-Newton (see Fig. 1). We used the good time intervals for the high and low modes selected by the strictly simultaneous XMM-Newton/EPIC light curves to extract the background-subtracted spectra, response matrices, and ancillary files separately for each mode.

The X-ray light curve extracted from NICER data on the second night, shown in the top panel of Fig. 3, reveals 22 instances of X-ray mode switching, without any obvious flaring episodes. We applied intensity thresholds to the light curve to select the time intervals corresponding to the high (> 2.5 counts s⁻¹) and low (< 2.5 counts s⁻¹) X-ray modes. Based on this selection, we estimate that J1023 spent ∼70% of the time in the high mode and ∼30% of the time in the low mode. We also detected two periods of X-ray low mode in the Swift/XRT data, one of which was also covered by the NICER observations (see the second panel of Fig. 3).

The UV light curves exhibit substantial variability. Specifically, the HST light curve on Night 1 displays flaring activity as well as dips in the count rate that coincide with the low-mode episodes observed in X-rays (Fig. 1; see also Miraval Zanon et al. 2022; Jaodand et al. 2021). A dip in the UV brightness by ∼0.7 mag is also clearly visible in the Swift/UVOT data during the third snapshot of Night 2. The decrease occurs over a span of several minutes during an episode of low mode (Fig. 3).

3.2. Optical/NIR photometric properties

On the first night, the source is brighter at the beginning than in the rest of the observation (reaching B ≃ 17 − 16.5 mag), displaying a smooth quasi-sinusoidal modulation at the orbital period with a semi-amplitude of ≃0.2 mag around an average value of B ≃ 17.3 mag. We fit a sinusoidal function to the light curve by fixing the period and reference phase of the sinusoid to the binary orbital period and the phase value obtained from the time of passage of the pulsar through the ascending node, respectively, as derived by Miraval Zanon et al. (2022). Then, we subtracted the best-fitting function from the observed count rates. The resulting detrended light curve is shown in Fig. 1. After initial flaring activity resembling the X-ray behaviour, the optical intensity appears to be flickering, but without any sign of mode switching.

Figure 1 also shows the NIR light curve. Variability in the form of dips and flickering is evident. The fractional rms magnitude deviation of the light curve, computed following the method described by Vaughan et al. (2003), is (10.1 ± 0.2)%. No NIR flare was observed at the switch from the low to the high X-ray mode, as previously reported (Papitto et al. 2019; Baglio et al. 2019). Furthermore, we did not detect any signs of mode switching in the light curve. The average magnitudes measured separately during the high and low modes are compatible within 1σ (J = 15.39 ± 0.16 mag and J = 15.37 ± 0.13 mag in the high and low modes, respectively), indicating that there is no significant correlation between the NIR variability and mode switching.

On Night 2, the optical light curves show clear variability on short timescales (Fig. 3, fourth panel). Unlike what was seen on Night 1, no modulation at the system orbital period is detected. The average magnitudes are g = 16.65 ± 0.02 mag, r = 16.45  ±  0.02 mag, i = 16.20  ±  0.03 mag, z = 16.09  ±  0.03 mag (in the AB photometric system), all marginally consistent with those reported by Coti Zelati et al. (2014), Baglio et al. (2016). To estimate the amount of intrinsic variability, we evaluated the fractional rms magnitude deviation of the light curves following Vaughan et al. (1994). We find a low-significance fractional rms of 6.2 ± 3.4%, 6.9 ± 2.3% and 12.7 ± 4.7% at a time resolution of 160 s in the g, r and z bands, respectively, while no excess intrinsic variability could be detected in the i band. This fractional rms is comparable to what has been observed for a few black hole XRBs in the optical band on a similar timescale, such as GX 339-4 (Gandhi 2009; Gandhi et al. 2010) and MAXI J1535−571 (Baglio et al. 2018). In these systems, the variability has been attributed to the emission of a flickering jet that contributes up to optical frequencies. In contrast, when accretion dominates, lower values of fractional rms are typically measured, as seen in the case of the black hole XRB GRS 1716−249 (Saikia et al. 2022).

Unfortunately, most of the K-band images are of bad quality, resulting in an empty FOV. We stacked together three of the combined images in which the field was visible and J1023 was detected (so a total of 15 × 15 s integration images, all taken during the high mode), and performed the calibration against a selection of stars in the field from the 2MASS point source catalogue. We obtained an average magnitude of K = 14.41 ± 0.29 mag, which is consistent with previous values reported in the literature (see e.g. Coti Zelati et al. 2014; Baglio et al. 2016).

3.3. Optical polarimetric properties

Figure 3 shows the evolution of the optical LP level on the second night of observations. To improve the S/N of the observations without losing too much information on the short-term variability of the LP, we combined the measurements so that each LP data point covers 80 s (four HWP angles ×10 s of integration × two sets). We define an LP measurement as a detection when it reaches a confidence level of greater than 2σ. Using this threshold, we obtain a total of 8 LP detections and 11 upper limits at a 99.97% confidence level (see Table 3). The average level of LP, calculated as the average of all 2σ-significant detections, is 0.92 ± 0.39%. This is consistent with previous results in the same band (R band; Baglio et al. 2016). The polarisation angle is well constrained at a value of (corrected for the polarisation angle of the standard star).

When detected, the LP appears fairly constant throughout the observations. We investigated any potential changes in the LP level as a function of the X-ray mode (low/high) by averaging all images taken during the high and low modes. We found that the LP was 0.65 ± 0.43% (with a 99.97% confidence level upper limit of 2.12%) in the high mode and (with a 99.97% confidence level upper limit of 1.58%) in the low mode.

The measured levels of LP are not corrected for interstellar contribution, which is expected to be low according to Serkowski’s law (< 0.52%; Serkowski et al. 1975). Figure 4 shows the Q and U Stokes parameters on the Q − U plane for the whole dataset, including J1023 (black dots), a comparison star in the field (red squares), and the non-polarised standard star (green star). The non-polarised standard star is located very close to the origin of the Q − U plane, indicating that the instrumental polarisation of VLT/FORS2 is low, as expected. The Stokes parameters of the comparison star, on the other hand, are scattered around a common point, consistent with the position of the non-polarised standard within ∼1.5δ, where δ is the standard deviation of the points. Any deviation from the position of the non-polarised standard on the plane is a good estimate of the interstellar polarisation of the field (assuming that the comparison star is intrinsically non-polarised). Interestingly, the Stokes parameters of J1023 are clustered around a common point in a different region of the Q − U plane, with a higher scatter than the comparison star due to larger uncertainties. For the two time series of the Q and U Stokes parameters, we are unable to claim any detection of short-timescale variability.

Fig. 4. Stokes Q and U measurements for the whole VLT/FORS2 dataset. Black dots, red squares, and the green star indicate the Stokes parameters of J1023, a comparison star in the field, and the non-polarised standard star, respectively.

3.4. Millimetre emission

Figure 5 shows the ALMA 3.1 mm continuum emission towards J1023. Two sources were detected within the FOV of ALMA. J1023 was detected at a significance level of ∼15σ. The other source, which we designate ALMA J102349.1+003824 following the naming convention of new sources discovered by ALMA, is located ∼22″ south-east of J1023 and was detected at a significance level of ∼9σ. The panels on the right of the figure display the zoom-in maps around each source, with the continuum image obtained with natural weighting displayed in colour scale and the black contours derived from the image with robust=0. Table 4 reports the peak coordinates, the peak intensity, the integrated flux, and the deconvolved size obtained by fitting a 2D Gaussian function to the high-resolution image (robust=0).

Fig. 5. ALMA images of the field around J1023. Left: ALMA 3.1 mm continuum image showing the two detected sources, J1023 and ALMA J102349.1+003824. The white dotted boxes display the FOV shown in the right panels. The ALMA synthesised beam (0.87″ × 0.53″, PA = −74.5°) is shown at the bottom-right corner. Right: Close-up images of J1023 (top) and ALMA J102349.1+003824 (bottom). In both panels, the colour scale corresponds to the image with natural weighting, with the 3σ contour indicated in white (σ values are listed in Table 2). The black contours are derived from the ALMA image obtained with robust=0. Contours range from 3σ to 9σ in steps of 1σ, where σ = 9.1 μJy beam−1 for J1023 and σ = 15 μJy beam−1 for ALMA J102349.1+003824. The white cross depicts the radio position of J1023 (Deller et al. 2012). A faint emission tail is detected at the 3σ level extending south-west of J1023. The synthesised beams for the robust=0 and natural images are displayed in the bottom-left and bottom-right corners, respectively (see Table 2). In all panels, the colour scale is in units of Jy beam−1.

The 3.1 mm continuum emission from J1023 has a compact shape with some elongation towards the south-east (shown in the top right panel of Fig. 5 using black contours). The millimetre continuum emission peaks ∼0.32″ south of the optical position. The deconvolved size obtained from the 2D Gaussian fit in the robust=0 image is 0.305″ × 0.205″, which corresponds to approximately 418AU × 281AU (∼6.3 × 10¹⁵ cm × 4.2 × 10¹⁵ cm) at a distance of 1.37 kpc (this is approximately 7 × 10⁴ times the orbital separation of the system; Archibald et al. 2013). The emission from ALMA J102349.1+003824 also appears compact (see the bottom-right panel of Fig. 5). The 2D Gaussian fit yields a point-like source with a full with at half maximum as large as 0.39″ × 0.22″ (see Table 4). The results of the searches for counterparts of this source at other wavelengths are reported in Appendix A.

To determine the spectral index of J1023 within the ALMA bandwidth, we imaged the continuum emission from the lower and upper sidebands, which have central frequencies of 91.47886 GHz and 103.5115 GHz, respectively. The resulting spectral index is −0.1 ± 0.8, which is consistent with previous reports for the average cm emission (Deller et al. 2015; Bogdanov et al. 2018).

To probe variability in the millimetre emission, we obtained images scan by scan (typically with a duration of ∼4–5 min), using a natural weighting. We then computed the peak intensities and flux densities considering the emission above 3σ level. The bottom panel of Fig. 3 shows the evolution of the flux density over time. This quantity was found to vary by a factor of ∼5 within the range ∼30–160 μJy on a timescale of about 5 min. In some cases, the flux density remained below the detection sensitivity, resulting in upper limits of the order of a few tens of microjanskys. Figure 3 also shows the ALMA time series extracted using time intervals of variable length corresponding to the durations of the X-ray-mode episodes. Two flaring episodes can be clearly seen. We measured average flux densities of 89  ±  13 μJy in the high mode and 134 ± 23 μJy in the low mode, excluding the two flaring episodes (see Table 5). Overall, the anti-correlated variability pattern appears to be more pronounced between the cm and X-ray emission than between the millimetre and X-ray emissions.

The integration times adopted above are effective for studying the variability of the millimetre emission on minute timescales, but may average out large enhancements of the millimetre flux on much shorter timescales (of the order of seconds). Therefore, we performed a dedicated search for short flaring episodes by extracting images of different durations using the entire dataset. This search resulted in the detection of the two above-mentioned episodes of bright millimetre emission (Fig. 6). No other remarkable events were detected. The first episode (hereafter ‘Flare 1’) occurred on 2021 June 26 within the time range 22:19:29–22:19:44 (UTC) at the switch from high mode to low mode and reached a peak flux density of 1.01 ± 0.08 mJy. The second event (‘Flare 2’) took place during the time range 23:48:29–23:48:44 (UTC) and reached a peak flux density of 0.7 ± 0.1 mJy. The lack of continuous coverage in the millimetre band before and during the switch from the high to the low mode does not allow us to determine the exact epoch of the flare onset. We note that these values refer to the flux at the peak of the flares and are therefore higher than those obtained in the X-ray mode-resolved analysis (see Fig. 3 and Table 6), where the emission is averaged over longer time intervals. We observed a significant increase in the peak intensity above the average level: a factor of ∼4 for Flare 1 and a factor of ∼7 for Flare 2.

Fig. 6. X-ray (top) and millimetre (bottom) time series collected around the epochs when the two millimetre flares were detected. The X-ray light curve is binned at a time resolution of 10 s, and the millimetre light curve is binned at 15 s. Arrows indicate 3σ upper limits. Note the lack of coverage in the millimetre band for Flare 2 before and during the switch from the high to the low mode in the time interval between two consecutive ALMA scans. Note also the different scales on the horizontal axis in the panels on the left and on the right.

We observed a clear change in the morphology of the millimetre continuum emission during the two flaring episodes (see the white contours in Fig. 7). The 2D Gaussian fitting (see Table 6) indicates that the emission is elongated along the north-east to south-west direction during Flare 1 (roughly perpendicular to the ALMA 3.1 mm emission obtained with all the visibilities), while it is marginally resolved in only one direction during Flare 2. The nominal rms absolute positional accuracy for standard ALMA observations can be estimated as beam_FWHM × (0.9 × S/N)⁻¹¹², where beam_FWHM is the full width at half maximum synthesised beam in arcsec and S/N is the signal-to-noise ratio of the image target. The positional accuracy for Flare 1 is ≃0.12 − 0.15″, while for Flare 2 it is slightly smaller at ≃0.1″. However, the actual absolute positional accuracy may be worse than the nominal values, potentially by a factor of 2 or more, depending on the atmospheric phase conditions during the observation. The offset between the ALMA peak positions and the radio position (Deller et al. 2012) is ≃0.4″ for Flare 1 and ≃0.3″ for Flare 2. Given the stable weather conditions during the observation, it is possible that this offset is genuine for both flares. However, caution must be exercised, as this offset is comparable to the positional accuracy of ALMA.

Fig. 7. Contour levels for the millimetre emission detected during 15 s time intervals at the peak of Flare 1 (left panel) and Flare 2 (right panel) shown with solid white lines overlaid on the naturally weighted ALMA 3.1 mm image extracted using the whole integration time (colour scale). These contour levels start at 3σ and increase in steps of 1σ, where σ is the rms of the image, ∼90 μJy beam−1. Both panels display the same colour scale range in units of Jy beam−1. Time intervals (in UTC) are displayed at the top of each panel. The red cross depicts the radio position of J1023 (Deller et al. 2012). The synthesised beams, with extension 0.696″ × 0.540″ and position angle PA = −79.92° for Flare 1 and 0.931″ × 0.525″ and PA = −74.21° for Flare 2, are shown at the bottom-right corner of each panel.

3.5. Radio emission

3.5.1. Continuum emission

The radio light curve extracted from VLA data is shown in Fig. 1. J1023 is not detected most of the time, especially during the high mode intervals. However, it is detected during almost all episodes of low mode (here the uncertainties on the flux densities are reported at the 1σ confidence level). The average radio flux densities measured separately in the high and low mode are reported in Table 7. These values are consistent with those reported by Bogdanov et al. (2018) within the uncertainties. This indicates that the radio flux densities measured separately in the two modes are remarkably stable over timescales of years.

3.5.2. Non-detection of pulsed emission

The limit on the radio pulsed flux density derived from the FAST observations is a factor of > 5000 smaller than the flux density measured at similar frequencies when J1023 was active as a radio pulsar (Archibald et al. 2009) and a factor of a few tens more constraining than the limits previously evaluated (again at similar frequencies) as soon as J1023 transitioned to the sub-luminous X-ray state (Stappers et al. 2014; Patruno et al. 2014). The limit on the radio pulsed flux density in the high mode is also at least four orders of magnitude smaller than the value predicted by extrapolating the best-fitting power-law model for the SEDs of the optical-to-X-ray pulsed emission (Papitto et al. 2019; Miraval Zanon et al. 2022).

3.6. The low-mode ingress and egress timescales

We first examined the low mode ingress and egress times in the XMM-Newton observation. For each low-mode episode, we used a simple model for the count rate light curve consisting of four components: a constant level for the high mode (allowing for a different count rate at the entrance and at the exit of the low-mode episode), a linear ingress in the low mode, a constant count rate in the low mode, and a linear egress. This model is similar in some aspects to what is used to model eclipses in binary stars or planetary transits. While we do not expect to accurately describe each low-mode episode individually, we consider the results in a statistical sense.

We applied this model to each low-mode episode in the XMM-Newton/EPIC count rate light curve, which was binned at 10 s, excluding the initial flaring mode of the observation. We performed the same procedure on the NICER data from the second night, using again a 10 s time bin. We fit the model to 56 low-mode episodes for XMM-Newton and 10 for NICER.

As shown in Fig. 8, the ingress time is shorter than the egress time from the low mode (∼80%). Interestingly, we observe that the two low-mode instances detected in the NICER data and coinciding with short-duration flares detected by ALMA both have comparable ingress and egress times (see the left panel of Fig. 8). This is uncommon and may be related to the detectability of flares at millimetre wavelengths. The distributions of XMM-Newton low-mode ingress and egress times are shown in the right panel of Fig. 8. The average logarithmic ingress time is (10.4  ±  6.9) s, which is consistent with the time binning. The average logarithmic egress time is (30.6  ±  8.2) s. We also investigated whether there was any relationship between ingress time, egress time, duration of the low mode, count rate at the beginning and end of the low mode, and count rate in the low mode. However, we did not find any statistically significant correlations.

Fig. 8. Results of the statistical analysis of the X-ray mode switching timescales. Left: Duration of the low-to-high mode switch (egress) as a function of the duration of the high-to-low mode switch (ingress) computed using data from XMM-Newton/EPIC (in red) and NICER (green). The circled green points represent the values for the ingress associated with the two millimetre flares detected with ALMA and the subsequent egress. The solid black line indicates where the two quantities are equal. Right: Distribution of the low-mode ingress (in green) and egress (in red) times in logarithmic scale from XMM-Newton and NICER data.

Although the model we used to describe the low-mode episodes is simple, it is clear that the duration of the ingress time in the low mode is short and consistent with the adopted time binning, and the egress time from the low mode is measurable and longer.

4. Discussion

4.1. The physical picture

Our campaign, along with previous findings (see e.g. Ambrosino et al. 2017; Bogdanov et al. 2018; Papitto et al. 2019; Miraval Zanon et al. 2022; Illiano et al. 2023 and references therein), allows us to draw a comprehensive physical picture of the high-low mode switches in J1023, which involves a permanently active rotation-powered pulsar, an accretion disc and episodes of discrete mass ejection on top of a compact jet. In the following, we outline the physical scenario while in Sect. 4.2 we describe in detail the SED modelling in support of our scenario. Figure 9 shows a visual representation of our scenario. We do not discuss the flaring activity here as it was not covered over a wide energy range.

Fig. 9. Schematic visual representation of the evolution of the inner flow into an outflow at the switch from the high to the low X-ray emission mode. When the system is in the high mode (top panel), a small-size inner flow is present, together with a faint steady jet that is launched along the pulsar rotational axis and gives rise to the observed low-level radio and millimetre emission. As the pulsar rotates, the pulsar wind (marked with solid green lines) wobbles around the equatorial plane (see e.g. Bogovalov 1999) and shocks off the electrons in the inner flow at two opposite sides (red spots) at a distance that is slightly larger than the light cylinder radius (≃80 km). At each pulsar rotation, synchrotron emission at the shock at X-ray, UV, and optical frequencies is modulated at the spin period at one side (bright red spot), while it is absorbed by material in the inner flow at the other side (light red spot; Papitto et al. 2019). When the system makes a switch to the low mode (bottom panel), a bright discrete ejection is launched along the pulsar rotational axis on top of the compact jet, the inner flow disappears, and the shock emission is quenched. For displaying purposes, only the inner regions of the compact jet are represented in both panels.

We assume that during the high mode the accretion flow is kept just beyond the light cylinder radius (≃80 km) by the radiation pressure of the particle wind from an active rotation-powered pulsar (Papitto et al. 2019). This causes the innermost region of the truncated, geometrically thin accretion disc to be replaced by a radiatively inefficient, geometrically thick flow, as expected for low-accretion rates (Narayan & Yi 1994). In this configuration, synchrotron radiation at a shock front that forms where the pulsar wind interacts with the inner flow produces the bulk of the X-ray emission as well as X-ray, UV, and optical pulsed emission (Papitto et al. 2019; Veledina et al. 2019). The collision of high-energy charged particles from the pulsar wind with the in-flowing matter in the disc increases the ionisation level of the disc. As a result, the in-flowing matter is drawn into the pulsar magnetic field and accelerated, producing a compact jet of plasma that streams out. According to our SED modelling, the corresponding jet synchrotron spectrum is optically thick up to the mid-infrared band and breaks to optically thin above a frequency of Hz, giving only a small contribution at X-ray energies (∼8–18% in the 0.3–10 keV band; see Fig. 10 and Table 8).

Fig. 10. Broadband (radio to X-rays) SEDs of J1023 during the high and low modes, along with the best-fitting models. Both SEDs were corrected for extinction effects. The red-dashed, blue-dotted, purple-dashed-dotted and orange solid lines represent the contributions of the accretion disc, companion star, discrete ejection and shock region, respectively. The green dashed line represents the compact jet, whose contribution does not vary between the two modes. The green solid line is the sum of all contributions in both SEDs. The ratio between the data points and the best-fitting model is shown at the bottom of each panel.

The switch from the high to the low mode can be connected to the occurrence of short millimetre flares. There is strong evidence that at least one such flare is directly involved in triggering the switch (see Fig. 3). We interpret these flares as the fingerprint of a discrete mass ejection that marks the removal of the inner flow (see Sect. 4.1.1), similar to what has been observed in other black-hole and NS XRBs (e.g. Cyg X-3, Circinus X-1, V404 Cygni, MAXI J1535-571, and MAXI J1820+091; Egron et al. 2021; Fender et al. 1998, 2023; Russell et al. 2019; Homan et al. 2020), albeit in different accretion regimes. As soon as the inner flow is expelled in the discrete ejection, synchrotron emission from the shock no longer takes place, resulting in a drop of the X-ray flux and the quenching of X-ray, UV and optical pulsations. The system enters the low mode.

During the low mode, the pulsar wind can still penetrate the accretion disc, collide with the in-flowing matter, and launch the compact jet, contributing to the multi-band emission in the same way as in the high mode. The X-ray emission is likely to be (almost) entirely from optically thin synchrotron radiation from this jet (Fig. 10, right panel), with a power-law index of α = −0.78 ± 0.04 (where the flux density F_ν scales with the observing frequency ν as F_ν ∝ ν^α; see Table 8). At the same time, the emission associated with the discrete ejection progressively shifts towards lower frequencies as the ejected material expands adiabatically during the low mode and quickly evolves from optically thick to optically thin (Bogdanov et al. 2018). This emission component has a steep spectrum () and provides an additional flux contribution above that of the underlying optically thick emission from the compact jet at frequencies below the millimetre band (see Figs. 1 and 3).

Eventually, the flow from the disc starts to refill the inner regions just outside the NS light cylinder on a thermal timescale (tens of seconds) due to its advection-dominated nature. The refilling of the disc re-establishes the region of shock with the pulsar wind, leading to increased emission and pulsations at X-ray, UV, and optical frequencies (Papitto et al. 2019; Miraval Zanon et al. 2022; Illiano et al. 2023). The system thus makes a switch back to the high mode.

As shown in Sect. 3.6, the low mode ingress and egress times show a statistically significant difference, with the average egress time being longer than the average ingress time. In our picture, the ingress in the low mode is caused by the ejection of the innermost regions of the disc. This will occur on a dynamic timescale. We can estimate the distance at which the innermost edge of the disc in J1023 must recede away from the NS at the high-to-low mode switch by assuming that most of the X-ray emission in the low mode is due to the jet and hence imposing that the luminosity of the synchrotron emission at the shock front is at most, say, 10% of that observed in the low mode. This gives a radial distance of r_low = [L_sd/(0.1L_X, low)]^0.5r_LC ≈ 20 r_LC ≈ 1600 km (L_sd = 4.43 × 10³⁴ erg s⁻¹ is the spin-down luminosity; Archibald et al. 2013). This value should be considered as a lower limit: at larger radii, the X-ray luminosity due to the shock would be even smaller. The dynamic timescale is s at 20 r_LC. Observations from XMM-Newton and NICER show that current X-ray instruments are capable of probing mode switching in J1023 on timescales as short as ∼10 s. In the future, proposed X-ray missions that offer instruments with much larger collecting areas, high throughput, and excellent time resolution, such as the New Advanced Telescope for High-ENergy Astrophysics (NewAthena; Nandra et al. 2013), the enhanced X-ray Timing and Polarimetry mission (eXTP; Zhang et al. 2016) and the Spectroscopic Time-Resolving Observatory for Broadband Energy X-rays (STROBE-X; Ray et al. 2019), would be able to sample the switch on a timescale of < 10 s and test our predictions.

The ejection and replenishment of the inner flow in J1023 during the switches from high to low mode and from low to high mode are reminiscent of the ejection of the X-ray corona into a jet and subsequent replenishment observed in the Galactic microquasar GRS 1915+105 (Méndez et al. 2022). A previous study has demonstrated that, in the case of GRS 1915+105, the outer disc steadily refills the central part of the system on a viscous timescale (Belloni et al. 1997). The rise time is the duration it takes for the heating front to propagate through the central disc. By assuming a subsonic velocity for the heating front, with the sound speed c_S calculated in the cool (≲6  ×  10⁶ K) disc state, we can estimate the rise time to be r_low/c_S ≈ 20 s for J1023.

In the above picture, we have assumed that the emission from the compact jet remains relatively constant over time, but it is also possible that the jet could be partially disrupted by the discrete ejection of plasma at the switch from high to low mode. However, if the ejection is not energetic enough to overcome the magnetic fields and pressures within the jet, or if it is not dense or well aligned enough to effectively interact with the jet, then the jet would quickly re-establish itself and appear as if it had never turned off. In either case, the jet would still contribute the same amount of flux in both modes, especially in the radio band where the emission is produced at greater distances from the NS (of the order of light-minutes to light-hours; Chaty et al. 2011) and therefore any changes are on much longer timescales than the mode switching.

4.1.1. Millimetre flares as the fingerprint of discrete ejecta

Starting from the timescales of the first millimetre flare that occurred on the second night (which has more accurate start and end times than the second flare) and the typical timescales of the radio flares that occurred on the first night, we attempted to reproduce the light curves of similar flares by assuming they are caused by the motion and expansion of a blob of plasma accelerated near the pulsar. Following the work of Fender & Bright (2019), we assumed that a typical discrete ejection expansion is well described by the van der Laan model (van der Laan 1966)¹³. Since no flare in this study was observed simultaneously at two different frequencies, we are unable to properly constrain the model parameters. However, we are able to check whether we could reproduce flares with peak flux densities and durations that are consistent with what we observed. We assumed that the blob is launched a few hundred kilometres away from the NS with a mildly relativistic velocity (v ≈ 0.05 c), and expands linearly over time. We also assumed a magnetic field that scales with the distance from the compact object as B ∝ R⁻¹, with a magnetic field at the launching site of the order of 10⁵ G (see Eq. (6) in Papitto et al. 2019 for an expression of the post-shock magnetic field of an isotropic pulsar wind). Additionally, we considered an electron number density of 10¹¹–10¹² cm⁻³, which is consistent with the lower limit estimated from high-resolution X-ray spectroscopy (Coti Zelati et al. 2018). Our calculations show that we can reproduce flares that last for a few seconds and reach a peak flux density of the order of 1 mJy at the ALMA Band-3 centroid frequency. Furthermore, these flares would result in flares that last for a few tens of seconds and reach a peak flux density of the order of 0.2 mJy at the VLA Band-X centroid frequency, as observed. These results are consistent with the predictions of the van der Laan model for an expanding blob of plasma.

Although we cannot firmly prove the true origin of the flares, we can ascribe the non-detection of flares in five out of the seven high-to-low-mode switches covered by ALMA to the variation in the intrinsic properties of the ejected material. These include the size, density, temperature, electron energy distribution, and opening angle of the ejected material. In general, smaller, less dense, and cooler ejecta that travel at slower bulk speeds tend to produce fainter millimetre emission (see e.g. Tetarenko et al. 2017; Wood et al. 2021). Additionally, the distribution of the ejected material along the line of sight can also affect the detectability of millimetre flares through absorption and scattering effects (this is likely the main reason for the non-detection of pulsed radio emission as well; see Sect. 4.1.2). Overall, a combination of these factors can explain the varying emission properties observed at millimetre and radio frequencies during (and following) different high-to-low-mode switches. These variations are likely to affect the time it takes for the inner flow to restore itself and produce the shock emission, resulting in the observed differences in duration among the individual low-mode episodes.

4.1.2. The non-detection of radio pulsations

The non-detection of radio pulsations could be due to ‘intrinsic’ phenomena within the pulsar magnetosphere (i.e. irrespective of the presence of external drivers such as accreting material) such as intermittency and/or nulling (see e.g. Weng et al. 2022 for recent similar considerations for the system LS I +61 303). However, these phenomena were not observed during extensive radio timing campaigns when the system was in the radio pulsar state (e.g. Archibald et al. 2013; Stappers et al. 2014), so they are unlikely to be the cause of the non-detection. Alternatively, the non-detection of pulsations may be due to the presence of matter enshrouding the system (see also Coti Zelati et al. 2014; Stappers et al. 2014). In this scenario, we can estimate a lower limit on the mass transfer rate from the companion star by requiring that the optical depth due to free-free absorption is larger than unity throughout the duration of the FAST observation. The analytical expression for the free-free optical depth is given by

(2)

where ṁ_C is the mass transfer rate in units of M_⊙ yr⁻¹, M_NS and M_C are the mass of the NS and the companion star, respectively, X and Y are the mass fraction of hydrogen and helium, respectively, T₄ is the plasma temperature in units of 10⁴ K, v₈ is the velocity of the particle wind in units of 10⁸ m s⁻¹ and P_orb is the orbital period in hours (Burgay et al. 2003; Iacolina et al. 2009). We assume M_NS = 1.65 M_⊙, M_C= 0.24 M_⊙, P_orb = 4.75 h (Archibald et al. 2009), v₈ = 1, X = 0.7 and Y = 0.3¹⁴, and v₈ = 1. In order to have τ_{ff, 1.25 GHz} ≃ 1, and therefore not detect the radio signal due to free-free scattering, we estimate ṁ_C ≳ 10⁻¹¹ M_⊙ yr⁻¹ ≈ 6 × 10¹⁴ g s⁻¹. This value is a factor of ≳10–60 larger than the mass accretion rate estimated based on a systematic analysis of the X-ray aperiodic variability of J1023 in the sub-luminous X-ray state (Linares et al. 2022). This implies that a large fraction of the accreting mass does not reach the compact object’s surface and is instead ejected from the system. It is worth mentioning that the predicted mass transfer rate for J1023, based on the transfer of mass through the loss of orbital angular momentum, is ∼ 10⁻¹¹ M_⊙ yr⁻¹ (Verbunt 1993).

4.2. Spectral energy distributions

4.2.1. Data preparation

We extracted the SEDs separately in the high and low mode and fit them using a model that reflects the physical scenario proposed in Sect. 4.1.

The NIR, optical and UV data obtained on Night 2 were de-reddened assuming A_V = 0.17 ± 0.01 mag, calculated from the XMM-Newton estimate of N_H using the N_H/A_V relation reported by Foight et al. (2016). The absorption coefficients in the other bands were estimated using the equations reported by Cardelli et al. (1989), and are listed in the last column of Table 7.

For the low-mode SED, we considered optical data acquired between MJD (Modified Julian Day) 59392.01085 and 59392.01434, totaling three consecutive 60 s integration images in each band; the J-band flux is the average of the fluxes measured with NTT/SOFI while the source was in the low mode, according to our XMM-Newton monitoring; the UV (uvm2) flux was measured on MJD 59392.01302; radio flux densities are taken from the VLA observation from 2021 June 3–4; the millimetre flux density was derived from stacking all the ALMA images acquired during the low modes. No K-band observations have been acquired during the low mode. The unabsorbed X-ray flux was calculated as the average flux during all periods of low mode.

For the high-mode SED, we considered K-band data acquired between MJD 59392.00237 and 59392.00289, and then from MJD 59392.00593 to 59392.00832 (i.e. all the images in which J1023 is detected, 15 × 15 s); the J-band flux is the average of all the fluxes measured with NTT/SOFI while the source was in the high mode, according to our XMM-Newton monitoring; for the optical griz bands, we considered one 60 s image for each band, at MJD 59392.00528; the UV (uvm2) flux was acquired on MJD 59392.00413; radio flux densities are taken from the VLA observation on 2021 June 3–4; the millimetre flux density was derived after stacking all the ALMA images acquired during the high mode episodes. The unabsorbed X-ray flux was calculated as the average flux during all periods of high mode.

4.2.2. Model setup

Low mode. To perform the fit to the SED, we used a model comprising a blackbody of an irradiated star and a multi-colour accretion disc, together with a broken power law to describe the optically thick and thin components of the synchrotron emission from the compact jet. In addition, we considered a power law to model the optically thin emission from the discrete ejecta at all frequencies. For the stellar component, we considered an irradiated star model (see Eqs. (8) and (9) of Chakrabarty 1998) where the surface temperature of the star (T_*) is allowed to vary and the distance to J1023 (D), the radius of the companion star (R_c), the binary separation (a), the irradiation luminosity (L_irr) and the albedo of the star (η_*) are fixed to known (or reasonable) values: D = 1368 pc (Deller et al. 2012), R_c = 0.43 R_⊙ (Archibald et al. 2009), η_* = 0.1, and , where M_NS = 1.65M_⊙ is the NS mass, M_C = 0.24 M_⊙ is the mass of the companion star, P_orb = 4.75 h is the binary orbital period and G is the gravitational constant. L_irr was fixed to 6.5 × 10³³ erg s⁻¹, the value estimated by Shahbaz et al. (2019) from the difference of the temperatures of the illuminated and dark sides of the companion star of J1023. This value corresponds to 14% of the pulsar spin-down luminosity (Archibald et al. 2013), and is of the same order of magnitude as the X-ray and gamma-ray luminosities in the sub-luminous X-ray state. For the multi-colour accretion disc model (Eqs. (10)–(15) by Chakrabarty 1998), we allowed the inner radius of the accretion disc where the optical/UV emission becomes relevant (r_in, opt/UV) to vary, while we fixed the X-ray albedo of the disc to 0.95 (Chakrabarty 1998) and the mass transfer rate to 10⁻¹¹ M_⊙ yr⁻¹. This is the value predicted for an XRB system with a short orbital period hosting a main-sequence star and transferring mass through loss of angular momentum (Verbunt 1993), and is consistent with the lower limit estimated from the non-detection of radio pulsations (Sect. 2.6). The broken power-law component has the normalisation F₀, the two slopes α₁ and α₂ and the break frequency ν_break as free parameters; the power law associated with the emission from the discrete ejecta has the normalisation F_{thin , 10} (calculated at a frequency of 10 GHz) and the slope α_thin as free parameters.

High mode. We used a more complicated model following the physical picture explained above. In particular, we considered that the X-rays are produced by the sum of the optically thin component of the synchrotron emission from the base of the compact jet and the synchrotron emission from the region of shock between the pulsar wind and the inner flow. For the latter, we used the analytical expression for the synchrotron spectrum by Chevalier (1998; see his Eq. (1)). For the emission at lower frequencies, we included the blackbody of the irradiated star (fixing L_irr to 6.5 × 10³³ erg s⁻¹), the multi-colour blackbody of the accretion disc, and the optically thick component of the synchrotron emission from the compact jet. The parameters of the fit are T_*, ν_break, α₁, α₂, F₀, r_in, opt/UV, ν_sync, α_sync and F_sync, where ν_sync, α_sync, F_sync are the peak frequency of the synchrotron emission, the slope of the optically thin synchrotron emission and the normalisation of the synchrotron emission, respectively.

4.2.3. Model fit

Employing both datasets (high and low mode), we performed an MCMC sampling of the posterior probability density function for the parameter space of the two models simultaneously (since some of the free parameters are shared by both models). We used a Gaussian prior for the temperature of the star centred on the value measured by Shahbaz et al. (2022) during the disc state, T_* = (6128 ± 33) K. Gaussian priors have also been used for the parameters α₁ (centred on 0.2 ± 0.2) and α₂ (centred on −0.7 ± 0.2), following the results of Bogdanov et al. (2018). We note that we also tried wider Gaussian priors for these parameters, resulting in wider posterior distributions with medians compatible with the results below within 1–2σ. For all other parameters, we used uniform priors in very wide intervals to sample a large region of the parameter space (see the last column of Table 8). This is of great importance to minimise the probability that the actual source parameters fall outside the region sampled by the chosen priors.

4.2.4. Results

The results of the fit are reported in Table 8 and are summarised in the following paragraphs. All the parameters have been estimated as the median of the marginal posterior distributions, with 1σ credible intervals coming from the 16th–84th percentiles of the posterior distributions (see Fig. 11 for the histograms of the parameter samples). Most of the marginal posterior distributions of the parameters shown in Fig. 11 are nearly Gaussian, meaning that they are almost symmetrical. Figure 11 also shows a strong correlation between ν_sync and F_sync, and between F₀ and α₁. This is however not surprising considering the scarce number of points at the lowest frequencies and the impossibility to constrain quantities such as ν_break and ν_sync, due to numerous overlapping components at infrared and optical frequencies. In particular, we note that F_sync, ν_sync and α_thin have a non-symmetrical posterior distributions, with long tails towards the upper or the lower limit of the chosen interval for the prior, depending on the parameter. However, the peak of their distributions is well determined (see Fig. 11). To estimate the goodness of the fit, we used the Bayesian p-value (reported in the last line of Table 8; Lucy 2016). The p-value is 0.78 for the high mode and 0.07 for the low mode, meaning that the results of the fit are indeed acceptable (with a slight indication of overestimated and underestimated uncertainties for the high and low mode fit, respectively).

Fig. 11. Corner plot showing the posterior probability distributions of the parameters obtained from the MCMC sampling algorithm. The solid lines in the panels along the diagonal represent the probability density of each parameter. The vertical dashed lines indicate the 16th, 50th, and 84th percentiles of the distributions. The median of the parameter posterior and the 1σ error bars are displayed above each panel. The off-diagonal panels show the 2D posterior distributions, with contour lines indicating the 1σ, 2σ, and 3σ equivalent bounds.

The synchrotron emission caused by the shock has a higher flux density in the optical band than the average pulsed flux density fraction of ∼1% estimated by Papitto et al. (2019), which makes our scenario plausible. In particular, we find that the shock contributes ∼3% of the emission in the optical band in the high mode, which is similar to the contribution of the optically thin component of the compact jet emission. In the XMM-Newton/EPIC 0.3–10 keV range, the shock accounts for ∼83% of the emission, while the jet only contributes ∼17% of the flux in the high mode.

4.3. Optical polarisation measurements and predictions of X-ray polarisation

The average optical polarisation reported in this work is similar to levels previously measured for this source (Baglio et al. 2016; Hakala & Kajava 2018). Since the polarisation level is < 1%, it could have multiple origins. Given the lack of modulation of the Stokes parameters with the orbital period, it is unlikely that the source of the polarisation is scattering of radiation in any region of the system (such as the surface of the accretion disc), in contrast to previous suggestions (Baglio et al. 2016). We identify two possible main emission mechanisms for the polarised emission, both of which are non-thermal: optically thin synchrotron emission from the compact jet, or synchrotron emission at the shock between the pulsar wind and the accretion inflow, which also gives rise to the observed optical pulsations (as also suggested by Hakala & Kajava 2018). The jet flux density in the r band is ∼0.02 mJy, corresponding to a fractional contribution to the total flux in the same band of ∼2.4% and ∼2% in the high and low modes, respectively. In the high mode, where the polarisation level measured in the R band is 0.65 ± 0.43%, this translates into a polarisation level of the jet emission of ∼26% (or ∼37% if the average polarisation level of 0.92% is considered instead). In the low mode, the jet emission is instead intrinsically polarised at a level of ∼17%. Following Shahbaz (2019) and references therein, we estimate a level of ordering of the magnetic field lines of ∼30% (considering the case of an optically thin synchrotron emission with a spectral slope of ∼–0.78, as obtained from our fit and averaging across modes; see Table 8). This indeed indicates a highly ordered magnetic field and suggests that tMSPs might have a more ordered magnetic field at the base of the jet than other NS and black hole XRBs (Russell & Fender 2008; Baglio et al. 2014, 2020; Russell et al. 2016). Based on our picture, this could be ascribed to the presence of a continuous, steady flow of plasma near the pulsar, which may help maintain a more ordered magnetic field structure at the base of the jet. Assuming that the observed optical polarisation is due to this emission component and that the same population of accelerated electrons in the jet produces optically thin synchrotron radiation at all frequencies above the break frequency, we can predict a polarisation level in the X-ray band – in the 2–8 keV energy range of the Imaging X-ray Polarimetry Explorer (IXPE; Weisskopf et al. 2016; Soffitta et al. 2021) due to the jet of ≈17% in the low mode, where the jet contributes to ∼100% of the flux according to our SED modelling.

Alternatively, the observed optical polarisation might originate from synchrotron radiation at the region of shock between the pulsar wind and the accretion flow. In this scenario, polarisation should only be detected during the high mode. Unfortunately, the significance of our polarisation results obtained in the high and low modes separately is not high enough to verify this, as the results are compatible with each other within 1σ ( in the low mode and 0.65 ± 0.43% in the high mode). The fractional flux contribution of the shock in the r-band according to the SED modelling in the high mode is ∼4.8%, which translates into a LP of the emission of the shock of ∼15–20%, with peaks of up to 65–90% at the pulsations (in fact, the fraction of pulsed radiation is typically of ∼1%; Papitto et al. 2019). This prediction could be tested using next-generation optical high-time resolution polarimetric instruments (Ghedina et al. 2022). In the X-rays instead, we predict a fractional flux contribution of the shock emission of ∼90% in the IXPE energy range, which converts into an intrinsic LP level of ≈12–17%. This could be detected in an IXPE observation of ∼460 ks (assuming J1023 spends 70% of the time in the high mode).

In both scenarios, the validity of these predictions depends on the assumption that the optical polarisation is solely caused by synchrotron radiation from the jet or shock. However, if the observed polarisation is due to a combination of these factors or if other processes are at play that are not being considered (such as scattering on the disc surface contributing to the modulation with a low S/N in our measurements), these predictions may need to be revised.

4.4. The millimetre continuum emission of J1023 and other X-ray binary systems

We have reported the detection of millimetre continuum emission for the first time from J1023 (and from a tMSP in general) using ALMA observations at ≈100 GHz. The millimetre emission is unresolved, implying a size for the emission region smaller than ≈6 × 10¹⁵ cm (this corresponds to ≈7 × 10⁴ times the orbital separation of the system; Archibald et al. 2013). Millimetre/sub-millimetre continuum emission has been firmly detected from a number of transiently accreting stellar-mass black holes (Paredes et al. 2000; Ogley et al. 2000; Fender et al. 2000, 2001; Ueda et al. 2002; Russell et al. 2013a, 2020; van der Horst et al. 2013; Tetarenko et al. 2015, 2017, 2019, 2021; de Haas et al. 2021; Koljonen & Hovatta 2021; Díaz Trigo et al. 2021) and in three accreting NSs so far: the persistent systems 4U 1728−34 and 4U 1820−30 (Díaz Trigo et al. 2017) and the transient accreting millisecond pulsar Aql X-1 (Díaz Trigo et al. 2018). In these cases, however, such emission has been detected during accretion states that are seemingly different from the X-ray sub-luminous disc state attained by J1023 at the epoch of our ALMA observations, as detailed below.

In black hole systems, millimetre continuum emission is typically detected during the rising phase of an X-ray outburst or along the outburst decay, during the ‘hard’ X-ray spectral state. In this state, the X-ray emission is typically dominated by a power-law component that is believed to arise from Compton up-scattering of soft disc photons by a population of hot thermal electrons located either in the innermost region of the accretion flow or in an extended cloud above the disc (for a recent review, see Motta et al. 2021). At low frequencies, a flat to slightly inverted optically thick spectrum (α > 0) has been observed from radio through sub-millimetre frequencies and above and has been associated with partially self-absorbed synchrotron radiation from electron populations in different regions along a compact jet. This spectrum becomes optically thin (α < 0) typically around infrared frequencies, resulting in a spectral break (Russell et al. 2013a,b). On the other hand, at very high mass accretion rates (typically close to the outburst peak), these systems can be in a ‘soft’ X-ray spectral state where the X-ray emission is dominated by a thermal component produced in the hot inner regions of the accretion disc (see Motta et al. 2021 and references therein). In the state that immediately precedes or follows the soft state during the outburst (the so-called soft-intermediate state), the optically thick jet emission is observed to be quenched altogether, while discrete jet ejecta characterised by an optically thin steep spectrum (α < 0) in the radio band can be launched (Fender et al. 2004; Russell et al. 2020) and be accompanied by strong radio variability and bright multi-frequency flaring activity (even in the sub-millimetre band; Tetarenko et al. 2019).

The picture for NS systems is more diverse (see van den Eijnden et al. 2021 for a recent overview). All three of the above-mentioned NS systems were observed in the millimetre band with ALMA, along with observations at other frequencies, in regimes of high mass accretion rates (L_X > 5 × 10³⁶ erg s⁻¹). 4U 1728−34 was observed during a transition from a hard to a soft state. The low-frequency spectrum was compatible with synchrotron emission from a compact jet and also revealed the presence of a break from optically thick to optically thin synchrotron radiation at frequencies in the range 10¹³–10¹⁴ Hz. In the case of 4U 1820−30, enhanced millimetre emission was detected during a soft X-ray state (the first time millimetre emission has been observed for an accreting NS in the soft state). The spectrum from the radio to the sub-millimetre bands was slightly inverted and showed no evidence for the presence of a spectral break (Díaz Trigo et al. 2017). This indicates emission from a compact jet and supports earlier findings that jets in NS systems do not all appear to be completely suppressed in soft accretion states (Migliari et al. 2004). Aql X-1 was observed five times using ALMA over a time span of ∼1.5 months during the decay phase of an X-ray outburst in 2016, covering multiple transitions between distinct accretion states. The spectrum from the radio to the millimetre bands was again compatible with emission from a compact jet and, for the first time in an accreting NS in a low-mass XRB, a break in the SED was detected at a frequency that varies dramatically depending on the accretion state. Specifically, the break frequency initially decreased from ∼100 GHz to less than ∼5 GHz during a transition from a hard to a soft state, and then increased again up to a frequency within the range ∼30–100 GHz during a transition to a hard state at later stages of the outburst decay (Díaz Trigo et al. 2018).

Unlike all the cases mentioned above, the ALMA observations of J1023 were performed at an epoch when the system did not display any changes in its spectral state and was lingering in a sub-luminous X-ray state (≈10³³ erg s⁻¹). It is worth noting that both the average and the mode-resolved X-ray spectra of J1023 have always been hard and have not shown any substantial changes over the past ten years or so during the sub-luminous X-ray state.

We can compare the millimetre and X-ray luminosities of J1023 with those measured for black-hole and NS systems using quasi-simultaneous millimetre and X-ray observations. These luminosities are reported in Table 9 and plotted in Fig. 12. Despite the relatively low amount of existing millimetre/sub-millimetre data for XRBs, Fig. 12 shows that millimetre continuum emission has been detected in systems that attain a wide range of X-ray luminosities (∼10³³–10³⁹ erg s⁻¹) and that J1023 currently holds the record as the faintest XRB firmly detected in the millimetre band.

Fig. 12. Millimetre and X-ray luminosities for accreting compact objects in XRBs. Black hole systems are marked in black. The adopted values for the luminosities are listed in Table 9. The dashed black line and the grey shaded area show the best-fitting slope for the black hole correlation and the region on the plane at less than 3σ from the best-fitting slope, respectively (see the main text for details).

The mechanism responsible for the X-ray emission in the two modes in J1023 is clearly different from that at work in black hole systems in the hard state, so a proper comparison cannot be made. With this in mind, we assessed how the position of J1023 on the L_mm–L_x plane compares with those of black hole systems by parameterised the dependence of L_mm upon L_x as a relation of the form Log L_mm − Log L_mm, ave = α + β(Log L_X − Log L_X, ave)¹⁵. We then investigated the shape of the L_mm – L_x correlation for black hole systems by performing a Bayesian-based MCMC sampling of the black hole data (including upper limits) using the linear regression algorithm LINMIX₋ERR (Kelly 2007). There are three free parameters: the intercept α, the slope β and σ₀, a parameter that accounts for an intrinsic random (Gaussian) scatter of the values around the best-fit model. We calculated the median values of the fit parameters from 10⁴ draws from the posterior distribution, and determined the 1σ credible intervals of these parameters using the 16th–84th percentiles of the posterior distributions. We derived α = −0.04  ±  0.15, β = 0.46  ±  0.09 and σ₀ = 0.18  ±  0.13. The best-fitting slope is consistent with that derived for the L_radio–L_x correlation for black hole systems (Gallo et al. 2018) within the uncertainties. Figure 12 shows the best-fitting slope (black dashed line) and the region on the plane where the values of L_mm and L_X are located at a distance of less than 3σ from the best-fitting slope (grey shaded area). This figure shows that, during the low mode, J1023 follows the same correlation as black hole systems in the hard state (Bogdanov et al. 2018).

5. Summary and conclusions

We have presented the results of the largest multi-wavelength campaign ever conducted on the tMSP PSR J1023+0038. This campaign took place in June 2021 while the system was in an active low-luminosity X-ray disc state. Our key findings can be summarised as follows:

− On the first night of observations, flaring activity was detected at X-ray and optical wavelengths, followed by mode-switching behaviour in X-rays. The UV emission was significantly fainter in low mode, while the radio continuum emission was brighter. No evidence of pulsed radio emission was detected in either mode. Our observations provide the most stringent limits on the pulsed radio flux density to date.

− On the second night of observations, a low-significance variation in the optical polarised signal was detected during mode changes, with a slightly higher LP level in the high mode compared to the low mode. The millimetre emission did not strictly follow the X-ray low modes, except for two brief flares. One of these flares occurred at the high-to-low mode switch, while the exact onset time of the other flare could not be determined due to a lack of observational coverage.

− These results allow us to draw a physical picture of high-low mode switches in J1023, which involves a rotation-powered pulsar, an accretion disc, and discrete mass ejections on top of a relativistic compact jet. During the high mode, synchrotron radiation at the shock front between the pulsar wind and the inner accretion flow produces most of the X-ray emission, as well as X-ray, UV, and optical pulsations. The switch to low mode is triggered by discrete mass ejections, resulting in a drop in the X-ray flux and a quenching of pulsations. During the low mode, the pulsar wind still penetrates the accretion disc and launches a compact jet, contributing to the multi-band emission in the same way as in the high mode. Eventually, the flow from the disc refills the inner regions outside the NS light cylinder, leading to increased emission and pulsations at X-ray, UV, and optical frequencies. The system then switches back to the high mode. This scenario is supported by the modelling of the SEDs extracted separately during the high and low modes.

− Relativistic pulsar winds can shape the structure and dynamics of accretion flows in XRB systems that contain quickly spinning NSs and can be used to probe emission mechanisms in tMSPs and other XRB systems.

− Black hole systems and tMSPs share intriguing similarities in their phenomenological properties, which emphasises the need for further research to deepen our understanding of accretion physics in compact objects.

These results show how multi-wavelength campaigns that combine diverse observational techniques provide powerful tools for unravelling the nature of elusive objects such as tMSPs.

Acknowledgments

We thank the referee for helpful comments. We thank the ALMA staff, particularly Edwige Chapillon, for the help with scheduling the observations; the NICER PI, Keith Gendreau, for approving our ToO request and the operation team for executing the observations; Brad Cenko and the Swift duty scientists and science planners, for making the Swift Target of Opportunity observation possible; Nando Patat, for the help with scheduling the NTT/SOFI observations; Giulia Illiano, for providing checks on timing analysis of X-ray data; Rob Fender, for providing useful notebooks on synchrotron radiation from variable radio sources; Nanda Rea, for helpful comments on the manuscript. The results reported in this study are based on observations obtained with XMM-Newton (PI: Campana); NuSTAR (PI: Campana); the NASA/ESA Hubble Space Telescope (obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are associated with program 16061; PI: Campana); FAST (PI: Hou); the REM Telescope, INAF Chile (PI: Baglio). Part of the observations presented in this study were obtained at the European Southern Observatory under ESO programme 107.22RK.001 (PI: Baglio). This paper makes also use of the following ALMA data: ADS/JAO.ALMA#2019.A.00036.S (PI: Coti Zelati). XMM-Newton is an European Space Agency (ESA) science mission with instruments and contributions directly funded by ESA member states and NASA. NuSTAR is a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory and funded by NASA. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. FAST is a Chinese national mega-science facility, operated by National Astronomical Observatories of the Chinese Academy of Sciences (NAOC). NICER is a 0.2–12 keV X-ray telescope operating on the International Space Station, funded by NASA. The Neil Gehrels Swift Observatory is a NASA/UK/ASI mission. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The XMM-Newton SAS is developed and maintained by the Science Operations Centre at the European Space Astronomy Centre. The NuSTAR Data Analysis Software (NuSTARDAS) is jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (Caltech, USA). The data that support the findings of this study are publicly available at their respective online archive repositories (XMM-Newton: http://nxsa.esac.esa.int/nxsa-web/; NICER: https://heasarc.gsfc.nasa.gov/FTP/nicer/data/obs; NuSTAR: https://heasarc.gsfc.nasa.gov/W3Browse/all/numaster.html; Swift: https://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/w3browse.pl; HST: http://hst.esac.esa.int/ehst/; VLA: https://data.nrao.edu/portal/; VLT/FORS2: http://archive.eso.org/scienceportal/; ALMA: http://almascience.nrao.edu/aq/; REM: http://ross.oas.inaf.it/REMDB/public.php; NTT/SOFI: http://archive.eso.org/scienceportal/). FAST data can be accessed through the FAST Data Center: https://fast.bao.ac.cn/cms/article/11/. F.C.Z. is supported by a Ramón y Cajal fellowship (grant agreement RYC2021-030888-I) and Catalan grant SGR-Cat 2021 (PI: Graber). F.C.Z. and D.F.T. are supported by the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. F.C.Z., S.C.a., P.D.’A., F.A., P.C., D.d.M. and A.P. acknowledge financial support from INAF-Fundamental research astrophysics project “Uncovering the optical beat of the fastest magnetised neutron stars” (FANS; PI: AP). S.C.a., S.C.o. and P.D.’A. acknowledge support from ASI grant I/004/11/5. G.B. acknowledges support from the PID2020-117710GB-I00 grant funded by MCIN/ AEI /10.13039/501100011033. A.M.Z. is supported by PRIN-MIUR 2017 UnIAM (Unifying Isolated and Accreting Magnetars, PI S. Mereghetti). X.H. is supported by the National Natural Science Foundation of China through grant No. 12041303. J.L. is supported by the National Natural Science Foundation of China through grant No.12273038. D.d.M. and A.P. acknowledge financial support from the Italian Space Agency (ASI) and National Institute for Astrophysics (INAF) under agreements ASI-INAF I/037/12/0 and ASI-INAF n.2017- 14-H.0, from INAF ‘Sostegno alla ricerca scientifica main streams dell‘INAF’, Presidential Decree 43/2018 and from INAF ‘SKA/CTA projects’, Presidential Decree 70/2016. D.M.R. and D.M.B. acknowledge the support of the NYU Abu Dhabi Research Enhancement Fund under grant RE124. D.F.T. is supported by the grants PID2021-124581OB-I00 funded by MCIN/AEI/10.13039/501100011033, 2021SGR00426 of the Generalitat de Catalunya and by MCIN with funding from European Union NextGeneration EU (PRTR-C17.I1). M.M. acknowledges the research programme Athena with project number 184.034.002, which is (partly) financed by the Dutch Research Council (NWO) and thanks the Team Meeting at the International Space Science Institute (Bern) for fruitful discussions. F.M.V. acknowledges support from the Spanish Ministry of Science for the grant PID2020-120323GB-I00 and for the FJC2020-043334-I grant financed by MCIN/AEI/10.13039/501100011033 and Next Generation EU/PRTR. This work was also partially supported by the COST Action ‘PHAROS’ (CA 16124). Software: APLPY v.2.1.0 (Robitaille & Bressert 2012); Astropy v.5.3 (Astropy Collaboration 2013, 2018); astropy:2022; CASA v.6.5.2 (CASA Team et al. 2022); CIAO v.4.15 and CALDB v.4.10.4 (Fruscione et al. 2006); CORNER PLOT v.2.2.1 (Foreman-Mackey 2016); DAOPHOT (Stetson 1987), part of the Starlink software (Currie et al. 2014; Berry et al. 2022); HEASOFT v.6.31.1 (NASA High Energy Astrophysics Science Archive Research Center (HEASARC) 2014); IRAF v.2.17 (https://github.com/iraf-community/iraf); LINMIX₋ERR (https://github.com/jmeyers314/linmix); MATPLOTLIB v.3.7.1 (Hunter 2007); NICERDAS v.10a; NUMPY v.1.24.0 (Harris et al. 2020); NUSTARDAS, v.2.1.2; PRESTO v.4.0 (https://github.com/scottransom/presto); SAOImageDS9 v.8.4 (Joye & Mandel 2003); SAS v.20.0 (Gabriel et al. 2004); SCIPY v.1.10.1 (Virtanen et al. 2020); SRPAstro v.4.8.0 (https://pypi.org/project/SRPAstro/); XRONOS v.5.22 (Stella & Angelini 1992); XSPEC v.12.13.0g (Arnaud 1996). Upon request, the corresponding authors (M.C.B. and F.C.Z.) will provide the code used to produce the figures.

References

Appendix A: Searches for counterparts to ALMA J102349.1+003824

We searched for a counterpart to ALMA J102349.1+003824 at other wavelengths using both simultaneous and archival observations of the field. No UV emission was detected at the source position in the Swift/UVOT observations simultaneous with ALMA on the night of 2021 June 26–27, either in the images collected separately during the three snapshots or in the stack of these images. We set a 3σ upper limit on the magnitude of a putative counterpart of UVM2> 21 mag from the stacked image (total exposure of ∼2.4 ks). A few weeks before the ALMA observations were acquired, the FOV of ALMA J102349.1+003824 was also observed in the optical (g′, r′ and i′ SDSS filters) with the Las Cumbres Observatory (LCO) network of telescopes. No counterpart was detected in these data. We set the following 3σ upper limits on the optical magnitude using 60 s images acquired with the 1m telescope of the LCO network located in South Africa on 2021 May 8: i′> 17.8 mag, r′> 17.9 mag, g′> 18.6 mag. There is also no optical/NIR counterpart in either the third data release (DR3) from the Gaia space telescope or in the 2MASS point source catalogue. ALMA J102349.1+003824 was also not detected in the VLA observation carried out on 2021 June 3, down to a 3σ flux density upper limit of 12 μJy at 10 GHz.

Several observations of the field have been performed in the soft X-ray band over the past years, mainly using Swift and XMM-Newton. However, given the X-ray luminosity of J1023 and its small angular distance from ALMA J102349.1+003824 (∼22″ ), the point-spread function of the X-ray instruments on board these observatories encompasses the position of ALMA J102349.1+003824, preventing a conclusive assessment of the presence or absence of an X-ray counterpart to this source. We then analysed a 83-ks long Chandra/ACIS-S (Advanced CCD Imaging Spectrometer) observation performed in March 2011 (Obs. ID 11075) to search for X-ray emission at the position of ALMA J102349.1+003824. We reprocessed the data using the task CHANDRA₋ REPRO within the CIAO software and ran the WAVDETECT source detection tool. No evidence for an X-ray counterpart to ALMA J102349.1+003824 was found. Using SRCFLUX, we set a 3σ upper limit on the net count rate of 1.3 × 10⁻⁴ counts s⁻¹ over the 0.5–7.0 keV energy range. Assuming an absorbed power-law spectrum with an absorption column density of N_H = 5 × 10²⁰ cm⁻² (i.e. the one expected within the Galaxy in the direction of the source: Willingale et al. 2013) and a photon index in the range Γ = 1 – 4, the limit above translates into an unabsorbed flux of < (2 − 4)×10⁻¹⁵ erg cm⁻² s⁻¹ (depending on the assumed power-law slope) in the 0.5–7 keV range.

All Tables

All Figures

Fig. 1. Temporal evolution of the X-ray, UV, optical, NIR, and radio emissions of J1023 on the night of 2021 June 3–4. The light curves are shown in decreasing order of energy band from top to bottom, with the XMM-Newton/EPIC one as the reference for all others (NuSTAR, HST/STIS, XMM-Newton/OM, NTT/SOFI, and VLA). The XMM-Newton X-ray light curve is shown over the time interval in which the three EPIC instruments collected data simultaneously. The OM residual count rates were obtained after removing the orbital modulation from the time series (see the main text for details). For displaying purposes, only VLA detections are shown. The dark yellow segment in the bottom panel indicates the time range covered by the FAST observation. The orange shaded areas mark the time intervals of the low X-ray mode detected by XMM-Newton/EPIC covered by observations with other instruments.

In the text

	Fig. 2. Zoom-in of the time series collected on Night 1 over the time span of maximum overlap among the different telescopes. Arrows indicate 3σ upper limits.
In the text

Fig. 3. Temporal evolution of the X-ray, UV, optical, and millimetre emissions and of the optical LP of J1023 on the night of 2021 June 26–27. The light curves are shown in decreasing order of energy band from top to bottom (NICER, Swift/XRT, Swift/UVOT, REM, VLT/FORS2 polarisation, and ALMA). In the bottom panel, the ALMA time series are plotted using either bins of length 4–5 min (cyan) or bins of variable length corresponding to the duration of the X-ray mode episodes identified by the strictly simultaneous NICER and Swift data (blue). In some panels, the sizes of the markers are bigger than the uncertainties and/or the time bins. Arrows indicate 3σ upper limits. The orange shaded areas mark the time intervals of low X-ray mode detected by NICER and Swift/XRT covered by observations with other instruments. A plot of the X-ray and millimetre light curves around the epoch of the millimetre flares is shown in Fig. 6.

In the text

	Fig. 4. Stokes Q and U measurements for the whole VLT/FORS2 dataset. Black dots, red squares, and the green star indicate the Stokes parameters of J1023, a comparison star in the field, and the non-polarised standard star, respectively.
In the text

Fig. 5. ALMA images of the field around J1023. Left: ALMA 3.1 mm continuum image showing the two detected sources, J1023 and ALMA J102349.1+003824. The white dotted boxes display the FOV shown in the right panels. The ALMA synthesised beam (0.87″ × 0.53″, PA = −74.5°) is shown at the bottom-right corner. Right: Close-up images of J1023 (top) and ALMA J102349.1+003824 (bottom). In both panels, the colour scale corresponds to the image with natural weighting, with the 3σ contour indicated in white (σ values are listed in Table 2). The black contours are derived from the ALMA image obtained with robust=0. Contours range from 3σ to 9σ in steps of 1σ, where σ = 9.1 μJy beam−1 for J1023 and σ = 15 μJy beam−1 for ALMA J102349.1+003824. The white cross depicts the radio position of J1023 (Deller et al. 2012). A faint emission tail is detected at the 3σ level extending south-west of J1023. The synthesised beams for the robust=0 and natural images are displayed in the bottom-left and bottom-right corners, respectively (see Table 2). In all panels, the colour scale is in units of Jy beam−1.

In the text

Fig. 6. X-ray (top) and millimetre (bottom) time series collected around the epochs when the two millimetre flares were detected. The X-ray light curve is binned at a time resolution of 10 s, and the millimetre light curve is binned at 15 s. Arrows indicate 3σ upper limits. Note the lack of coverage in the millimetre band for Flare 2 before and during the switch from the high to the low mode in the time interval between two consecutive ALMA scans. Note also the different scales on the horizontal axis in the panels on the left and on the right.

In the text

Fig. 7. Contour levels for the millimetre emission detected during 15 s time intervals at the peak of Flare 1 (left panel) and Flare 2 (right panel) shown with solid white lines overlaid on the naturally weighted ALMA 3.1 mm image extracted using the whole integration time (colour scale). These contour levels start at 3σ and increase in steps of 1σ, where σ is the rms of the image, ∼90 μJy beam−1. Both panels display the same colour scale range in units of Jy beam−1. Time intervals (in UTC) are displayed at the top of each panel. The red cross depicts the radio position of J1023 (Deller et al. 2012). The synthesised beams, with extension 0.696″ × 0.540″ and position angle PA = −79.92° for Flare 1 and 0.931″ × 0.525″ and PA = −74.21° for Flare 2, are shown at the bottom-right corner of each panel.

In the text

Fig. 8. Results of the statistical analysis of the X-ray mode switching timescales. Left: Duration of the low-to-high mode switch (egress) as a function of the duration of the high-to-low mode switch (ingress) computed using data from XMM-Newton/EPIC (in red) and NICER (green). The circled green points represent the values for the ingress associated with the two millimetre flares detected with ALMA and the subsequent egress. The solid black line indicates where the two quantities are equal. Right: Distribution of the low-mode ingress (in green) and egress (in red) times in logarithmic scale from XMM-Newton and NICER data.

In the text

Fig. 9. Schematic visual representation of the evolution of the inner flow into an outflow at the switch from the high to the low X-ray emission mode. When the system is in the high mode (top panel), a small-size inner flow is present, together with a faint steady jet that is launched along the pulsar rotational axis and gives rise to the observed low-level radio and millimetre emission. As the pulsar rotates, the pulsar wind (marked with solid green lines) wobbles around the equatorial plane (see e.g. Bogovalov 1999) and shocks off the electrons in the inner flow at two opposite sides (red spots) at a distance that is slightly larger than the light cylinder radius (≃80 km). At each pulsar rotation, synchrotron emission at the shock at X-ray, UV, and optical frequencies is modulated at the spin period at one side (bright red spot), while it is absorbed by material in the inner flow at the other side (light red spot; Papitto et al. 2019). When the system makes a switch to the low mode (bottom panel), a bright discrete ejection is launched along the pulsar rotational axis on top of the compact jet, the inner flow disappears, and the shock emission is quenched. For displaying purposes, only the inner regions of the compact jet are represented in both panels.

In the text

Fig. 10. Broadband (radio to X-rays) SEDs of J1023 during the high and low modes, along with the best-fitting models. Both SEDs were corrected for extinction effects. The red-dashed, blue-dotted, purple-dashed-dotted and orange solid lines represent the contributions of the accretion disc, companion star, discrete ejection and shock region, respectively. The green dashed line represents the compact jet, whose contribution does not vary between the two modes. The green solid line is the sum of all contributions in both SEDs. The ratio between the data points and the best-fitting model is shown at the bottom of each panel.

In the text

Fig. 11. Corner plot showing the posterior probability distributions of the parameters obtained from the MCMC sampling algorithm. The solid lines in the panels along the diagonal represent the probability density of each parameter. The vertical dashed lines indicate the 16th, 50th, and 84th percentiles of the distributions. The median of the parameter posterior and the 1σ error bars are displayed above each panel. The off-diagonal panels show the 2D posterior distributions, with contour lines indicating the 1σ, 2σ, and 3σ equivalent bounds.

In the text

Fig. 12. Millimetre and X-ray luminosities for accreting compact objects in XRBs. Black hole systems are marked in black. The adopted values for the luminosities are listed in Table 9. The dashed black line and the grey shaded area show the best-fitting slope for the black hole correlation and the region on the plane at less than 3σ from the best-fitting slope, respectively (see the main text for details).

In the text