© The Authors 2025

1. Introduction

Millisecond pulsars (MSPs) are the fastest spinning neutron stars (NSs) in the Universe. They achieve their millisecond spin period during an accretion phase lasting billions of years. Throughout this time, a low-mass (≲1 M_⊙) donor star transfers matter and angular momentum to the NS surface, causing the system to shine as a bright low-mass X-ray binary (LMXB). Once the mass accretion rate ceases, the NS transitions to a radio pulsar state driven by the gradual loss of its rotational energy (Alpar et al. 1982; Radhakrishnan & Srinivasan 1982). The MSPs found in tight (i.e., with an orbital period P_orb < 1 d) binary systems provide crucial clues for investigating the complex mechanisms behind mass accretion and ejection as well as the processes of magnetospheric particle acceleration.

Among MSPs, an intriguing sub-class is represented by transitional MSPs (tMSPs; Papitto & de Martino 2022, and references therein). These unique objects have conclusively bridged the long-sought evolutionary gap between accreting NSs and radio MSPs in binary systems. Remarkably, tMSPs can swing back and forth between accretion-powered and rotation-powered states on timescales accessible within a human lifetime, from years to decades. All confirmed tMSPs and a few candidates have also been caught in a third state, which based on its characteristics appears to be an intermediate between the two standard states discussed above. This intermediary is referred to as the sub-luminous disk state, and it is characterized by the presence of an accretion disk, an X-ray luminosity that is fainter than in the fully accreting state but brighter than in the rotation-powered phase, and a γ-ray luminosity up to ten times brighter than that of MSPs in their state fed by rotation (see, e.g., Torres & Li 2022, and references therein). The most defining feature is the variable X-ray emission that unpredictably switches between two distinct intensity levels, dubbed ‘high’ and ‘low’ modes (Papitto et al. 2013; Linares 2014), along with occasional flares (e.g., Tendulkar et al. 2014; Bogdanov et al. 2015). For confirmed tMSPs, the high mode is dominant, typically occurring for 65–80% of the time (e.g., Archibald et al. 2015; Bogdanov et al. 2015; de Martino et al. 2013; Papitto et al. 2013; Linares 2014). A similar behavior has been observed in the candidate tMSP CXOU J110926.4−650224 (Coti Zelati et al. 2024). In contrast, some other candidates spend a smaller fraction of time in the high mode, such as around 45–50% for 3FGL J1544.6−1125 (hereafter J1544; Gusinskaia et al. 2025) and below 40% for 4FGL J0407.7−5702 (Miller et al. 2020), while 3FGL J0427.9−6704 seems to spend most of its time in an X-ray flaring mode (Li et al. 2020). However, continued long-term monitoring is needed to better establish the long-term behavior of these systems.

The archetype of tMSPs, PSR J1023+0038 (hereafter J1023; Archibald et al. 2009) has been continuously observed in all energy ranges since its turn-on from a long rotation-powered state. Its X-ray emission switches from high to low mode on a timescale of ≲10 s, whereas the switch from low to high mode takes ≈30 s (Baglio et al. 2023). Flickering and flaring activities are also observed across UV, optical, and near-infrared (NIR) bands. The UV mode switching, correlating with X-rays, was clearly identified through high-time resolution observations performed with the Hubble Space Telescope (HST; Jaodand et al. 2021a; Miraval Zanon et al. 2022; Baglio et al. 2023). High-cadence optical observations (Shahbaz et al. 2015, 2018; Hakala & Kajava 2018) have revealed sharp-edged, rectangular flat-bottomed dips strongly resembling X-ray bimodal behavior (see also Kennedy et al. 2018; Papitto et al. 2018). On the other hand, evidence of NIR bimodal variability is limited. Dips in this band appear less pronounced than those at higher energies, possibly due to the enhanced contribution from the outer disk regions and the irradiated companion star (Baglio et al. 2019, 2023; Papitto et al. 2019). Simultaneous radio and X-ray observations have shown an anti-correlated variability pattern (Bogdanov et al. 2018). Specifically, when the source switched from the high to the low X-ray mode, a brightening in radio emission was observed, and vice versa. Thisbehavior has been interpreted as being due to rapid ejections of plasmoids during switches from high to low modes (Bogdanov et al. 2018; see also Baglio et al. 2023). X-ray, UV, and optical pulsations at the NS spin period have been detected exclusively during the high mode (Archibald et al. 2015; Ambrosino et al. 2017; Papitto et al. 2019; Jaodand et al. 2021a; Miraval Zanon et al. 2022). This, along with the similar pulse shapes and the fact that the spectral energy distribution (SED) of the pulsed emission from optical to X-ray wavelengths follows a single power-law relation (Papitto et al. 2019; Miraval Zanon et al. 2022), suggests a common underlying mechanism for the optical/UV and X-ray pulsations.

Recently, the so-called mini-pulsar nebula scenario was proposed, and it suggests that during the high mode, X-ray/UV/optical pulsations originate from synchrotron radiation in the boundary region where the particle wind ejected from a rotation-powered pulsar meets the matter from the inner accretion disk at about 100 km from the pulsar (Papitto et al. 2019; Veledina et al. 2019; see also Illiano et al. 2023). Recent IXPE polarimetric observations have revealed polarized X-ray emission with a polarization degree of up to 15% during the high mode, strongly supporting the rotation-powered pulsar wind model (Baglio et al. 2025). The switch to the low mode is believed to be triggered by discrete mass ejections (Baglio et al. 2023), which remove material from the inner disk and cause a decline in the X-ray, UV, and optical flux. As the flow from the disk refills the inner regions, the boundary region front is restored, and the system switches back to the high mode.

The source J1544, initially identified as a γ-ray source in the third Fermi-LAT source catalog (Acero et al. 2015), has emerged as a robust candidate tMSP in the sub-luminous disk state due to its X-ray and γ-ray emission properties being strikingly similar to those of J1023 (Bogdanov & Halpern 2015). Specifically, the X-ray emission of J1544 exhibits a clear bimodal behavior, alternating between high- and low-intensity modes (Bogdanov & Halpern 2015; Bogdanov 2016), and it shows a clear dipping behavior in the optical band (Bogdanov & Halpern 2015). The spectrum of the X-ray emission is well described by an absorbed power law with a photon index of Γ ∼ 1.7 (Bogdanov 2016), which is consistent with the typical spectra observed for tMSPs during their intermediate states (e.g., Papitto & de Martino 2022). The broadband SED from the optical band to the γ-rays closely resembles that of J1023 (Fig. 4. from Bogdanov 2016). Further evidence for classifying J1544 as a tMSP comes from its position on the radio/X-ray luminosity plane (for an X-ray luminosity of 5.51 × 10³³ erg s⁻¹; Jaodand et al. 2021b), which is very close to the region occupied by the confirmed tMSPs J1023 and XSS J12270−4859. Additionally, J1544 shows variable radio continuum emission on timescales of a few days, akin to the behavior observed in J1023 (Jaodand et al. 2021b; Gusinskaia et al. 2025). Hints of an anti-correlated variability pattern between simultaneous radio and X-ray observations resembling that seen in J1023 have been reported by Gusinskaia et al. (2025). Optical spectroscopy has revealed that J1544 is nearly face-on and in a ∼5.8 h orbit around a likely mid-to-late K-type main-sequence star (Britt et al. 2017). Using graph theory algorithms, García et al. (2024) recently showed that given the system parameters and their associated uncertainties, J1544 predominantly clusters with confirmed tMSPs, thereby reinforcing its classification as a highly promising candidate.

Ultimately, to validate J1544 as a tMSP, it is necessary to observe a transition to either a rotation-powered or accretion-powered state with detectable millisecond pulsations. At the time of writing, IGR J18245−2452 is the only tMSP to have exhibited a bright accretion outburst, during which 3.9-ms X-ray pulsations were detected (Papitto et al. 2013). This source was previously identified as the rotation-powered radio pulsar PSR J1824−2452I in the M28 globular cluster (Manchester et al. 2005). While IGR J18245−2452 was observed during an X-ray outburst, where the high X-ray luminosity (L_X ∼ 10³⁶ − 10³⁷ erg s⁻¹) enabled the detection of coherent millisecond pulsations, such pulsed signals have never been blindly discovered in the sub-luminous disk state. The X-ray luminosity in this state (L_X ∼ 10³³ − 10³⁴ erg s⁻¹), which is much weaker than that of accreting MSPs in outburst, hinders the discovery of signals at the NS spin period. The detection of X-ray pulsations in the sub-luminous disk state for J1023 (Archibald et al. 2015) and XSS J12270−4859 (Papitto et al. 2015) effectively relied on prior knowledge of the pulsar ephemeris obtained during their rotation-powered states (Archibald et al. 2009; Roy et al. 2015).

In light of the prevailing scarcity of multi-wavelength observations of J1544, we coordinated the most extensive high-time resolution campaign ever conducted on this source, covering the radio band to the X-rays. The goal was to investigate the properties of J1544 in detail across the electromagnetic spectrum, connecting spectral variability and mode switching phenomena across multiple bands to provide a more comprehensive understanding of the source behavior.

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 Sects. 4 and 5, we discuss our findings and draw conclusions, respectively.

2. Observations and data analysis

Table 1 provides a list of the observations performed during the four-day-long multi-wavelength campaign. The first day (Day 1) involved XMM-Newton, the fast optical photometer SiFAP2 mounted on the 3.6-m Telescopio Nazionale Galileo (TNG), the HST, the Rapid Eye Mount (REM) telescope, and the Karl G. Jansky Very Large Array (VLA); the second and third days (Day 2–3) involved the Nuclear Spectroscopic Telescope Array (NuSTAR) and REM; the fourth day (Day 4) involved the Neutron Star Interior Composition Explorer Mission (NICER), the high-speed, quintuple-beam camera HiPERCAM mounted on the 10.4-m Gran Telescopio Canarias (GTC), and REM. To further characterize the source across different wavelengths, we also included observations acquired with the Near Infrared Camera Spectrometer (NICS) at the TNG in 2016, as well as additional radio data collected in November and December 2024 using the Australia Telescope Compact Array (ATCA). In the following, we detail these observations and outline the procedures used for data processing and analysis. Archival NIR observations are described in Sect. 2.9, while additional radio data acquired in November and December 2024 are detailed in Sect. 2.10.

2.1. XMM-Newton (Day 1)

We analyzed XMM-Newton (Jansen et al. 2001; Schartel et al. 2022) Target of Opportunity (ToO) observations of J1544 performed on 2024 February 8 (AO-22 Proposal 092346, PI: Miraval Zanon), using the European Photon Imaging Cameras (EPIC) and the Optical/UV Monitor telescope (OM; Mason et al. 2001). The EPIC-pn (Strüder et al. 2001) operated with a time resolution of 29.5 μs (fast timing mode), the two EPIC-MOS (Turner et al. 2001) with a time resolution of 0.3 s (small window mode), and the OM with a time resolution of 0.5 s (fast window mode) with the V filter (effective wavelength 5407 Å; full-width at half-maximum 684 Å). We processed and analyzed the Observation Data Files using the Science Analysis Software (SAS; v.21.0.0). We identified high background flaring activity in the 10–12 keV light curves of the whole field of view over a time span of ∼4 ks during the initial part of the observation. However, we checked that these flares did not impact the identification of X-ray modes. Hence we decided to keep all the data for the following analyses. We extracted the spectra as described below, either by discarding the high-background contaminated data or by including it, and verified that the results were entirely consistent with each other.

For the EPIC-pn we defined the source region as a 11-pixel-wide strip centered on the brightest pixel column (RAWX = 32–42), and the background region as a strip with the same width far from the source position (RAWX = 3–13). For each MOS, we extracted source photons within a circular region centered on the source position with a 40″ radius, and background photons from a 80″ wide, source-free circular region on one of the outer CCDs. Using the task epiclccorr, we extracted the background-subtracted light curves from the three EPIC instruments over the time interval of simultaneous coverage, binning it with a time resolution of 20 s. We selected the high modes as the time intervals when the count rate exceeded 0.8 counts s⁻¹, and the low modes as those when the count rate dropped below 0.4 counts s⁻¹ (see Fig. 1). No flaring mode was observed. The good time intervals corresponding to the high and low modes were then used to extract background-subtracted spectra, response matrices, and ancillary files separately for each mode. The average EPIC-pn spectrum was extracted with a minimum of 100 counts in each channel within the 0.7–7.0 keV energy range since EPIC-pn data in timing mode are not well calibrated below 0.7 keV¹ along with the background domination at energies above 7 keV. For the average EPIC-MOS1 and EPIC-MOS2 spectra extracted in the 0.3–10 keV band, a minimum of 50 counts per channel bin was selected. During high-mode intervals, EPIC spectra were extracted with a minimum of 50 counts per channel bin. Conversely, during low-mode intervals, where the EPIC-pn spectrum is background-dominated, only EPIC-MOS1 and EPIC-MOS2 spectra with a minimum of 20 counts per channel wereconsidered.

Fig. 1. Distribution of count rates obtained from the background-subtracted XMM-Newton/EPIC light curve acquired on 2024 February 8 binned with a time resolution of 20 s. We defined the high modes as the time intervals when the count rate exceeded 0.8 counts s−1 (red dashed line), and the low modes as those when the count rate dropped below 0.4 counts s−1 (blue dashed line). A “transition” region is also present, following the procedure adopted by Bogdanov et al. (2015) for J1023. No flaring mode was observed.

We extracted the optical background-subtracted light curve using the omfchain pipeline with default parameters. The light curve was binned at a time resolution of 1200 s to search for any orbital modulation as previously done by Bogdanov & Halpern (2015). We converted the observed count rates, R_OM, to magnitudes in the Vega system using the relation² mag = 17.963–2.5 log₁₀(R_OM). Following the procedure described by Baglio et al. (2023) to align the X-ray and optical light curve with those obtained from other observations at different wavelengths, we converted the photon arrival times from the Terrestrial Time (TT) to the Coordinated Universal Time (UTC) without applying a barycentric correction³

2.2. TNG/SiFAP2 (Day 1)

The fast optical photometer SiFAP2 (Ghedina et al. 2018), mounted on the 3.58-meter INAF TNG (Barbieri et al. 1994), observed J1544 on 2024 February 8 (PI: Illiano), simultaneously with XMM-Newton and HST observations. SiFAP2 observations were performed in white light spanning the 320–900 nm band and peaking between 400 and 600 nm (i.e., roughly corresponding to the B and V Johnson filters; see Supplementary Figure 1 in Ambrosino et al. 2017). Due to adverse weather conditions and gusty winds, we only retained data between 60348.21 MJD and 60348.26 MJD, corresponding to a total duration of 4320 s. To take into account spurious sky effects and monitor the sky condition, we simultaneously observed a V = 12.3 mag reference star TYC 5605-459-1 (Høg et al. 2000). We estimated the sky background by offsetting the telescope in a region located 10″ away from the target (i.e., outside the 7″ × 7″ region observed when pointing at the source position) toward the east direction for 31 s at the start of the observation. The background rates are 17249.8 counts s⁻¹ for the source channel and 24568.3 counts s⁻¹ for the reference channel, the difference being due to the presence of clouds. During the retained time interval, the mean source count rate was 13837.5 counts s⁻¹, with an extrapolated background count rate of 12848.3 counts s⁻¹. For the reference channel, the mean count rate was 581086.8 counts s⁻¹, with an extrapolated background of ∼2.3 × 10⁴ counts s⁻¹. To extract a background-subtracted light curve, the background counts are subtracted from the source counts. This process often resulted in counts consistent with zero within the errors, indicating that the source was not detected.

2.3. HST (Day 1)

The Space Telescope Imaging Spectrograph (STIS; Woodgate et al. 1998) aboard HST observed J1544 on 2024 February 8 for ∼35 minutes (Program 17587, PI: Illiano). These observations were performed using the near-UV Multi-Anode Micro-channel Array detector (NUV-MAMA) in TIME-TAG mode (time resolution of 125 μs), and collected with the G230L grating equipped with a 52″ × 0.2″ slit with a spectral resolution of ∼500 over the nominal range of 1570–3180 Å. Due to high background noise in the spectrum outside the 2000–3000 Å range, we restricted our spectral and timing analysis to this interval. However, we verified that the light curve does not change significantly when comparing the nominal and chosen ranges. We calibrate the UV data using the CALSTIS v3.4.2 pipeline, which includes dark subtraction and flat fielding. To extract the background-subtracted light curve, we employed the stistools Python package’s inttag function⁴ to obtain 20-s-long time series before applying CALSTIS.

2.4. REM (Day 1–2–3–4)

The REM telescope located in La Silla, Chile, observed J1544 over four consecutive nights: 2024 February 8, 9, 10, and 11 (MJD 60348–60351; PI: Baglio). For each night, three sets of observations were acquired using the NIR camera with the H-band filter (central wavelength: 1.662 μm). Each set comprised five dithered 30-second integration exposures, combined to enhance background subtraction. The target was not bright enough to be detected in individual exposure images. Additionally, seeing conditions were suboptimal, with the PSF varying between 2″ and 4″ across different epochs. To attempt a detection, we averaged all images for each epoch. After averaging, the target was consistently detected, except for the final epoch on February 11, when the seeing conditions were at their worst (∼4″). We conducted aperture photometry using the daophot tool (Stetson 1987), and an aperture approximately twice the average full-width at half-maximum (FWHM) of the flux distribution of field stars. Magnitudes were then calibrated against a group of four bright stars in the field, with magnitudes reported in the Two-Micron-All-Sky-Survey (2MASS; Skrutskie et al. 2006)⁵ catalog.

2.5. VLA (Day 1)

The Karl G. Jansky VLA observed J1544 (project code: SX222346, PI: Miraval Zanon) on 2024 February 8 for a total of 4 hours at the central frequency of 6 GHz (C-band) with a bandwidth of 4 GHz. The target and the phase calibrator J1543−0757 were observed in ten-minute cycles, with eight minutes on the former and two minutes on the latter. The distance between the target and the phase calibrator is about 3.5°. The observation included a scan on the flux and bandpass calibrator (J1331+3030). The total on-source time is ∼2.89 h (specifically 10395 s), corresponding to ∼72% of the total observation time, while the remaining 28% was allocated to calibrators and overheads. Data were calibrated using the custom casa pipeline (Version 6.5.4, CASA Team et al. 2022) and visually inspected for possible radio frequency interference (RFI). The final images were produced with the tclean task in casa (Version 6.5.4), using Briggs weighting with robust parameter 0.5.

2.6. NuSTAR (Day 2–3)

NuSTAR (Harrison et al. 2013) observed J1544 on 2024 February 9, for a total elapsed time of ∼94.8 ks and a net exposure of ∼45 ks (as part of the joint NuSTAR/XMM-Newton cycle-22 program 092346, PI: Miraval Zanon)⁶. Data were processed and analyzed using NuSTARDAS v2.1.2 along with the calibration database (CALDB) 20240506. We accounted for South Atlantic Anomaly passages using nupipeline task with SAAMODE = optimize and TENTACLE = yes. Source and background events were extracted from circular regions of 60″ radius centered at the source position and in a source-free area of the CCD, respectively. J1544 was detected up to energies of ∼30 keV in the two focal plane modules (FPMA and FPMB). Light curves were extracted separately for the two modules, combined to increase the signal-to-noise ratio (S/N), and binned at a time resolution of 100 s. Spectra for both detectors and the corresponding response files were extracted using the nuproducts command. We grouped the average spectra with a minimum of 50 counts per channel. Spectra were also extracted separately during high-mode intervals, ensuring a minimum of 20 counts per channel. For the low-mode intervals, the counting statistics was insufficient to extract meaningful spectra.

2.7. NICER (Day 4)

NICER (Gendreau et al. 2012) observed J1544 on 2024 February 11 (ObsID 6030120103, PI: Illiano). We analyzed the data using HEASoft version 6.33.2 and NICERDAS version 12, with CALDB version 20240206. After retrieving the latest geomagnetic data with the nigeodown task, we processed the observational data using nicerl2 task with default screening parameters. The obtained cleaned event files were used to generate the light curve in the 0.5–10 keV energy range with a time bin of 10 s using the nicerl3-lc pipeline with SCORPEON background model⁷, version 23. Due to visibility issues, only 240 s of NICER observations were conducted simultaneously with GTC/HiPERCAM, and no clear low-mode intervals were detected during this limited timeframe. Consequently, we opted not to include the NICER data in the subsequent analysis.

2.8. GTC/HiPERCAM (Day 4)

The high-speed, quintuple-beam camera HiPERCAM (Dhillon et al. 2021) mounted on the 10.4-m Gran Telescopio Canarias (GTC) on the island of La Palma (Canary Islands) observed J1544 for ≃2 h on 2024 February 11 (PI: Coti Zelati). HiPERCAM was operated in windowed mode to achieve single frame exposure times of ≃0.19 s for the Super SDSS g_S, r_S, i_S, and z_S filters and three times longer in the u_S filter (to compensate for the lower filter throughput). The average seeing during the observations was ≃1.7″.

We processed and analyzed the data using the HiPERCAM reduction software v. 1.5.2⁸. We processed each science frame by subtracting the bias and applying a flat-field correction using the median of twilight-sky frames. Additionally, we corrected for potential fringing effects in the z_S filter by using publicly available, pre-prepared fringe maps. We used variable aperture photometry to extract counts from the target and a bright comparison star, that is, for each frame we adjusted the aperture radius based on the FWHM of the fitted profile of the comparison star. The comparison star was Gaia DR3 6268482198463636224 (also referred to as USNO-B1.0 0785-0287404). We estimated the sky level around the target and the comparison star by averaging the counts within an annular region surrounding each of them. Next, we performed differential photometry by dividing the background-subtracted counts of J1544 by those of the comparison star. Finally, the time series were rebinned to a time bin of 20 s to increase the S/N while still allowing us to study the variability of the optical emission over timescales compatible with those expected for the X-ray mode switching.

2.9. TNG/NICS (2016)

To further characterize J1544 across different wavelengths, we included observations acquired with NICS (Baffa et al. 2001) at the TNG. NICS observations were conducted during two nights on 2016 March 29 and 30 (Prg. ID: A33DDT4A; PI: de Martino) for a total exposure of ∼4.8 h and ∼4.9 h respectively, under good seeing conditions. The J, H, and K filters were used cyclically, with 20-s exposures at 10 mosaic positions. The image reduction was performed using the jitter task of the ESO-eclipse software package⁹. Aperture photometry was carried out using the IRAF DAOPHOT package (Stetson 1987). The magnitudes were calibrated with the two nearby stable reference stars listed in the 2MASS catalog, 2MASS15444507-11284959 and 2MASS15443354-112930.

2.10. ATCA (Nov./Dec. 2024)

We obtained additional radio observations of J1544 with the Australia Telescope Compact Array (ATCA) under program CX588 (PI: Carotenuto). The first observation was conducted on 2024 November 25 from 01:40 UT to 07:00 UT, with the telescope in the 750D configuration. A second observation was carried out from 2024 December 3 at 20:40 UT to 2024 December 4 at 02:15 UT, with the telescope in the H214 configuration. For both epochs, data were recorded simultaneously at central frequencies of 5.5 GHz and 9.0 GHz, with 2 GHz of bandwidth at each frequency. We used PKS B1934−638 for the bandpass and flux density calibration, and PKS B1540−077 for the complex gain calibration. Data were flagged, calibrated, and imaged using standard procedures within CASA. When imaging, we included antenna CA06 for both epochs, and we used a natural weighting scheme to minimize the root mean square (RMS) noise of the image.

3. Results

3.1. X-ray emission

The XMM-Newton/EPIC background-subtracted light curve extracted during Day 1 revealed hundreds of low-mode intervals (see Fig. 2). By applying our threshold to discern different intensity levels in the EPIC net count rates (see Sect. 2.1), we determined that J1544 spent ≃58% and ≃33% of the observation time in the high and low modes, respectively (similar proportions were reported by Gusinskaia et al. 2025), and the remaining time switching between these modes. The shortest observed duration of a low-mode episode spanned ∼20 s (corresponding to a single time bin), whereas the longest recorded duration extended to ∼520 s (occurring between ∼60348.26 MJD and ∼60348.27 MJD; see Fig. 3).

Fig. 2. Temporal evolution of the X-ray, UV, and optical emissions of J1544 during the first three days of observations. For Day 1, the light curves are shown in decreasing order of energy band, from top to bottom, with the XMM-Newton/EPIC (20 s time bin), HST/STIS (20 s), and XMM-Newton/OM (1200 s). For Days 2–3, the 3–30 keV NuSTAR light curve (100 s time bin) is shown. The yellow-shaded areas denote the time intervals of the low X-ray mode identified by XMM-Newton/EPIC on Day 1 and by NuSTAR on Day 2–3. The error bars represent 1σ uncertainties.

Fig. 3. Zoom-in of the time series collected on Day 1 capturing the period of maximum overlap among different telescopes.

The NuSTAR light curve displays a lower S/N compared to the XMM-Newton/EPIC light curve (see Fig. 2). However, we identified high- and low-mode intervals and used the good time intervals to extract background-subtracted spectra, response matrices, and ancillary files separately for each mode.

We performed the X-ray spectral analysis on XMM-Newton and NuSTAR data (see Sects. 2.1 and 2.6) using the X-ray spectral fitting package XSPEC (Arnaud 1996) version 12.14.0. We adopted the interstellar medium abundance and the cross-section tables from Wilms et al. (2000) and Verner et al. (1996), respectively. Uncertainties on spectral parameters are given at 1σ confidence intervals unless otherwise stated. The average broadband spectrum is well described by an absorbed power-law model (constant * TBabs * powerlaw). We included a renormalization factor in our spectral modeling to address cross-calibration uncertainties between the two X-ray telescopes. We verified that these factors were consistent within 10% for both the average broadband spectrum and the spectra extracted separately in the high- and low-mode intervals. This condition was satisfied in all cases except for NuSTAR in the average spectrum, where the constant deviated by 11%. The discrepancy can be attributed to the NuSTAR observations being conducted one day after the XMM-Newton observations, resulting in differences in the number and duration of the high and low modes and thus the average intensity level. Supporting this interpretation, the renormalization factors for the spectra extracted separately in the high- and low-mode intervals remained consistent within 10%. The best-fitting parameters are an absorption column density of N_H = (1.52 ± 0.06)×10²¹ cm⁻², a photon index of Γ = 1.63 ± 0.01, and an unabsorbed 0.3–10 keV flux of F_unabs = (3.54 ± 0.02)×10⁻¹² erg cm⁻² s⁻¹ (χ² = 445.7 for 414 degrees of freedom, d.o.f.). The estimated value of N_H is broadly consistent with the average absorption column density in the direction of the source of ∼1.2 × 10²¹ cm⁻² (HI4PI Collaboration 2016), and both N_H and Γ are consistent within 3σ with values previously reported for this source from X-ray observations by Bogdanov & Halpern (2015), Bogdanov (2016), and Gusinskaia et al. (2025). Given that the reflection spectrum, and in particular the iron Kα complex, is often observed in LMXBs hosting a NS (see, e.g., Di Salvo et al. 2023, and references therein), we searched for similar features in J1544 by adding Gaussian terms to the power-law model. This approach is analogous to that used for the search in other tMSPs or candidate systems in the sub-luminous disk state. The absorption column density and the photon index of the power law were fixed to their best-fit values from the previous analysis. The line energies were set sequentially at 6.4 keV (for neutral or lightly ionized Fe I Kα₁ and Kα₂ transitions), then 6.7 keV (He-like Fe XXV), and 6.97 keV (H-like Fe XXVI). To simulate lines narrower than the instruments’ intrinsic energy resolution, we set the Gaussian line width to zero. No significant emission feature was detected. As a result, we determined 3σ upper limits for the equivalent widths of any Gaussian features at 6.4, 6.7, and 6.97 keV to be 20, 23, and 32 eV, respectively. The iron line has not been detected in this system or any other confirmed or candidate tMSPs in the sub-luminous disk state, including J1023 (see, e.g., Coti Zelati et al. 2014, 2019) and XSS J12270−4859 (de Martino et al. 2010), providing no evidence of fluorescence or disk reflection in these systems.

When fitting the spectra extracted separately in the high- and low-mode intervals, we fixed the N_H to the value obtained from the average spectral fitting. The high-mode spectrum includes all XMM-Newton and NuSTAR data, whereas the low-mode spectrum comprises only the two MOS data due to background dominance in the EPIC-pn and low statistics in NuSTAR data. For the high mode, we obtained Γ_H = 1.623 ± 0.008, and unabsorbed 0.3–10 keV flux of F_unabs, H = (5.58 ± 0.04)×10⁻¹² erg cm⁻² s⁻¹ (χ²/d.o.f. = 521.5/538), while for the low mode Γ_L = 1.65 ± 0.07 and an unabsorbed 0.3–10 keV flux F_unabs, L = (0.38 ± 0.02)×10⁻¹² erg cm⁻² s⁻¹ (χ²/d.o.f. = 37.9/45). The values for the photon indices in the different modes are compatible with each other and with the average value within 1σ. To improve the statistics and search for a potential steepening of the photon index in the low mode (Campana et al. 2019), we added all XMM-Newton archival observations performed on this source at the moment of writing. The results remained consistent with those presented above, with a slight improvement in parameter precision, as detailed in Appendix A.

3.2. UV emission

We conducted HST observations on Day 1 concurrently with XMM-Newton to investigate the presence of a bimodal pattern in the UV data. Notably, similar switches between high and low modes have been recorded in HST data of J1023 (Miraval Zanon et al. 2022; Jaodand et al. 2021a; Baglio et al. 2023). The UV light curve of J1544, shown in Fig. 3 (dark blue points), exhibits substantial variability, with some instances appearing to align with the varying intensity levels observed in X-rays, particularly between ∼60348.25 and ∼60348.26 MJD. However, the background level (light blue points in Fig. 3) varies in amplitude and timescale similarly to the total count rate, particularly in the final part of the HST observation. As a result, the expected high/low mode switches, likely weaker than this variability, may be obscured by background contamination. Future HST observations may help determine whether this background variability is persistent and, despite the faint nature of J1544, potentially enable us to identify bimodality within the count rate distribution.

3.3. Optical variability and its connection with mode switching

Attempts to bin the XMM-Newton/OM light curve acquired on Day 1 using time bins of a few hundred seconds were unsuccessful in detecting mode switching. As noted by Bogdanov & Halpern (2015), the faint nature of J1544 makes it challenging to identify variability on timescales of a few tens of seconds with the OM. Even for J1023, variability on such short timescales was not observed, likely masked by the intrinsic statistical uncertainties of the data (Baglio et al. 2019). Consequently, we binned the OM light curve using a time bin of 1200 s (see Fig. 2) to increase the S/N. We observed similar fluctuations in the V magnitude as reported by Bogdanov & Halpern (2015) in the U filter, recording a modulation semi-amplitude of ∼0.4 mag. To quantify the object’s inherent variability beyond the limitations imposed by observational uncertainties, we computed the excess variance (Nandra et al. 1997; Edelson et al. 2002; Vaughan et al. 2003; Yuk et al. 2022), which resulted to be of ( − 0.02 ± 0.03) mag². This value is consistent with zero within the uncertainty and suggests that the observed variability is not intrinsic to the source. The U filter data acquired in 2014 showed a sinusoidal pattern with a periodicity of ∼5.2 h (Bogdanov & Halpern 2015), differing from the orbital period of ∼5.8 h estimated through spectroscopic observations by Britt et al. (2017) using radial velocity curves. To further investigate this discrepancy and compare with data in the V filter, we computed the excess variance also on the U filter data, obtaining (0.01 ± 0.05) mag². This value is also consistent with zero within the uncertainty, indicating that we cannot draw definitive conclusions.

GTC/HiPERCAM light curves acquired on Day 4 displayed flickering and dipping activities (Fig. 4) reminiscent of the X-ray bimodal pattern. This optical behavior resembled that observed for the prototype tMSP J1023 (Shahbaz et al. 2015, 2018; Hakala & Kajava 2018) as well as for the candidate CXOU J110926.4−650224 (Coti Zelati et al. 2024). Coti Zelati et al. (2024) recently reported that the optical emission from the candidate tMSP CXOU J110926.4−650224 becomes redder during the low mode compared to the high mode, which aligns with the expectation of the proposed model that explains optical and X-ray pulsations from J1023 as a result of synchrotron radiation at the boundary region between the striped pulsar wind and the inner accretion disk (Papitto et al. 2019; Veledina et al. 2019). The findings from the extensive multi-wavelength campaign conducted by Baglio et al. (2023) on J1023 supported the hypothesis that during switches from high to low modes, the inner and hotter regions of the accretion flow are ejected, resulting in a decline of X-ray, UV and optical fluxes. This leaves a fainter optical emission from the cooler outer regions of the accretion disk. Taking into account these recent observations, along with the presence of at least two distinct rectangular flat bottom dips highlighted in Fig. 4, we plotted the magnitudes in the g_S filter against the magnitude colors g_S − r_S, g_S − i_S, g_S − z_S (see Fig. 5). To quantify the observed correlation, we rebinned our data to reduce the scatter and computed Spearman’s rank correlation coefficients, which were found to be ≃0.77, 0.78, and 0.93, respectively, for 32 d.o.f. The corresponding probabilities of observing such coefficients from uncorrelated data are ≃8.8 × 10⁻⁸, 1.3 × 10⁻⁸, and 1.9 × 10⁻¹⁶, respectively. These results indicate that the optical emission becomes significantly redder during the low mode compared to the high mode, particularly when measured using the g_S and z_S filters, which are separated by the largest wavelength intervals.

Fig. 4. GTC/HiPERCAM light curves acquired on Day 4 in the Super SDSS uS, gS, rS, iS, and zS filters. The yellow-shaded areas highlight the most prominent potential low modes. The error bars represent 1σ uncertainties.

Fig. 5. Color-magnitude diagrams for J1544 derived from the most highly sampled gS, rS, iS, and zS bands of the GTC/HiPERCAM time series (Fig. 4) and rebinned to emphasize correlation trends. The upper panels display the smoothed Kernel Density Estimation curves that illustrate the distribution of the color-magnitude values.

3.4. Hint of an optical flare

While analyzing archival XMM-Newton data, we identified peculiar features in the OM light curve acquired on 2018 January 28 (ObsID 0800280101). The OM was operated in fast window mode with a time resolution of 0.5 s in white light. We analyzed the data as detailed in Sect. 2.1. As shown in Fig. 6, the OM light curve reveals an increase in flux following an observational gap around 58147.1 MJD and a distinct flare-like feature around 58147.4 MJD. These findings are particularly intriguing, as no flares have previously been reported from J1544, which is peculiar given that confirmed and candidate tMSPs frequently exhibit flaring activity (e.g., Tendulkar et al. 2014; Bogdanov et al. 2015; Coti Zelati et al. 2024; Li et al. 2020). If confirmed, this would represent the first detection of a flare from this source. However, caution is warranted due to the lack of a corresponding feature in the X-ray light curve, as shown in the bottom panel of Fig. B.1 provides a detailed analysis to verify that the observed flare-like event is not spurious or instrumental and to explore its possible origin.

Fig. 6. XMM-Newton/OM light curve acquired in 2018 using 100-s bins. The error bars represent 1σ uncertainties.

3.5. Near-infrared variability

Over the four days of this multi-wavelength campaign, we performed NIR observations using the REM telescope with the H-band filter. As outlined in Sect. 2.4, due to the source faintness and suboptimal seeing conditions, we were only able to detect J1544 by averaging all exposure images for each observational epoch. On the fourth night, the poor seeing conditions prevented any detection. The H-band magnitude of J1544 was found to vary to some extent between epochs. Specifically, we obtained H = (16.86 ± 0.26) mag, H = (16.34 ± 0.19) mag, H = (16.73 ± 0.26) mag in the first three epochs of observations, respectively, while we determined a 3σ upper limit of > 16.88 mag during the last epoch. Although the REM observations do not allow for study of the short-term variability on timescales of a few tenths of a second, the reported average magnitudes still indicate that the NIR emission is variable.

To further enhance the multi-band characterization of J1544, we analyzed NIR observations acquired in 2016 with TNG/NICS. Due to limited X-ray coverage, we cannot definitively confirm that the source was in the intermediate sub-luminous disk state during our NIR observations. However, X-ray flux estimates derived from Swift/XRT observations on 2015 October 2 and 2016 May 7 (Obs. IDs: 00084762008 and 00092235001; ≃(3 − 4)×10⁻¹² erg cm⁻² s⁻¹ in the 1 − 10 keV energy range) align with the values reported by Jaodand et al. (2021b) and fall within the expected range for the sub-luminous disk state of tMSPs. Although state transitions can occur on timescales as short as a few weeks (e.g., Papitto et al. 2013), the sub-luminous disk state in tMSPs like J1544 has shown remarkable stability over periods exceeding a decade (e.g., Stappers et al. 2014; Baglio et al. 2023). This sustained stability suggests that J1544 has likely remained in this state since its discovery.

Figure 7 shows the light curves alternately acquired in the J, H, and K filters. The NIR time series exhibit a variability that may resemble the X-ray bimodal pattern. To highlight this, we plotted the differential magnitudes by subtracting the weighted average value in each filter for each night. The weighted average magnitudes of J1544 during the first night were J = (17.083 ± 0.005) mag, H = (16.659 ± 0.006) mag, and K = (16.32 ± 0.01) mag, and during the second night J = (17.100 ± 0.005) mag, H = (16.644 ± 0.006) mag, and K = (16.29 ± 0.01) mag. The consistency between these average magnitudes in the H band and those observed by REM during our 2024 multi-wavelength campaign further supports that J1544 was in the sub-luminous disk state during the TNG/NICS observations in 2016. However, the lack of simultaneous X-ray coverage prevented us from confirming the presence of high and low modes in the NIR band. We computed the excess variance, as outlined in Sect. 3.3, to quantify the intrinsic variability beyond measurement uncertainties. The estimated values for the J, H, and K bands during the first night are (0.01 ± 0.02) mag², (0.03 ± 0.03) mag², and (0.04 ± 0.03) mag², respectively. During the second night of NIR observations, we obtained (0.01 ± 0.02) mag², (0.02 ± 0.02) mag², and (0.07 ± 0.04) mag² for the J, H, and K bands, respectively. While these results indicate that we cannot draw definitive conclusions about the intrinsic variability of J1544 in the NIR band, potentially related to the bimodality between high and low modes, visual inspection of the light curves in Fig. 7 and the slightly higher excess variance in the K band may suggest that further investigations are required to better characterize the variability pattern at these wavelengths. To our knowledge, these are the first NIR light curves reported for this source, and future high-time resolution NIR observations with simultaneous X-ray coverage could provide valuable insights.

Fig. 7. Near-infrared light curves of J1544 observed with TNG/NICS on 2016 March 29 (top panel) and 30 (bottom panel). We plotted differential magnitudes by subtracting the weighted average value in each filter for each night (see Sect. 3.5 for details). Photometric data were acquired in the J (blue points), H (green points), and K (red points) filters. The error bars represent 1σ uncertainties.

3.6. Radio variability

J1544 appeared pointlike in all our radio observations, with an angular resolution ranging from 3.5″ to 6″. The source was not detected during our VLA observation on 2024 February 8, with a 3σ flux density upper limit of ∼8 μJy at 6 GHz. While we observed an excess of emission at the target’s position, the low signal-to-noise ratio (∼3) and the presence of nearby noise peaks with comparable flux densities prevent us from confirming a detection. Considering the radio/X-ray anti-correlation observed for J1023 – where the radio flux increases as the system switches from the high to the low X-ray mode – (Bogdanov et al. 2018; Baglio et al. 2023) and the hint of a similar trend for J1544 (Gusinskaia et al. 2025), we performed a deeper analysis separately extracting radio images during X-ray high and low modes simultaneously observed with XMM-Newton. Considering the total VLA on-source time of ∼2.89 h (see Sect. 2.5) and that J1544 spent ∼58% and ∼33% of the XMM-Newton observation time in the high and low modes, respectively (see Sect. 3.1), we estimated that the radio images extracted during simultaneous X-ray high modes correspond to ∼1.68 h, while those extracted during low X-ray modes correspond to ∼0.95 h. However, the source remained undetected in both modes, with a 3σ flux density upper limit of 13 μJy during the X-ray low modes and 11 μJy during the high modes. Notably, during the X-ray low modes, we observed an excess emission with a signal-to-noise ratio of ∼4 (peak of ∼17 μJy/b). However, a definitive detection remains challenging due to the previously mentioned limitations, particularly the presence of noise peaks with comparable flux densities. For comparison, Gusinskaia et al. (2025) reported a VLA flux density of (28.6 ± 2.1) μJy during a 2019 observation in the C-band (4–8 GHz), which increased to (56.6 ± 8.3) μJy during simultaneous X-ray low-mode intervals. However, the total on-source time during their VLA observation was ∼4.2 h. Similarly to our findings, they also reported that J1544 spent ∼33% of the time in the X-ray low mode, meaning that their radio images extracted in the low mode corresponded to a longer total exposure time than ours.

J1544 was not detected in either of the ATCA observations, with a 3σ upper limit of 27 μJy at 7.25 GHz (combining the C and X bands) in the first epoch, and 21 μJy at 7.25 GHz in the second epoch. For both observations, the RMS noise was higher than expected due to the presence of strong RFI in our spectral windows. Stacking together the two ATCA epochs yielded a 3σ flux density upper limit of 15 μJy at 7.25 GHz, which is slightly higher than the VLA one. Figure 8 extends Fig. 1 from Jaodand et al. (2021b), integrating results from Gusinskaia et al. (2025) and this work, to illustrate the radio variability of J1544 against a relatively stable X-ray flux over the years.

Fig. 8. Radio variability of J1544 compared to its relatively stable X-ray flux over the years. This figure extends Fig. 1 from Jaodand et al. (2021b) by incorporating observations from Gusinskaia et al. (2025) and this work. Top panel: Unabsorbed X-ray flux in the 1–10 keV range (estimated using the entire observations, including both high and low modes) from Jaodand et al. (2021b) (dark blue squares), Gusinskaia et al. (2025) (light blue hexagons), and this work (green point). For consistency, we extracted the flux in the 1–10 keV band from our XMM-Newton observation, obtaining FX = (2.79 ± 0.04)×10−12 erg cm−2 s−1, by following the procedure detailed in Sect. 3.1. The low X-ray flux around MJD 57179 is likely due to the source remaining in the X-ray low mode throughout this short (418 s) observation, as noted by Jaodand et al. (2021b). Bottom panel: VLA radio flux density in the X band (8–12 GHz) from Jaodand et al. (2021b) (dark red squares are for detections and downward triangles are for 3σ upper limits) and from Gusinskaia et al. (2025) at 8–12 GHz (X band) around 58224 MJD (red hexagon) and at 4–8 GHz (C band) around 58491 MJD (orange hexagon). The 3σ upper limit of ∼9 μJy at 4–8 GHz (C band) from this work is indicated by the brown downward triangle. For consistency with Jaodand et al. (2021b), the plot shows the radio flux density from Gusinskaia et al. (2025) and the upper limit derived in this work both extracted using the entire VLA observation, rather than restricting to intervals simultaneous with X-ray low modes. The dashed vertical lines mark the epochs of X-ray observations for reference.

3.7. Spectral energy distribution

To derive the SED from the UV to the X-rays, we analyzed the HST NUV spectrum between 2000 and 3000 Å to avoid high background noise outside this range. We corrected for interstellar extinction using an absorption column density of N_H = (1.52 ± 0.06)×10²¹ cm⁻² as estimated from X-ray spectral fitting (see Sect. 3.1), and the empirical relation for the optical extinction A_V = N_H/[(2.87 ± 0.12)×10²¹ cm⁻²] (Foight et al. 2016). Using the commonly adopted Milky Way average value for the total-to-selective extinction ratio, R ≡ A_V/E(B − V) = 3.1 (e.g., Johnson 1965; Schultz & Wiemer 1975; Whittet & van Breda 1980; Fitzpatrick & Massa 1999), we determined a color excess of E(B − V) = (0.17 ± 0.1) mag. To account for interstellar reddening and correct the flux values across different frequencies, we employed the dust_extinction Python package with extinction curves from Gordon (2024) (see also Gordon et al. 2009, 2021; Fitzpatrick et al. 2019; Decleir et al. 2022).

The left panel of Fig. 9 shows the SED of the total emission from UV to X-rays for J1544, compared to J1023 as reported by Miraval Zanon et al. (2022). The striking similarity strongly supports the classification of J1544 as a very promising candidate tMSP in the sub-luminous disk state.

Fig. 9. Unabsorbed broadband SED of J1544. Left panel: SED of J1544 from UV to X-rays compared with that of J1023 from Miraval Zanon et al. (2022). For J1023 and J1544, HST data are shown in orange and red, XMM-Newton data are in light blue and dark blue, and NuSTAR data are in light green and dark green, respectively. The HST spectrum is plotted from 165 to 310 nm (i.e., ∼(1.0 − 1.8)×1015 Hz) for J1023 and from 200 to 300 nm (i.e., ∼(1.0 − 1.5)×1015 Hz) for J1544. The UV data were rebinned using the coronagraph Python package (Robinson et al. 2016; Lustig-Yaeger et al. 2019) with a low-resolution wavelength grid of width 4 for J1023 and width 20 for J1544. Right panel: The SED of J1544 from the optical band to X-rays, extracted during X-ray high modes observed with XMM-Newton, along with the best-fitting model (see text for details). Due to the lack of simultaneous X-ray observations, we include all available optical data from GTC/HiPERCAM. The red dashed, blue dotted, and orange solid lines represent the contributions from the accretion disk, companion star, and boundary region, respectively, while the green solid line shows their combined emission. The ratio between the data points and the best-fitting model is shown in the bottom panel.

We then extracted the SED only during the high modes (right panel of Fig. 9), where the source spends most of its time, ensuring better statistical quality compared to the low modes. This choice is further justified by the absence of a clear bimodal behavior in the UV light curve as well as by the challenges in modeling the SED of J1544 due to its poorly constrained parameters. We fitted the SED adopting a model based on the mini-pulsar nebula scenario (see Sect. 1 and Baglio et al. 2023 for a detailed discussion). In this framework, the X-ray emission primarily originates from synchrotron radiation produced at the boundary region between the pulsar wind and the inner accretion flow. In the SED extracted during the high modes of J1023 and analyzed by Baglio et al. (2023), an additional contribution to the X-ray emission was attributed to an optically thin synchrotron component from the base of a compact jet. However, in our analysis of J1544, we did not include this component as the available data do not provide sufficient constraints to model it reliably. Since this component contributed only ∼4% of the X-ray emission, its omission does not significantly impact our results. At lower frequencies, the model includes blackbody emission from the irradiated companion star and a multi-color blackbody component from the accretion disk. Following similar reasoning as for the optically thin synchrotron component, we did not include the optically thick synchrotron emission from the compact jet, which was considered for J1023 by Baglio et al. (2023). To better constrain the contributions from the accretion disk and companion star, we included all available optical data from GTC/HiPERCAM, as no strictly simultaneous X-ray and optical observations were available.

The free parameters in the model are the surface temperature of the companion star (T_*), the irradiation luminosity (L_irr), the inner radius of the accretion disk where optical/UV emission becomes relevant (r_in, opt/UV), the peak frequency of the synchrotron emission (ν_sync), the slope of the optically thin synchrotron emission component (α_sync), and the normalization of the synchrotron emission component (F_sync). The fixed parameters in the model include the source’s distance (D = 3.8 kpc; this value was estimated by Britt et al. 2017, but recent parallax measurements of the proposed optical counterpart by Koljonen & Linares 2023 suggest it should now be considered as an upper limit), the companion star radius (R_c = 0.75 R_⊙, assuming the upper limit on the companion mass of M_C = 0.7 M_⊙ from Britt et al. 2017 and the main–sequence mass–radius relationship), the albedo of the companion star (η_* = 0.1), and the binary separation (a). The latter is given by a = [G(M_NS + M_C)P_orb²/(4π)²]^1/3, where M_NS = 1.4 M_⊙ is assumed for the NS mass, P_orb = 5.8 h is the binary orbital period (Britt et al. 2017), and G is the gravitational constant. We also fixed the X-ray albedo of the disk to 0.95 (Chakrabarty 1998) and the mass transfer rate to 10⁻¹⁰ M_⊙ yr⁻¹ (Verbunt 1993; see Baglio et al. 2023 for further details). We performed a Markov Chain Monte Carlo (MCMC) sampling to investigate the posterior probability distribution of the model’s parameter space. The best-fit values for the free parameters are provided in Table 2, where they are also compared with the results from the modeling of the SED of J1023 in the high mode, as reported by Baglio et al. (2023). Each parameter is determined as the median of its marginal posterior distribution, with the 1σ credible intervals derived from the 16th to 84th percentiles of the posterior samples. The selected prior distributions are provided in the last column of Table 2. Figure 10 displays the corresponding histograms of the sampled parameters.

Fig. 10. Corner plot displaying the posterior probability distributions of the parameters obtained from the MCMC sampling algorithm. The diagonal panels show the probability density of each parameter, with solid lines representing the distribution and vertical dashed lines marking the 16th, 50th, and 84th percentiles. The median values and corresponding 1σ uncertainties are indicated above each panel. The off-diagonal panels illustrate the 2D posterior distributions, with contour lines representing the 1σ, 2σ, and 3σ equivalent bounds.

4. Discussion

We have presented the most extensive high-time resolution multi-wavelength campaign on the candidate tMSP J1544. The observed bimodality in X-ray intensity between high and low modes, along with an X-ray spectrum well described by an absorbed power law with a photon index of ∼1.6 – characteristic of tMSPs in the intermediate, sub-luminous disk state (see Papitto & de Martino 2022, and references therein) – further supports the classification of this source as a promising candidate. We also aimed to investigate a bimodal variability pattern in the UV data through simultaneous XMM-Newton and HST observations. However, due to the faint nature of the source and the variability in the background, we could not draw definitive conclusions regarding the presence of high and low modes in the UV band, which were observed for the brighter prototype of tMSPs, J1023 (Jaodand et al. 2021b; Miraval Zanon et al. 2022; Baglio et al. 2023).

The GTC/HiPERCAM light curves show flickering and dipping activities that mirror the X-ray variability. Moreover, they provide clear evidence of potential low modes, during which the optical emission is significantly redder than in high modes. This finding aligns with the mini-pulsar nebula scenario proposed by Papitto et al. (2019), which explains optical, UV, and X-ray pulsations from J1023 as due to synchrotron emission in the boundary region where the striped pulsar wind interacts with the inner accretion disk (see also Veledina et al. 2019). During the switch from high to low mode, characterized by the simultaneous disappearance of pulsations across all three bands, the inner flow is ejected, resulting in a decline in the optical, UV, and X-ray fluxes. Subsequently, as the inner flow begins to replenish, the system switches back to the high mode (Papitto et al. 2019; Baglio et al. 2023; see also Bogdanov et al. 2018 for a similar interpretation). The reddening of the optical emission at low fluxes, observed both in this work and for the candidate tMSP CXOU J110926.4−650224 in the sub-luminous disk state (Coti Zelati et al. 2024), is interpreted as a residual, fainter optical emission originating from the cooler, outer regions of the disk during the low modes. This result – along with other evidence such as the similar pulse shapes observed in the optical and X-ray bands, the simultaneous detection of pulsations occurring only during the high modes, and the pulsed SED described by a single power-law relationship (Papitto et al. 2019) – indicates that during the high modes, the majority of the optical emission originates from the boundary region between the pulsar wind and the accretion disk, just like the X-rays.

While analyzing archival XMM-Newton/OM data acquired in 2018, we identified hints of an optical flare around 58147.4 MJD with no corresponding feature in the X-ray light curve, marking the first such detection from this source. Specifically, after an observation gap around 58147.1 MJD, we observe a rise in optical flux followed by the flare-like feature. Similar optical spikes have been observed from J1023 with XMM-Newton/OM (although they always occur in correspondence with an X-ray flare; see, e.g., the bottom panel of Fig. 7(c) from Jaodand et al. 2016) as well as with other instruments (e.g., Shahbaz et al. 2015). A study of an 80-day-long uninterrupted Kepler monitoring campaign of J1023 in 2017 by Papitto et al. (2018) revealed that optical flares occurred for ∼15.6% of the time, a noticeably higher fraction than previously reported based on X-ray and optical observations. Kennedy et al. (2018) later reported an even higher occurrence rate of ∼22%, likely due to differences in the definition of a flare. Both studies concluded that the orbital dependence of these flares remains uncertain and is strongly influenced by the way flares are classified. However, Papitto et al. (2018) found that optical flares were more frequently detected when the companion star was at superior conjunction in its orbit. If confirmed, this would suggest that at least some of the flares originate from the reprocessing of X-ray emission off the companion star’s surface. Given that the 2018 OM observations were conducted in white light (spanning from the UV to the redder V-band in the optical range), we speculate that the rise in the optical flux and the hint of a flare may originate from the outer regions of the accretion disk or the companion star. Future long-term optical monitoring with optical telescopes, ideally with as much simultaneous X-ray coverage as possible, will be necessary to identify additional optical flares from J1544, determine whether they have an X-ray counterpart, and ultimately shed light on the physical mechanisms driving these events.

To further probe the multi-band variability of J1544, we analyzed two NIR observations acquired in 2016 with TNG/NICS. The resulting light curves displayed significant variability, possibly mimicking the X-ray bimodality between high and low modes. However, the estimated values of excess variance indicate that we cannot definitively confirm intrinsic NIR variability in this source. Investigating NIR variability in tMSPs seems to be a challenging task. The NIR light curves of the confirmed tMSPs J1023 and XSS J12270−4859 typically exhibit substantial variability across all filters in the sub-luminous disk state, including several flares and just a few, mildly pronounced dips that may correspond to the X-ray mode switching (de Martino et al. 2010, 2014; Saitou et al. 2011; Hakala & Kajava 2018; Shahbaz et al. 2018; Papitto et al. 2019; Baglio et al. 2019). In contrast, recent NIR observation of the candidate CXOU J110926.4−650224 in the sub-luminous disk state revealed a prolonged 20-minute dip (Coti Zelati et al. 2024). Future NIR observations simultaneous with high-time resolution X-ray observations are needed to better characterize the NIR variability for various confirmed and candidate tMSPs and to explore potential connections with X-ray high and low modes. Such studies could indeed help constrain the origin of the NIR emission, and thus the expected variability pattern, by determining whether it arises from a jet (e.g., Baglio et al. 2019), the outer accretion disk, the companion star, or some combination of these components (see, e.g., Papitto et al. 2019 for tentative evidence of correlated NIR and X-ray variability for J1023; see also Baglio et al. 2023 for subsequent lack of flares at the low-to-high mode switch and of any clear signs of modeswitching).

We performed radio observations with the VLA simultaneously with XMM-Newton and later, in November and December 2024, with ATCA. In the following, we focus on the VLA observations, which provided more stringent upper limits because the ATCA observations were affected by higher-than-expected RMS noise due to the presence of significant RFI in our spectral windows. Even when VLA radio images were extracted separately during the X-ray high and low modes – an analysis motivated by the anti-correlated variability observed between the radio and X-ray emissions in J1023 (Bogdanov et al. 2018; Baglio et al. 2023) and a similar trend suggested for J1544 (Gusinskaia et al. 2025) – we did not detect the radio counterpart to J1544 with VLA. Our non-detection supports the radio variability previously observed (Jaodand et al. 2021b; Gusinskaia et al. 2025), despite the relatively stable X-ray flux observed over time (see Fig. 8). To quantify the X-ray flux stability, we used the estimates from Gusinskaia et al. (2025), as they are, to our knowledge, the only ones in the literature derived exclusively during the high modes. Following the procedure outlined in Sect. 3.1, we extracted the 0.5–8 keV unabsorbed flux during the high mode and found F_unabs, H = (4.67 ± 0.03)×10⁻¹² erg cm⁻² s⁻¹, showing stability within a factor of 1.5 compared to the Chandra 2018 and 2019 observations reported byGusinskaia et al. (2025).

Jaodand et al. (2021b) observed J1544 with the VLA in the X-band (8–12 GHz) during four epochs in 2015. These observations were quasi-simultaneous with Swift/XRT observations, which revealed unabsorbed 1–10 keV X-ray flux values ranging from (0.55 ± 0.34)×10⁻¹² erg cm⁻² s⁻¹ to (4.21 ± 0.43)×10⁻¹² erg cm⁻² s⁻¹. During this period, the VLA radio flux density varied significantly, from (47.7 ± 6.0) μJy to < 13.8 μJy (3σ upper limit). This variability was further highlighted in subsequent VLA observations reported by Gusinskaia et al. (2025). A 2018 X-band (8–12 GHz) observation yielded a flux density of (11.9 ± 1.6) μJy, increasing up to (14.2 ± 2.9) μJy during simultaneous X-ray low-mode intervals. Similarly, a 2019 C-band (4–8 GHz) observation measured (28.6 ± 2.1) μJy, increasing up to (56.6 ± 8.3) μJy during X-ray low-mode intervals. In both cases, the unabsorbed 1–10 keV flux detected by Chandra remained consistent with the range reported by Jaodand et al. (2021b). Interestingly, during our VLA observations in February 2024, conducted in the C-band, J1544 was not detected, even though the 1–10 keV X-ray flux during our XMM-Newton observation was (2.80 ± 0.04)×10⁻¹² erg cm⁻² s⁻¹, consistent with the 2015 Swift/XRT values. In the bottom panel of Fig. 8, we compare our stringent upper limit on the radio flux density with the 2019 measurement in the same band by Gusinskaia et al. (2025). This comparison shows that if J1544 had maintained a similar flux density in our recent observation, we would have detected it. Our non-detection thus underscores the pronounced radio variability of the source. Radio variability without significant changes in X-ray flux or spectral state has been observed in other accreting NSs at higher accretion rates, underscoring the still-mysterious connection between accretion flow and the launching of outflows in these systems (see, e.g., Marino et al. 2023; Panurach et al. 2023; Pattie et al. 2024; also see the recent evidence of compact jet formation in tMSPs by Koljonen et al. 2025). Future, longer radio observations of J1544 – ideally simultaneous with XMM-Newton – will be highly valuable not only for understanding its radio variability and potential correlation with X-ray high/low modes but also for providing new insights into the broader, poorly understood mechanisms that drive mass accretion and ejection in NS X-raybinaries.

A significant step toward confirming J1544 as a tMSP is the striking similarity of its broadband (from UV to X-rays) SED to that of J1023 as reported by Miraval Zanon et al. (2022) (see left panel of Fig. 9). The SED modeling in the high mode (see right panel of Fig. 9), along with the reddening of optical emission at lower fluxes observed with GTC/HiPERCAM and discussed above, supports the mini-pulsar nebula scenario. Table 2 presents the best-fit parameter values for J1544, compared with those of J1023 from Baglio et al. (2023), further reinforcing the resemblance between the two sources. The main difference is that, while for J1023 the irradiation luminosity was fixed at the value of L_irr = 6.5 × 10³³ erg/s estimated by Shahbaz et al. (2019), the limited knowledge of the system’s parameters for J1544 required treating this parameter as a free variable. In the following, we discuss the obtained value of L_irr.

It is important to note that Britt et al. (2017) found no significant evidence of phase-dependent temperature variations, suggesting minimal irradiation of the companion star. Therefore, we adopted the effective temperature of Britt et al. (2017) as the prior bound on T_* (see Table 2) to account for the emission of the non-irradiated companion star. They also noted that the absence of irradiation signatures in their light curves could be attributed to the system’s nearly face-on inclination. However, in our analysis, we find that the contribution of irradiation is necessary to model the broadband flux distribution, both in the star and the accretion disk components. Keeping these considerations in mind, we now examine the relationship between the observed irradiation luminosity and the spin-down power. Assuming isotropy, the energy flux from the pulsar wind at the location of the companion star is given by F_irr = Ė/(4πa²), where Ė is the spin-down power and a is the orbital separation. The fraction of this energy intercepted by the companion star depends on its projected area, leading to f = R_c²/(4a²), where R_c is the companion star’s radius. Accounting for an efficiency factor η that quantifies the reprocessing of the absorbed energy into observable radiation, the irradiation luminosity can be expressed as $L_{\mathrm{irr}}=\eta [R_{\mathrm{c}}^2/(4a^2)]\dot{E}$. Considering that η ≤ 1 and the value reported in Sect. 3.7, the previous relation implies that the spin-down power for J1544 should be Ė ≤ (6 ± 1) × 10³⁴ erg/s. While this approach may be valuable – since it could, in principle, constrain the irradiation luminosity and, by extension, the spin-down power for confirmed and candidate tMSPs in the sub-luminous disk state that have not been observed in the rotation-powered state – it should be interpreted with caution. Notably, the upper limit obtained for J1544 is an order of magnitude lower than the estimate derived using the empirical correlation between X-ray luminosity and spin-down power recently proposed by Xu et al. (2025), combined with our estimation of the X-ray luminosity. In the future, should J1544 transition to the rotation-powered state, deep radio observations with single-dish telescopes could enable the detection of its rotational parameters and unambiguously provide a measure of its spin-down power.

5. Conclusions

With this paper, we have presented the most extensive high-time resolution multi-wavelength campaign conducted on the candidate tMSP J1544 in the sub-luminous disk state. Our key findings are summarized as follows:

• X-ray/UV emission: The X-ray data clearly exhibit a bimodal pattern between high and low modes. This is consistent with previous studies and strongly supports the classification of J1544 as a promising tMSP candidate in the sub-luminous disk state. The UV emission showed marked variability. In some cases, these variations appear to correspond to changes in the X-ray intensity level. However, due to the source faintness and background variability, we were unable to definitively establish bimodal behavior in the UV data.

• Optical/NIR variability: High-time resolution optical observations showed variability patterns similar to the X-ray modes, with notable reddening of optical emission at low fluxes, supporting the mini-pulsar nebula scenario. This suggests that during low modes, residual optical emission comes from the cooler outer regions of the disk after the inner flow is ejected, while in high modes, most optical emission originates from the boundary region between the pulsar wind and the accretion disk, as is the case for the X-rays. The NIR variability was too noisy to draw firm conclusions, emphasizing the need for further multi-wavelength monitoring.

• Optical flare candidate: We report the first candidate optical flare in archival XMM-Newton/OM data, which lacks an X-ray counterpart. This suggests the presence of complex variability mechanisms and motivates future coordinated optical and X-ray monitoring over extended periods.

• Radio flux variability: Although the X-ray flux of J1544 has remained relatively stable in recent years, the source has exhibited significant radio variability, alternating between detections and non-detections under similar observational conditions. This underscores the complexity of the mass ejection processes in tMSPs.

• SED modeling and spin-down power estimate: The broadband SED of J1544 is strikingly similar to that of the tMSP J1023 in the sub-luminous disk state, supporting the mini-pulsar nebula scenario. By treating the irradiation luminosity as a free parameter, we showed that irradiation is essential for modeling the broadband flux from both the companion star and the accretion disk. Using the relationship between irradiation luminosity and spin-down power, we estimated J1544’s spin-down power to be an order of magnitude lower than predicted by the empirical correlation between X-ray luminosity and spin-down power for rotation-poweredpulsars.

Acknowledgments

The research leading to these results has received funding from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158). This work is based on observations made with the Gran Telescopio Canarias (GTC), installed at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, on the island of La Palma (program ID: GTC129-23B). This work is based on data obtained with the instrument HiPERCAM, built by the Universities of Sheffield, Warwick, and Durham, the UK Astronomy Technology Centre, and the Instituto de Astrofísica de Canarias. Development of HiPERCAM was funded by the European Research Council, and its operations and enhancements by the Science and Technology Facilities Council. SiFAP2 and NICS observations were made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundación Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work is also based on observations acquired with the NICER mission, a 0.2–12 keV X-ray telescope operating on the International Space Station; the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by NASA; and XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Space Science Data Center (SSDC, Italy) and the California Institute of Technology (Caltech, USA). We also used software and tools provided by the High Energy Astrophysics Science Archive Research Center (HEASARC) Online Service. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Australia Telescope Compact Array is part of the Australia Telescope National Facility (https://ror.org/05qajvd42), which is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Gomeroi people as the Traditional Owners of the ATCA observatory site. G.I. thanks the HST program coordinator, A. Vick (STScI), for constant support in the observation planning and A. Fullerton (STScI) for checking the scheduling processes. G.I. thanks S. Dieterich (STIS Team) for the support in the scientific data analysis. G.I. and A.M.Z. also thank M. Cadelano for the useful discussion on HST data analysis. G.I., A.P., F.C.Z., S.C., D.d.M., C.M., R.L.P. are supported by INAF (Research Grant ‘Uncovering the optical beat of the fastest magnetised neutron stars 620 (FANS)’) and the Italian Ministry of University and Research (MUR) (PRIN 2020, Grant 2020BRP57Z, ‘Gravitational and Electromagnetic-wave Sources in the Universe with current and next-generation detectors (GEMS)’). A.P. acknowledges support from the Fondazione Cariplo/Cassa Depositi e Prestiti, grant no. 2023-2560. F.C.Z. acknowledges support from a Ramon y Cajal fellowship (grant agreement RYC2021-030888-I). N.R. and A.M. are supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (ERC Consolidator Grant “MAGNESIA” No. 817661). F.C.Z. and N.R. acknowledge support from grant SGR2021-01269 from the Catalan Government. This work was also supported by the Spanish program Unidad de Excelencia Marıa de Maeztu CEX2020-001058-M and by MCIU with funding from European Union NextGeneration EU (PRTR-C17.I1). T.D.R. is an INAF research fellow. D.F.T. has been supported by PID2021-124581OB-I00 funded by MCIU/AEI/10.13039/501100011033 and 2021SGR00426 as well as by funding from the European Union NextGeneration EU (PRTR-C17.I1) program.

References

Appendix A: Spectral analysis of archival XMM-Newton data

To enhance the statistics of the spectra separately extracted in high and low modes and to investigate significant changes in the photon index between these states, we analyzed XMM-Newton archival observations of J1544 performed on 2014 February 16 (ObsID 0724080101; PI: Turriziani) and on 2018 January 28 (ObsID 0800280101; PI: Bogdanov). ObsID 0724080101 is detailed in Bogdanov & Halpern (2015), whereas during ObsID 0800280101 the EPIC-pn operated in timing mode, and the two EPIC-MOS in small window mode. We followed the data analysis procedure described in Sect. 2.1, carefully filtering out the high background flaring activity observed in the 10–12 keV light curves for ObsID 0800280101 (see middle panel of Fig. B.1).

The high-mode spectrum (left panel of Fig. A.1) includes XMM-Newton and NuSTAR data as described in Sects. 2.1 and 2.6, together with EPIC-pn and EPIC-MOS1 data from ObsID 0724080101, and EPIC-pn, EPIC-MOS1, and EPIC-MOS2 data from ObsID 0800280101. The low-mode spectrum (right panel of Fig. A.1) comprises EPIC-MOS1 data from ObsID 0724080101, and EPIC-MOS1 and MOS2 from ObsID 0800280101, excluding the background-dominated EPIC-pn data and low-statistic data from NuSTAR. The resulting total exposure is ∼467.97 ks during the high-mode intervals and ∼112.38 ks during the low-mode intervals.

We fitted the spectra with an absorbed power law, keeping the absorption column density fixed to the value obtained from the average spectral fitting (Sect. 3.1). We included a renormalization factor in the model to account for cross-calibration uncertainties between the two X-ray telescopes and different observations. For the high mode, these values remained consistent within 10%, while for the low mode, they showed a larger discrepancy up to 14%, likely due to the limited statistics in this mode. For the high mode, we obtained Γ_H = 1.627 ± 0.006 and an unabsorbed 0.3-10 keV flux F_unabs, H = (9.4 ± 0.1)×10⁻¹² erg cm⁻² s⁻¹ (χ²/d.o.f. = 1100.12/1100), while for the low mode Γ_L = 1.66 ± 0.06 and F_unabs, L = (0.38 ± 0.03)×10⁻¹² erg cm⁻² s⁻¹ (χ²/d.o.f = 139.0/136). No significant change in the photon index between high and low modes was reported by Bogdanov & Halpern (2015), and our results likewise show consistency within 1σ across modes.

Fig. A.1. X-ray spectra extracted separately in the high-mode (top panel) and low-mode (bottom panel) intervals, along with the best-fitting models (see Appendix A for details). Top panel: Various shades of blue denote the three EPIC-pn spectra, while different shades of red represent the spectra obtained with the two EPIC-MOS cameras. Light and dark green indicate NuSTAR spectra. Bottom panel: As discussed in the text, the low-mode spectrum includes only data for the two EPIC-MOS cameras.

Appendix B: Inspection of the optical flare candidate

To investigate the nature of the flare-like event detected in the archival XMM-Newton/OM data from 2018 January 28 (ObsID 0800280101; see Sect. 3.4), we first confirmed that the observed feature was absent in the OM background light curve. During the OM processing using omfchain, the background level is derived from imaging-mode data using the recommended setting bkgfromimage=yes, rather than from the fast-mode window. This approach uses the average background obtained during the whole observing window as image mode and thus is assumed to be constant in the standard processing. For this reason, we instead used the background light curve as extracted in a region within the fast window (see top panel of Fig. B.1). Since the flare-like feature is only present in the data and not in the background light curve, this supports its genuine nature. We then carefully examined the XMM-Newton/EPIC data for any intervals of flaring particle background. We used the evselect command with the following selection expressions: “#XMMEA_EM && (PI> 10000) && (PATTERN= = 0)” for EPIC-MOS, and “#XMMEA_EP && (PI> 10000&&PI< 12000) && (PATTERN= = 0) for EPIC-pn”, using 100-s bins. As shown in the middle panel of Fig. B.1, the optical feature does not coincide with periods of enhanced particle background in the X-rays. Finally, we visually inspected all individual OM images and verified that the source consistently remained within the Field of View throughout the observation.

Fig. B.1. Temporal evolution of the optical emission as observed with XMM-Newton/OM in 2018. Top panel: OM light curve (black) and background light curve (purple), both binned at 100 s. Middle panel: Comparison between the XMM-Newton/OM light curve (100-s bins, shown in black) and the EPIC high-energy light curves in the 10-12 keV band (100-s bins; light blue for EPIC-pn, blue for EPIC-MOS1, and dark blue for EPIC-MOS2) extracted from the event files to identify intervals of flaring particle background. The y-axis is shown on a logarithmic scale for visual purposes. Bottom panel: Overlap of the XMM-Newton/OM light curve (100-s bins, shown in black) and the XMM-Newton/EPIC background-subtracted light curve (20-s bins, shown in green).

All Tables

All Figures

Fig. 1. Distribution of count rates obtained from the background-subtracted XMM-Newton/EPIC light curve acquired on 2024 February 8 binned with a time resolution of 20 s. We defined the high modes as the time intervals when the count rate exceeded 0.8 counts s−1 (red dashed line), and the low modes as those when the count rate dropped below 0.4 counts s−1 (blue dashed line). A “transition” region is also present, following the procedure adopted by Bogdanov et al. (2015) for J1023. No flaring mode was observed.

In the text

Fig. 2. Temporal evolution of the X-ray, UV, and optical emissions of J1544 during the first three days of observations. For Day 1, the light curves are shown in decreasing order of energy band, from top to bottom, with the XMM-Newton/EPIC (20 s time bin), HST/STIS (20 s), and XMM-Newton/OM (1200 s). For Days 2–3, the 3–30 keV NuSTAR light curve (100 s time bin) is shown. The yellow-shaded areas denote the time intervals of the low X-ray mode identified by XMM-Newton/EPIC on Day 1 and by NuSTAR on Day 2–3. The error bars represent 1σ uncertainties.

In the text

	Fig. 3. Zoom-in of the time series collected on Day 1 capturing the period of maximum overlap among different telescopes.
In the text

	Fig. 4. GTC/HiPERCAM light curves acquired on Day 4 in the Super SDSS uS, gS, rS, iS, and zS filters. The yellow-shaded areas highlight the most prominent potential low modes. The error bars represent 1σ uncertainties.
In the text

	Fig. 5. Color-magnitude diagrams for J1544 derived from the most highly sampled gS, rS, iS, and zS bands of the GTC/HiPERCAM time series (Fig. 4) and rebinned to emphasize correlation trends. The upper panels display the smoothed Kernel Density Estimation curves that illustrate the distribution of the color-magnitude values.
In the text

	Fig. 6. XMM-Newton/OM light curve acquired in 2018 using 100-s bins. The error bars represent 1σ uncertainties.
In the text

	Fig. 7. Near-infrared light curves of J1544 observed with TNG/NICS on 2016 March 29 (top panel) and 30 (bottom panel). We plotted differential magnitudes by subtracting the weighted average value in each filter for each night (see Sect. 3.5 for details). Photometric data were acquired in the J (blue points), H (green points), and K (red points) filters. The error bars represent 1σ uncertainties.
In the text

Fig. 8. Radio variability of J1544 compared to its relatively stable X-ray flux over the years. This figure extends Fig. 1 from Jaodand et al. (2021b) by incorporating observations from Gusinskaia et al. (2025) and this work. Top panel: Unabsorbed X-ray flux in the 1–10 keV range (estimated using the entire observations, including both high and low modes) from Jaodand et al. (2021b) (dark blue squares), Gusinskaia et al. (2025) (light blue hexagons), and this work (green point). For consistency, we extracted the flux in the 1–10 keV band from our XMM-Newton observation, obtaining FX = (2.79 ± 0.04)×10−12 erg cm−2 s−1, by following the procedure detailed in Sect. 3.1. The low X-ray flux around MJD 57179 is likely due to the source remaining in the X-ray low mode throughout this short (418 s) observation, as noted by Jaodand et al. (2021b). Bottom panel: VLA radio flux density in the X band (8–12 GHz) from Jaodand et al. (2021b) (dark red squares are for detections and downward triangles are for 3σ upper limits) and from Gusinskaia et al. (2025) at 8–12 GHz (X band) around 58224 MJD (red hexagon) and at 4–8 GHz (C band) around 58491 MJD (orange hexagon). The 3σ upper limit of ∼9 μJy at 4–8 GHz (C band) from this work is indicated by the brown downward triangle. For consistency with Jaodand et al. (2021b), the plot shows the radio flux density from Gusinskaia et al. (2025) and the upper limit derived in this work both extracted using the entire VLA observation, rather than restricting to intervals simultaneous with X-ray low modes. The dashed vertical lines mark the epochs of X-ray observations for reference.

In the text

Fig. 9. Unabsorbed broadband SED of J1544. Left panel: SED of J1544 from UV to X-rays compared with that of J1023 from Miraval Zanon et al. (2022). For J1023 and J1544, HST data are shown in orange and red, XMM-Newton data are in light blue and dark blue, and NuSTAR data are in light green and dark green, respectively. The HST spectrum is plotted from 165 to 310 nm (i.e., ∼(1.0 − 1.8)×1015 Hz) for J1023 and from 200 to 300 nm (i.e., ∼(1.0 − 1.5)×1015 Hz) for J1544. The UV data were rebinned using the coronagraph Python package (Robinson et al. 2016; Lustig-Yaeger et al. 2019) with a low-resolution wavelength grid of width 4 for J1023 and width 20 for J1544. Right panel: The SED of J1544 from the optical band to X-rays, extracted during X-ray high modes observed with XMM-Newton, along with the best-fitting model (see text for details). Due to the lack of simultaneous X-ray observations, we include all available optical data from GTC/HiPERCAM. The red dashed, blue dotted, and orange solid lines represent the contributions from the accretion disk, companion star, and boundary region, respectively, while the green solid line shows their combined emission. The ratio between the data points and the best-fitting model is shown in the bottom panel.

In the text

Fig. 10. Corner plot displaying the posterior probability distributions of the parameters obtained from the MCMC sampling algorithm. The diagonal panels show the probability density of each parameter, with solid lines representing the distribution and vertical dashed lines marking the 16th, 50th, and 84th percentiles. The median values and corresponding 1σ uncertainties are indicated above each panel. The off-diagonal panels illustrate the 2D posterior distributions, with contour lines representing the 1σ, 2σ, and 3σ equivalent bounds.

In the text

Fig. A.1. X-ray spectra extracted separately in the high-mode (top panel) and low-mode (bottom panel) intervals, along with the best-fitting models (see Appendix A for details). Top panel: Various shades of blue denote the three EPIC-pn spectra, while different shades of red represent the spectra obtained with the two EPIC-MOS cameras. Light and dark green indicate NuSTAR spectra. Bottom panel: As discussed in the text, the low-mode spectrum includes only data for the two EPIC-MOS cameras.

In the text

Fig. B.1. Temporal evolution of the optical emission as observed with XMM-Newton/OM in 2018. Top panel: OM light curve (black) and background light curve (purple), both binned at 100 s. Middle panel: Comparison between the XMM-Newton/OM light curve (100-s bins, shown in black) and the EPIC high-energy light curves in the 10-12 keV band (100-s bins; light blue for EPIC-pn, blue for EPIC-MOS1, and dark blue for EPIC-MOS2) extracted from the event files to identify intervals of flaring particle background. The y-axis is shown on a logarithmic scale for visual purposes. Bottom panel: Overlap of the XMM-Newton/OM light curve (100-s bins, shown in black) and the XMM-Newton/EPIC background-subtracted light curve (20-s bins, shown in green).

In the text