© ESO 2021

1. Introduction

Binary millisecond pulsars (MSPs) are believed to be the end product of the accretion process in low-mass X-ray binary systems (LMXBs; Alpar et al. 1982; Radhakrishnan & Srinivasan 1982). MSPs attain short spin periods (P_s ≲ 10 ms) by accreting matter and angular momentum from a companion star in a process called “recycling”, which can last up to a few gigayears (Tauris et al. 2013). During this prolonged accretion phase, LMXB systems are observed as bright X-ray sources (Campana et al. 1998), mainly due to the presence of an accretion disk around the neutron star (NS). The optical spectra are characterized by strong, double-peaked emission lines of H and He produced in the optically thick accretion disk (van Paradijs & McClintock 1995). When the mass transfer rate decreases such a X-ray bright phase ends, double-peaked emission lines disappear, and emission in the radio and γ-ray bands, powered by the NS’s rotation, becomes predominant. The recycling model has been strongly supported by observations of X-ray millisecond pulsations in transient accreting millisecond X-ray pulsars (AMXPs; Wijnands & van der Klis 1998) and coherent X-ray pulsations during Type I X-ray bursts (Strohmayer et al. 1996). The ultimate proof of this evolutionary scenario came with the discovery of transitional millisecond pulsars (tMSPs, see Papitto & de Martino 2020 for a review), which are systems that have shown transitions between LMXB and radio MSP states on very short timescales (Archibald et al. 2009; Papitto et al. 2013; Bassa et al. 2014).

“Spider” pulsars are a subclass of binary MSPs in which the pulsar is in a very tight orbit (≲1 day) around a low-mass companion star that is losing mass (Roberts 2013). The ionized gas lost by the companion can intercept and thus distort or even totally absorb the radiation coming from the NS, causing regular eclipses of the radio pulsar signal. After the discovery of the first eclipsing binary pulsar, PSR B1957+20 (Fruchter et al. 1988a,b), many similar systems have been found both in globular clusters¹ and in the Galactic field². Spider pulsars are divided into two further subpopulations (Roberts 2013): black widows (BWs) and redbacks (RBs). BW pulsars are characterized by a very low-mass (⩽0.1 M_⊙), semi-degenerate companion star. The companion mass loss is typically driven by the evaporation of the outer layer of the star, caused by the irradiation of the surface by the strong pulsar wind. RB systems have more massive nondegenerate companions (0.1–0.4 M_⊙) that are only partially ablated. Most systems of both classes show regular eclipses of their radio signal, which cover a small fraction (< 20%) of the orbit in the BWs and a more significant fraction (up to 80%) of the orbit in the RBs (Broderick et al. 2016; Polzin et al. 2020).

These spiders are key systems to study the evolution of binary MSPs and even more to investigate the formation of isolated MSPs. All the tMSPs discovered so far, which represent an important breakthrough in the knowledge of binary MSPs, are RB systems. Isolated MSPs might be originated from BW systems when their original companion star disappears, getting completely ablated by energetic particles produced by the pulsar wind or by γ-ray particles powered by the pulsar spin-down energy (Levinson & Eichler 1991; Tavani & London 1993). However, simulations indicate that the direct evaporation alone cannot be responsible for the formation of most isolated MSPs (Ginzburg & Quataert 2020). Ginzburg & Quataert (2021) propose that the ablation wind when coupled to the companion magnetic field is able to remove orbital angular momentum from the binary through magnetic braking, maintaining a stable Roche-lobe overflow and thus leading to the total evaporation of the companion star. These systems are also excellent laboratories to study a matter-radiation interaction in the low-density and high magnetic field regime. These conditions are typical of the shocks that originated in the interaction between matter released by the companion star and the relativistic pulsar wind. Moreover, RBs and BWs are crucial to constrain the maximum mass of a NS (see Linares 2020 and Strader et al. 2019 for a review).

In this paper we investigate both these aspects, constraining the shock geometry and the NS mass in a RB system. The Doppler tomography method (Marsh & Horne 1988) allows us to derive two-dimensional velocity maps of the Hα emission region, by combining spectra taken at different orbital phases along the whole orbit. This technique enables us to distinguish between an accretion disk pattern and the more complicated geometry and kinematics that are expected in the pulsar shock emission scenario. This method has not yet been exploited for RB systems, because only a few of them clearly show emission lines in their optical spectra during the radio pulsar state (for example PSR J1628–3205, Cho et al. 2018).

PSR J1048+2339 (hereafter J1048) was discovered with the Arecibo telescope in a search for high-latitude Fermi unidentified sources (Cromartie et al. 2016). Deneva et al. (2016) classified it as a RB system using multiwavelength observations. J1048 is a 4.66-ms-spinning NS with a spin-down power of 1.2 × 10³⁴ erg s⁻¹. It is part of a tight binary system with an orbital period of 6.01 hours and it is located at a distance of about 1.7 kpc (Deneva et al. 2021). Its companion star has a mass of ∼0.4 M_⊙ and an effective temperature T_eff ∼ 4200 K (Yap et al. 2019). The inclination of the system is not well-constrained. Yap et al. (2019) fix the inclination to a theoretical upper limit of 76° because X-ray eclipses in Chandra data are not observed. Contrastingly, spectroscopic studies from Strader et al. (2019) show that either the system must be close to edge-on ($i\ge {83}_{-10}^{+7}$ degrees), or the pulsar mass must be higher than 2 M_⊙. The system exhibits radio eclipses for a large fraction of the orbit (∼50%, Deneva et al. 2016, 2021). The radio pulsar has a γ-ray counterpart, 4FGL J1048.6+2340, which shows γ-ray pulsations at a statistical significance of ∼8σ (Deneva et al. 2021). Optical and X-ray studies (Cho et al. 2018; Yap et al. 2019) reveal a strong orbital modulation with rapid variations: the optical light curve changes from an ellipsoidal to a sinusoidal-like profile in less than 14 days. The average optical magnitudes when the pulsar wind heating dominates are r′∼19.4 mag and g′∼20.7 mag (Yap et al. 2019). Instead, when the orbital modulation is dominated by ellipsoidal modulation, the system is fainter by 0.3–0.4 mag. The X-ray light curve is also strongly modulated with a minimum centered at ϕ = 0.0, which corresponds to the inferior conjunction of the companion star in our convention. Optical spectra acquired with the 4-m SOAR telescope show the most intense variability yet observed for any RBs (Strader et al. 2019). The strong emission lines in some instances appear double-peaked, sometimes single-peaked and are absent at other orbital phases. The multiepoch poor orbital phase sampling precludes a secure origin of the emission lines region.

In this paper we present a multiwavelength campaign of J1048. Using Very Large Telescope (VLT), Sardinia Radio Telescope (SRT), Low-Frequency Array (LOFAR) telescope, Galileo telescope, and Swift simultaneous observations carried out during the radio pulsar state, we study the kinematics and geometry of matter in this interacting binary system. In Sect. 2 we describe our data set and data reduction steps. In Sect. 3 we present the average spectrum, the variability of Hα emission line and the optical light curve, while the radial velocity curves are presented in Sect. 4. In Sect. 5 we show the Hα Doppler tomography. Sections 6–7 are dedicated to radio and optical timing analysis. In Sect. 8, we discuss the results and we present future perspectives.

2. Observations

In March 2020, we performed a multiwavelength campaign focused on J1048 with simultaneous observations in the X-ray, UV, optical and radio bands. Table 1 lists the observations analyzed and discussed in this paper. In the following, we describe the data analysis for the different data sets acquired.

2.1. Optical and ultraviolet observations

2.1.1. Very Large Telescope

Optical spectroscopic observations of J1048 were carried out on March 19, 2020 at the Observatorio Monte Paranal, using the ESO 8.2 m Antu Telescope (UT1) equipped with the FORS2 Spectrograph (Appenzeller et al. 1998). During the night, the seeing ranged between 0.6 − 0.9″, with most of the time stable between 0.7 − 0.8″. A total of sixteen spectra were obtained, covering a complete orbital cycle, using 1200 s long exposures. We used a 1″ slit, always oriented along the parallactic angle, and the grism 1200R+93, centered at 6500 Åand with a dispersion of 0.76 Å/pixel. In order to extract the optical light curve, we acquired one or two R-band images with 20 s exposure time, once every two spectroscopic acquisitions, for a total of 10 images.

Image reduction was carried out following standard procedures: subtraction of an averaged bias frame, division by a normalized flat frame. Astrometry was performed using the USNOB1.0 catalog³. Aperture photometry was made with the PHOTOM software part of the STARLINK⁴ package. The photometric calibration was done against Stetson standard stars (Stetson 2000). In order to minimize any systematic effect, we performed differential photometry with respect to a selection of local isolated and nonsaturated reference stars. The reduction and extraction of the spectra was performed with the ESO-FORS pipeline⁵ and the ESO-MIDAS⁶ software package. Wavelength and flux calibration of the spectra were achieved using helium-argon lamp and observing spectrophotometric standard stars, respectively, at the beginning and end of the observing run. We further checked stability of the wavelength calibration using night-sky emission lines. For the cross-correlation of the spectra with late-type templates and the Doppler tomography analysis we employed the MOLLY package developed by Tom Marsh⁷ and the pydoppler⁸ package (Spruit 1998), respectively.

2.1.2. Galileo/IFI+Iqueye telescope

We observed J1048 with IFI+Iqueye mounted at the 1.2 m Galileo telescope in Asiago, Italy. Iqueye+⁹ is a fast photon counter with a field of view of 12.5 arcsec and the capability of time tagging the detected photons with sub-ns time accuracy (Naletto et al. 2009). The instrument is fiber-coupled with the Galileo telescope by means of a dedicated instrument, the Iqueye Fiber Interface (IFI; Zampieri et al. 2019). A total of 3 acquisitions (obs. ID 20200319-002631, 20200319-020229, 20200319-032617) were performed on March 18–19, 2020, between 23:26:33.4 and 02:56:18.0 UTC. The first two acquisitions lasted 3600 s and the third 1800 s. The sky background was regularly monitored between on-target observations. The average rate at the position of the pulsar varied between ∼1300 and ∼1600 counts s⁻¹ because of variations in the sky background. The companion star is not significantly detected.

The data reduction is performed with a dedicated software¹⁰. The whole acquisition and reduction chain ensures an absolute accuracy of ∼0.5 ns relative to UTC (Naletto et al. 2009). The photon arrival times are barycentered using TEMPO2 in TDB time units (Hobbs et al. 2006) and the JPL ephemerides DE421. The position of the pulsar is RA = 10:48:43.4183, Dec = +23:39:53.411 at MJD 56897.0 (Deneva et al. 2016). We corrected for the motion of the pulsar along the orbit using the orbital period P_b = 0.25051904499812022 days and projected semi-major axis x_p = 0.836122 light-seconds of Deneva et al. (2016), while the time of passage at the ascending node T_asc = 58721.666409 MJD (Barycentric Dynamical Time) was derived from a SRT observation on August 19, 2019 (see Sect. 6).

2.1.3. Swift-UVOT

The Ultraviolet Optical Telescope (UVOT; Roming et al. 2005) on-board the Neil GehrelsSwift Observatory (Gehrels et al. 2004) observed J1048 with the U (λ_c = 346.5 nm), B (λ_c = 439.2 nm), V (λ_c = 546.8 nm), UVW1 (λ_c = 260.0 nm), UVM2 (λ_c = 224.6 nm), and UVW2 (λ_c = 192.8 nm) filters for ∼3.4 ks on March 18–19, 2020 (see Table 2). We performed aperture photometry by using a circular region centered on the source with a radius of 5″ for the UV filters and 3″ for the optical filters. The background emission was extracted from a nearby source-free region. The target was not detected in any of the images. Using the tool uvotsource of the Swift analysis software in the HEASoft package version 6.28¹¹, we estimated the following 3σ upper limits on the optical/UV magnitude: U > 19.9 mag, B > 20.3 mag, V > 19.3 mag, UVW1> 19.5 mag, UVM2> 19.6 mag, and UVW2> 20.1 mag (Vega system).

2.2. Radio observations

We observed J1048 on two different occasions in the radio band.

On August 19, 2019 we used the 64-m SRT to carry out a 1.5-hour-long observation with the dual-band coaxial L-P receiver (Valente et al. 2010). In the L band (∼1.3 − 1.8 GHz), the total-intensity signal was recorded every 125 μs with the Australia Telescope National Facility (ATNF) Digital Filterbank Mark III backend (DFB¹²) in search mode. The nominal bandwidth was chosen to be of 1024 MHz, centered at 1548 MHz and divided into 1024 frequency channels. This setup was chosen to prevent aliasing effects within the actual receiver observing band, which is ∼500 MHz wide. The P band was instead acquired in baseband mode with the Reconfigurable Open Architecture Computing Hardware¹³ (ROACH-1) backend, recording 80 MHz of bandwidth centered at a frequency of 336 MHz. Thanks to the presence of a Gregorian cover, which shields the receiver from radio signals generated by the telescope instrumentation, this observation was mostly free of Radio Frequency Interference (RFI) in both bands.

The second radio observation was carried out in the night between March 18–19, 2020. This time, J1048 was observed simultaneously with the SRT and the Low-Frequency Array telescope (LOFAR, Stappers et al. 2011; van Haarlem et al. 2013), for more than a full orbit. The SRT observation lasted ∼8.3 h, and was carried out using exactly the same setup as in the previous observation, with the only difference being that the Gregorian cover was not mounted. This resulted in a higher presence of RFI in these data. LOFAR started the simultaneous observation of J1048 in the VHF band (110 − 188 MHz) with the high-band antennae of the LOFAR core about 50 min later than SRT, and recorded data for 7 h. The data were collected as beamformed, 8-bit complex voltages from the two linear polarizations, with a sampling time of 5.12 μs and divided into 400 frequency sub-bands, with a bandwidth of 0.195 MHz each.

2.3. X-ray observations

We analyzed data sets acquired on March 18–19, 2020 using the Swift X-ray Telescope (XRT; Burrows et al. 2005) configured in the photon counting mode (timing resolution of 2.5 s; obs. IDs 00034285004 and 00034285005). The source was observed for 3.6 ks, summing up the two pointings. The source is barely detected, with a signal-to-noise ratio of S/N ∼ 2.2 as evaluated by the sosta tool in the HEASoft package. Its net count rate is (2.0 ± 0.9) × 10⁻³ counts s⁻¹ over the 0.3–10 keV energy band (uncertainty at a confidence level of 1σ). Assuming an absorbed power law model with absorption column density of N_H = 2.5 × 10²⁰ cm⁻² and photon index of Γ = 1.5 (Cho et al. 2018), the above count rate translates into an observed flux of (9 ± 4) × 10⁻¹⁴ erg cm⁻² s⁻¹ and an intrinsic unabsorbed flux of (9.5 ± 4.0) × 10⁻¹⁴ erg cm⁻² s⁻¹, or a luminosity of (3.3 ± 1.4) × 10³¹ erg s⁻¹ (0.3–10 keV) for a distance of 1.7 kpc (Deneva et al. 2021). These values are compatible within the uncertainties with those derived by Cho et al. (2018) and Yap et al. (2019) using Chandra deeper observations, and indicate that, at the epoch of our multiband campaign, the system remained in a rotation-powered radio pulsar state.

3. Average spectrum and orbital variability

The average optical spectrum of J1048 is dominated by absorption lines typical of late-type stars (for instance TiO λ6161, FeII λ6361 and FeII λ6494) and by the Hα emission line. It is displayed in Fig. 1 together with a K8 V template, HD 154712B, that well matches the absorption features (see also Sect. 4). One of the typical features of the accretion disk, the HeI λ6678 emission line, is absent. The Hα profile in the average spectrum is complex and appears single-peaked, casting doubts on the presence of an accretion disk.

Fig. 1. Average optical spectrum of PSR J1048+2339 (top spectrum) compared with a K8 V template (HD 154712B, bottom spectrum). Arbitrary vertical offsets have been applied for clarity. Main emission and absorption lines are indicated. The vertical black dashed line represents the rest wavelength of He I λ6678, which is absent in our spectra.

The Hα emission shows strong variations along the orbit (Fig. 2): at some orbital phases it disappears (for example at ϕ = 0.42), while at other epochs it appears with a single component (ϕ = 0.63) and sometimes it displays a weak secondary peak (ϕ = 0.84). Binary phases are computed using the extremely accurate radio pulsar ephemeris (see Sect. 6). Given the nearly circular orbit, the choice of the time of the ascending node T_asc is arbitrary and in the paper we opted to decrement that by P_b/4, where P_b is the orbital period, with respect to the one computed from radio ephemeris, so that phase 0 corresponds to the inferior conjunction of the companion star.

Fig. 2. Hα emission variability along the whole orbit. The orbit phase increases from left to right, from top to bottom. The vertical red dashed line represents the Hα λ6561 rest wavelength. We use the convention where ϕ = 0 is the point where the companion star is closer to the observer.

Differently from what was found by Strader et al. (2019), in our case the Hα emission is weak and asymmetric, with stronger blue or red components depending on the orbital phase. Near the superior conjunction (ϕ = 0.4–0.5), the Hα emission has almost completely vanished or is weakly detected as an absorption feature, while near the inferior conjunction (ϕ = 0) two peaks emerge (see Fig. 3, left plot). The blue component is weaker and associated with the companion star that is approaching the Earth. The red component vice versa is emitted by material receding from the Earth along our line-of-sight. This could be ascribed to the material flowing from the donor star toward the pulsar. We constructed the trailed spectrogram, shown in the right panel of Fig. 3, to investigate the orbital line variability in more detail. It displays the orbital evolution of the Hα emission line in 16 phase bins. Between phases 0.75 and 0.90, an emission spot with positive velocities that does not trace the orbital motion of the companion star is clearly visible. This is a first indication for the presence of material overflowing from the Roche lobe, confirmed by the Doppler tomography study (see Sect. 5). The gas from which the Hα emission originates together with the irradiated companion star are probably responsible for the radio eclipses of J1048 that cover about half of the orbit (Deneva et al. 2016, 2021).

Fig. 3. Left: a schematic panel of PSR J1048+2339 showing the orbital motion of the companion with Hα variability. The vertical red dashed line represents the Hα λ6561 rest wavelength. Right: trailed spectrum showing the orbital evolution of the Hα emission line. The spectra have been normalized to the continuum. Two cycles are shown for clarity.

Measurements of the equivalent widths (EWs) of the Hα emission line show a hint of modulation at the orbital period although at low statistical significance (∼1.9σ, see Fig. 4, bottom panel). The mean EW of Hα is −3.1 ± 1.3 Å. In the trailed spectra of Hα in Fig. 3 the emission is at maximum at orbital phases 0.1–0.2 and 0.6–0.7. For comparison, we measured the EWs of metallic lines (in the spectral ranges 5970–6291 Å and 6421–6522 Å). The EW curve of metallic lines shows almost a sinusoidal profile with a single peak significance at a 3.9σ, while the EW curve of Hα has a noisy profile. This might be due to the fact that the Hα line is contaminated by extended emission as hinted at by the trailed spectra (Fig. 3), while the metallic lines originate only from the companion star.

Fig. 4. Equivalent width of metallic (red dots) and Hα (blue triangles) absorption and emission lines, respectively, as a function of the orbital phase. Two cycles are shown for clarity.

The simultaneous R-band light curve (Fig. 5) shows an asymmetric sinusoidal-like profile with a minimum around ϕ =  0 and an average magnitude R ∼ 19.3 mag (Vega). The peak-to-peak amplitude of the light curve is of about 0.9 magnitude. The light curve shape and the magnitude variations along the orbit suggest that the optical emission of the system is dominated by the emission from the hot face of the companion, heated by the pulsar wind. The curve does not have a pure sinusoidal profile, the minimum is around phase zero as expected, while the maximum precedes phase 0.5.

Fig. 5. Phase-folded R-band light curve of PSR J1048+2339 companion star. The light curve shows a sinusoidal-like modulation with a peak-to-peak amplitude of ∼0.9 mag. The minimum is around phase zero as expected, while the maximum slightly precedes phase 0.5. Two cycles are shown for clarity.

Differently from what Linares et al. (2018) observed in the RB PSR J2215+5135, we do not find drastic changes in the spectral absorption lines of the irradiated and the cold face of the companion star (Fig. 6). In PSR J2215+5135, Balmer and Mg-I triplet absorption lines show strong EW variation as a function of the orbital phase due to the irradiation by the pulsar. The EWs of Balmer and Mg-I triplet absorption lines change by a factor ≃7 and ≃3, respectively (Linares et al. 2018). Indeed, optical photometry of PSR J2215+5135 revealed single humped modulation with magnitude variations of ∼1.5 mag along the orbit, typical of strongly irradiated system (Breton et al. 2013). In J1048, the EWs of Hα emission line and of metallic absorption lines vary by a small factor (≃3 and ≃2 respectively) along the orbit and we do not observe large temperature changes between the dark and bright sides of the companion star. The ∼0.9 mag variation in amplitude observed in the R-band light curve could therefore be attributed to a moderate pulsar irradiation.

Fig. 6. Main metallic absorption lines as observed at the inferior conjunction (ϕ ∼ 0.02) and at the superior conjunction (ϕ ∼ 0.55).

4. Radial velocity

Following Yap et al. (2019) we cross-correlated the 16 VLT-FORS2 spectra with six K V template stars in the spectral ranges 5970–6291 Å and 6421–6522 Å (taken from Casares et al. 1996), masking the main telluric and interstellar features and the Hα emission line, obtaining six radial velocity curves for each template star. We then performed a (least-squares) sine fit to our radial velocity data, allowing the projected radial velocity semi-amplitude (K₂) and the systemic velocity of the companion star (γ) to vary. We rescaled the error bars of the FORS2 measurements by a factor of 4.3, in order to obtain a reduced χ² = 1. The best match was obtained for the K2 V and K8 V templates, with equal significance (see Table 3).

We employed two independent approaches to classify the spectral type of the companion star. In the first one, we compared the ratio of the EW of the absorption lines in J1048 with different stellar templates. We considered the ratios λ6361(FeII)/λ6161 (TiO) and λ6494(FeII)/λ6161 (TiO). Using this criterion, the most likely spectral types are K7 V-K8 V (see Table 4).

In an alternative approach, we subtracted different templates appropriately degraded from the target average spectrum. We broadened stellar templates from 10 to 200 km s⁻¹, in steps of 10 km s⁻¹ and subtracted the resultant broadened spectra from the J1048 Doppler-corrected average spectrum in the spectral range 6432–6477 Å, using a limb darkening coefficient μ = 0.5. This subtraction method allowed us to estimate the companion star rotational broadening v sin i. In the optimal subtraction, the stellar template is scaled by a factor 0 ⩽ f_frac ⩽ 1, which represents the fraction of flux from the companion star. We used this quantitative approach to match the observed absorption lines from J1048 to set of templates used before. The best fit was obtained for the template K8 V with a star contribution of f_frac = 0.56 ± 0.02 and a projected rotational velocity for the companion of v sin i = 105 ± 15 km s⁻¹ (see Table 5). We verified that this result is not significantly affected by the variation of the limb darkening μ in the range of [0.0,0.8].

Based on the results of these three analyses, we conclude that the most likely spectral type for J1048 is a K8 V star (T_eff ∼ 4000 K). This determination is consistent with what was found by Yap et al. (2019) and has an uncertain of a couple of subspectral types. We found K_2, obs = 343.3 ± 4.4 km s⁻¹ and γ = −10.5 ± 3.1 km s⁻¹ (see Table 3 and Fig. 7). However, our constraint on K_2, obs does not seem to depend on the fine determination of the spectral type.

Fig. 7. Radial velocities of the companion star in PSR J1048+2339 computed using a K8 V template and folded on the ephemeris derived using SRT. The best-fitting sinusoidal radial velocity curve is superposed.

4.1. System parameters

Our determination of the mass ratio q and the NS mass relies on the determination of K₂ obtained by the measurement of the Doppler shifts of the absorption lines along the orbital period. This measurement is affected by a systematic error due to the pulsar irradiation, which determines a nonuniform absorption distribution on the stellar surface. The absorption lines formed on the companion photosphere are partially quenched in its irradiated face so that the center of light of these lines is shifted toward the outer face of the companion, causing a systematic shift in the velocity determined at each binary phase. Consequently, the radial velocity curve is not exactly sinusoidal and the K₂ obtained by fitting a sine curve is different from that corresponding to the center of mass of the donor. In our case, considering only absorption lines in the cross-correlation with stellar templates, the observed K₂ can be overestimated.

As shown in Wade & Horne (1988), we corrected the velocity using the so-called K-correction. The corrected observed radial velocity semi-amplitude is:

$$\begin{array}{c}\hfill K_{2,\mathrm{corr}}(q)=K_{2,\mathrm{obs}}-\mathrm{\Delta }K(q)=K_{2,\mathrm{obs}}-f_K(q)v\sin i.\end{array}$$(1)

Here, f_K is a geometrical correction factor smaller than unity (0 < f_K < 1). f_K = 0 indicates uniform emission over the donor star and thus no K-correction. If the nonirradiated hemisphere has a uniform absorption line strength and the other hemisphere does not contribute to the observed stellar absorption lines at all, f_K = 4/3π ≈ 0.42 (Wade & Horne 1988). Hence, applying the K-correction (Eq. (1)), we get K_2, corr = 298.7 ± 7.7 km s⁻¹. This is an unrealistic assumption in the case of J1048 that leads to obtain a firm lower limit of K₂, because both hemispheres contribute to the observed absorption lines (see Fig. 6). On the other hand, the measured value of K₂ = 343.3 ± 4.4 km s⁻¹, without K-correction, is an upper limit.

Using the NS mass function from Deneva et al. (2016), we derived the semi-amplitude of the pulsar projected radial velocity K₁ = (2π c x_p)/P_b = 72.7396 ± 0.0003 km s⁻¹, where x_p is the projected semi-major axis of the pulsar orbit and P_b is the orbital period. By combining K₁ and K₂ (291 < K₂ < 348 km s⁻¹, at 1σ c.l.), we found a mass ratio q = M₂/M_NS = K₁/K_2, corr of 0.209 < q < 0.250. This implies that the NS mass M_NS = $P_b{K}_2^3{(1+q)}^2/(2\pi G{\text{sin}}^{3}i)$ is (1.0 < M_NS < 1.6) sin⁻³ i M_⊙. Fixing the inclination to the upper limit of 76° due to the absence of X-ray eclipses (Yap et al. 2019), we constrained the NS mass to the range of 1.08 < M_NS < 1.76 M_⊙ and the mass of the companion M₂ = (K₁ M_NS)/K₂ to the range of 0.27 < M₂ < 0.37 M_⊙. We favor values of M_NS closer to the more massive scenario, because the irradiation seems to weakly affect the K₂ measurement. On the other hand, assuming the NS mass upper limit of ∼2.3 M_⊙ (Margalit & Metzger 2017), we can derive a firm lower limit inclination of 62° for a nonirradiated system (f_K = 0) and of 53° for a strongly irradiated system (f_K = 4/3 π).

We note that the values of K₂, M_NS and M₂ obtained here are systematically lower than those obtained from Strader et al. (2019). However, the data of Strader et al. (2019) are spread over four nonconsecutive nights (whereas our data were collected over the same night) and they are not adjusted for the K-correction. Without the K-correction, we get an upper limit of M_NS⩽1.76 M_⊙ (at 1σ c.l.), compatible within the errors with the value of Strader et al. (2019).

5. Hα Doppler tomography

We constructed a trailed spectrogram covering the 6250–6700 Å interval to compare the absorption and emission lines structures (see Fig. 8). The absorption features are visible throughout the orbital cycle and trace the companion star motion, while the Hα emission lines fade at some orbital phases (see Sect. 3). The derived Doppler image of the Hα emission line computed with the maximum entropy method (Marsh & Horne 1988) is shown in velocity coordinates in Fig. 9. The Roche lobe of the companion and the gravitational free-fall path of the gas stream are superposed over the map. Figure 9 shows the map constructed assuming an inclination i ∼ 76° and using the K₂ upper limit with M_NS ⩽ 1.76 M_⊙ and q ⩾ 0.209.

Fig. 8. Trailed spectra averaged in 16 orbital phase bins revealing the orbital evolution of several absorption features and Hα emission lines. Three orbital phases are shown for clarify.

Fig. 9. Doppler map of the Hα emission lines. We plot the gravitational free-fall gas stream trajectory (lower curve), the Keplerian disk velocity along the stream (upper curve) and the Roche lobe of the donor star using K1 = 72.7639 km s−1 (Deneva et al. 2016), i = 76° (Strader et al. 2019; Yap et al. 2019), MNS = 1.76 M⊙ and q = 0.209. The plus signs represent, from top to bottom, the companion and the NS, while the ‘×’ sign represents the center of mass. The map was computed for a systemic velocity γ = −10.5 km s−1.

This image cannot be interpreted as a time-average map because the data were acquired over one night and the system displays strong variability on very short timescale (few days; Wang et al. 2013; Yap et al. 2019). Therefore, this Doppler map represents the 2D velocity distributions of the Hα emission on March 19, 2020.

The map does not show evidence of an accretion disk, but shows a relatively faint spot emission coincident in phase and velocity with the companion star. The S-wave profile component from the companion is detected in the trailed spectra (see Fig. 3), although the emission is not continuous along the orbit. In addition, the Hα map reveals an extended bright spot in the upper left velocity quadrant close to the companion star, suggesting an intra-binary shock created from the interaction between the pulsar relativistic wind and matter leaving the companion star. The proximity of the intra-binary shock to the companion star indicates the presence of material in that region. If the secondary is filling its Roche-lobe the material streams away from L₁, otherwise the companion mass loss mechanism is driven by the ablation of the outer layers of the star by the energetic pulsar wind. Given the orbital separation and the mass ratio q, we can determine the approximate Roche-lobe radius R_L (Eggleton 1983). Assuming i ∼ 76° we obtain 0.50 R_⊙ < R_L < 0.55 R_⊙, comparable to the radius of a K8 V star (Pasinetti Fracassini et al. 2001). Yap et al. (2019) derived a Roche-lobe filling factor of 0.86±0.02 during a phase of single-peaked optical light curve, casting doubts on whether the secondary is truly filling its Roche-lobe. Given the uncertainty in K₂ and in the dimension of the companion star irradiated by the pulsar wind, our analysis does not favor one mechanism for the mass loss over the other. There is no evidence for emission from an accretion disk, although its formation could take place on a very rapid timescale, as in IGR J18245–2452 (Papitto et al. 2013).

6. Radio timing analysis

The raw data taken in baseband mode with the SRT in the P band, and with the LOFAR in the VHF band, were first coherently dedispersed (Bassa et al. 2017) at J1048’s nominal dispersion measure (DM) of 16.65 pc cm⁻³. The four Stokes parameters were summed to form total-intensity search-mode files. For the SRT P-band data, these retained a time resolution of 64 μs and 640, 125-KHz-wide frequency channels. For the LOFAR VHF-band data, the search-mode files kept a time resolution of 81.92 μs and 1600, 48-kHz-wide channels. As explained in Sect. 2.2, the SRT L-band data were already taken in search mode, so no prior manipulation was needed for them.

Once all the observations were available as search-mode files, we conducted the analysis as follows. For each observation, we first created a mask to filter out the strongest RFI present in the data, using the rfifind routine, part of the PRESTO¹⁴ pulsar searching package (Ransom et al. 2002). Taking the mask into account, we then used PRESTO’s prepdata to generate an RFI-free, frequency-summed dedispersed time series. These were folded using PRESTO’s prepfold and the ephemeris published by Deneva et al. (2016). However, it is well known that many RBs and BWs show a strong orbital variability (Shaifullah et al. 2016; Bak Nielsen et al. 2020; Hebbar et al. 2021), often causing an old ephemeris not to be able to predict the actual orbital phase at a future time. This results in a wrong correction for the orbital Rømer delay (Blandford & Teukolsky 1976) and, consequently, the pulsed emission can be missed when folding. Following Ridolfi et al. (2016), we overcame this issue by doing a brute-force search in the time of the ascending node, T_asc, using SPIDER_TWISTER¹⁵ and a step size of dT_asc = 4.655 s. We obtained a clear detection of J1048 in the L-band observation taken in August 2019 when using a T_asc MJD value of 58721.666409(54). This is ∼28 s later than the T_asc value that can be extrapolated from the Deneva et al. (2016) ephemeris. Folding P band data taken simultaneously with the same parameters, however, resulted in a nondetection.

The simultaneous SRT+LOFAR observations of J1048 taken in March 2020 were analyzed in the same way. The only difference was that, given the very long duration of these observations (which covered one full orbit), we also splitted them into 1-hr-long chunks and searched for J1048’s pulsations separately in each of them. Still, we were unable to detect J1048 in any of the three radio bands.

The nondetection in the L-band observation of March 2020 can likely be ascribed to interstellar scintillation, which is typically affecting pulsars with very low DMs (such as J1048). This effect randomly enhances or suppresses the observed flux density of the pulsar, with a dependence on the observing frequency. The confirmation that scintillation is indeed important for J1048 is given by the detection of the August 2019 L-band observation (Fig. 10): we can see that the J1048’s pulsed signal is boosted only in a very narrow band around ∼1.4 GHz, whereas it is barely detectable in the rest of the band. In addition to scintillation, another mechanism responsible for the nondetection in the P and VHF bands, can be the absorption due to the intra-binary eclipsing material. Partial or total eclipses seen in many BW and RB systems (Rasio et al. 1989) can be ascribed to the frequency-dependent optical depth of the eclipsing material. The dense layers of ionized material ablated from the donor star are less transparent at low frequencies, determining eclipses in the P and VHF bands (Polzin et al. 2018). As a side effect, at low frequencies, the scattering (Lorimer & Kramer 2004) can smear the individual pulses, making them broader than the pulsar spin period, but not preventing their detection.

Fig. 10. Intensity as a function of pulse phase (x-axis) and frequency (y-axis) for the SRT observation performed in the L band on 19 August, 2019. Signal strength is not constant at all frequencies. At ∼1.4 GHz there is a peak of intensity, while below 1.3 GHz and over 1.8 GHz the pulsed signal is almost absent.

From the nondetections, we can calculate an upper limit on the mean flux density, S^mean, that J1048 must have had at the time of the observations. To do so, we use the modified radiometer equation (Manchester et al. 1996; Lorimer & Kramer 2004)

$$\begin{array}{c}\hfill S^{\mathrm{mean}}=\beta {\displaystyle \frac{(S/N)T_{\mathrm{sys}}}{G\sqrt{n_p{t}_{\mathrm{int}}\mathrm{\Delta }\nu }}}\sqrt{{\displaystyle \frac{\delta }{1-\delta }}},\end{array}$$(2)

where β ≳ 1 is a correction factor due to digitization; n_p is the number of polarizations summed; S/N is the pulse signal-to-noise ratio; T_sys is the system noise temperature; G is the telescope gain; t_int is the integration time; Δν is the bandwidth; δ is the pulse duty cycle. Table 6 lists the values of each parameter for all the radio observations that we used in the S^mean estimations. For the August 2019 P band nondetection we estimate $S_P^{\mathrm{mean}}\lesssim 0.57$ mJy, whereas in the case of the March 2020 nondetections we find $S_L^{\mathrm{mean}}\lesssim 0.05$ mJy, $S_P^{\mathrm{mean}}<0.24$ mJy and $S_{\mathit{VHF}}^{\mathrm{mean}}\lesssim 9.21$ mJy for the P, L and VHF observations, respectively. For the only radio detection available to us (August 2019 L-band observation) we measure $S_L^{\mathrm{mean}}\sim 1.22$ mJy.

7. Optical timing analysis

We tested the acquired data with IFI+Iqueye for pulsations folding all three observations together with the spin period of the pulsar (0.00466516294153(3) s, extrapolated from Deneva et al. 2016 at MJD 58927.0) and 10 phase bins. The time of passage at the ascending node was fixed at T_asc = 58721.666409 MJD (see Sect. 6). A χ² test with a constant term corresponding to the average count rate was applied to the folded profile. The χ² returns a value of 11.5 for 9 degrees of freedom, corresponding to a false positive probability of 24% and a significance of pulsed emission of only 1.2σ. Different choices for the number of phase bins provide similar results.

A conservative upper limit to the optical pulsed emission is derived from the power spectra of the three observations. The maximum excess power in the 0.5–500 Hz frequency range with a time bin of 1 ms corresponds to a root mean square fractional variability of ∼0.1% on the non background subtracted light curve or ∼2 counts s⁻¹ for the observed average background count rate. This upper limit refers to a single trial for the time of passage at the ascending node.

We performed also an epoch folding search varying the value of T_asc for our epoch in an interval between −250 s and +250 s around the value reported above and in steps of 1 s, and folding the event lists with 16 bins per period. No significant signal is detected. The maximum of the χ² is ∼42 for 15 degrees of freedom, corresponding to a significance of only 1.6σ (for 501 trials).

8. Discussion and conclusions

The multiwavelength observations acquired in March 2020 confirmed the extreme emission variability of J1048 detected by Deneva et al. (2016), Cho et al. (2018), Strader et al. (2019) and Yap et al. (2019). We observed rapid and strong emission changes in the Hα emission line along the orbit (Fig. 2). Using the optical photometry we revealed an average R ∼ 19.3 mag counterpart with orbital variability typical of moderate irradiated systems (Fig. 5). The shape of the R-band light curve showed that the modulation was dominated by the pulsar heating and not by ellipsoidal variations. We measured with radio timing analysis a strong orbital variation as in other RB systems (Jaodand et al. 2016; Ridolfi et al. 2016); the T_asc lagged the predicted value of Deneva et al. (2016) ephemeris by ∼28 s.

We measured the radial velocity curve through the cross-correlation with late-type template stars and we set the following constraints to the mass ratio, to the mass of the NS and of the companion star: 0.209 < q < 0.250, $(1.0<M_{\mathrm{NS}}<1.6){\sin }^{-3}i{M}_{\odot }$ and (0.24 < M₂ < 0.33) sin⁻³i M_⊙. Assuming the inclination upper limit of 76° (Yap et al. 2019), we derived a NS mass of 1.08 < M_NS < 1.76 M_⊙. Given the uncertainty on the inclination, we cannot exclude the presence of a high-mass NS. Assuming the NS mass upper limit of ∼2.3 M_⊙, the inclination of J1048 must be > 62° if the system is nonirradiated, and > 53° if it is strongly irradiated. We constrained the spectral type of the donor to a K8 V star (∼4000 K). The companion star contributes 56% to the total light observed at 6400Å, suggesting the presence of an additional diluting source such as an intra-binary shock.

We presented the first Hα Doppler tomography of a RB system in a rotation-powered, radio pulsar state. Doppler tomograms of RBs were only performed for the two transitional systems PSR J1023+0038 and XSS J12270-4859, but during the disk state (Hakala & Kajava 2018; de Martino et al. 2014). The reconstructed Hα Doppler map revealed a significant emission close to the inner Lagrangian point along the gravitational free-fall trajectory of the gas and at the companion star surface (see Fig. 9). The extended structure seen in the Hα line is a first and clear evidence of an intra-binary shock emission caused by the interaction of the pulsar wind with material lost from the companion star. Such shocks have been earlier proposed to be responsible for the X-ray orbital variability in MSP binaries (Bogdanov et al. 2005) and their presence has been indirectly confirmed by several observations in different domains (Wadiasingh et al. 2018). Our detection of a localized emission of Hα line makes J1048 the first RB showing direct evidence of an intra-binary shock emitting at other wavelengths. Whether the secondary star is overflowing its Roche lobe, the observed emission could be produced by the interaction of the gas streaming the L₁ and the pulsar wind. Due to the lack of strong constraints on Roche-lobe filling factor (see Yap et al. 2019), the shock emission could also be due to the interaction of ablated material from the companion and the strong pulsar wind. In this last configuration, and given that the Hα emission comes close to L₁, this could be a signature of matter from the intra-binary shock that are ducted to the companion surface along the magnetic field of the companion star as hypothesized by Sanchez & Romani (2017). The fact that this system shows variability on different epochs may be compatible with this scenario.

The presence of the intra-binary shock emission can explain the prolonged radio eclipses (Deneva et al. 2016, 2021) and the sporadic appearance of Hα double-peaked emission lines (Strader et al. 2019) observed in J1048. This map confirms the hypothesized scenario in which RBs are characterized by a large amount of material being stripped off the companion star due to the pulsar ablation process (Breton et al. 2013). The intra-binary shock in J1048 appears to be closer to the companion than to the pulsar. This evidence contrasts with the recently proposed scenario for the X-ray shock emission in RB systems (Wadiasingh et al. 2018; van der Merwe et al. 2020) in which the intra-binary shock is assumed to be closer to the pulsar. This assumption is supported by the large radio eclipse fractions observed in RBs at pulsar superior conjunction (Archibald et al. 2009; Roy et al. 2015; Deneva et al. 2016), as well as the double-peaked X-ray orbital modulation centered at the pulsar inferior conjunction (Romani & Sanchez 2016; Al Noori et al. 2018; de Martino et al. 2020). However, the Doppler map of J1048 suggests that X-ray and Hα emissions do not trace the same region of the intra-binary shock material. Hα emission may trace a portion of intra-binary shock closer to the companion while the X-rays could originate closer to the NS.

The absence of the pulsed radio signal in March 2020 observations cannot be ascribed to a state transition such as the one observed in the tMSPs (Archibald et al. 2009; Papitto et al. 2013), but rather to scintillation or absorption phenomena (Gedalin & Eichler 1993; Thompson et al. 1994; Roy et al. 2015). At very low frequency (0.1–0.2 GHz) the radio signal could be totally eclipsed by absorption due to the presence of a low-density highly-ionized gas cloud spilling off the companion (Broderick et al. 2016). Indeed, no increase in the X-ray and optical fluxes were observed and the Hα Doppler map clearly showed the presence of an intra-binary shock and not an accretion disk. The presence of an extended shock front could lead in the near future to the formation of an accretion disk, but as we have observed in the three tMSPs discovered so far, the state transition is not predictable and takes place in a very rapid (few weeks) timescale (Archibald et al. 2009; Papitto et al. 2013; Bassa et al. 2014).

Although J1048 is one of only a handful of RBs that occasionally displays emission lines in their optical spectra, similar behavior is observed in other RBs and BWs, such as 3FGL J0838.8–2829 (Halpern et al. 2017a,b), PSR J1628–3205 (Cho et al. 2018; Strader et al. 2019), and PSR J1311–3430 (Romanova & Owocki 2015). Further spectroscopic observations of these RBs and BWs therefore seem warranted. It would be interesting to monitor these systems for several contiguous orbital phases to study the evolution of the Doppler map and confirm the rapid variability of the gas distribution.

Acknowledgments

We thank the referee for very useful comments and suggestions. A.M.Z. acknowledges the support of the PHAROS COST Action (CA16214). A.M.Z. would like to thank G. Benevento for comments on draft and for discussions. We thank Michele Fiori and Giampiero Naletto (part of the Aqueye+Iqueye team) for the support and service activities carried out within the framework of this project. We thank Jorge Casares for providing the stellar templates spectra. S.C. and P.D.A. acknowledge support from ASI Grant I/004/11/3. A.R. gratefully acknowledges financial support by the research Grant “iPeska” (P.I. Andrea Possenti) funded under the INAF national call PRIN-SKA/CTA approved with the Presidential Decree 70/2016. F.C.Z. is supported by a Juan de la Cierva fellowship. A.P. and F.C.Z. acknowledge the International Space Science Institute (ISSI- Beijing), which funded and hosted the international team ‘Understanding and Unifying the Gamma-rays Emitting Scenarios in High Mass and Low Mass X-ray Binaries’. 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 PRIN-SKA/CTA: Toward the SKA and CTA era” (PI: Giroletti) Presidential Decree 70/2016; and partial support from PHAROS COST Action (CA16214). T.M.D acknowledges support from the Spanish Ministerio de Ciencia e Innovación under grant AYA2017-83216-P and the Ramón y Cajal Fellowship RYC-2015-18148. L.Z. acknowledges 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 and from INAF “Sostegno alla ricerca scientifica main streams dell’INAF” Presidential Decree 43/2018. The Sardinia Radio Telescope is funded by the Department of Universities and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS), and is operated as a National Facility by the National Institute for Astrophysics (INAF). This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. The results presented in this paper are based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme ID 0104.D-0589(A), on the observations collected at the Galileo telescope (Asiago, Italy) of the University of Padova and on data obtained with the International LOFAR Telescope (ILT) under project code DDT13_001. LOFAR is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefitted from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK.

References

All Tables

All Figures

	Fig. 1. Average optical spectrum of PSR J1048+2339 (top spectrum) compared with a K8 V template (HD 154712B, bottom spectrum). Arbitrary vertical offsets have been applied for clarity. Main emission and absorption lines are indicated. The vertical black dashed line represents the rest wavelength of He I λ6678, which is absent in our spectra.
In the text

	Fig. 2. Hα emission variability along the whole orbit. The orbit phase increases from left to right, from top to bottom. The vertical red dashed line represents the Hα λ6561 rest wavelength. We use the convention where ϕ = 0 is the point where the companion star is closer to the observer.
In the text

	Fig. 3. Left: a schematic panel of PSR J1048+2339 showing the orbital motion of the companion with Hα variability. The vertical red dashed line represents the Hα λ6561 rest wavelength. Right: trailed spectrum showing the orbital evolution of the Hα emission line. The spectra have been normalized to the continuum. Two cycles are shown for clarity.
In the text

	Fig. 4. Equivalent width of metallic (red dots) and Hα (blue triangles) absorption and emission lines, respectively, as a function of the orbital phase. Two cycles are shown for clarity.
In the text

	Fig. 5. Phase-folded R-band light curve of PSR J1048+2339 companion star. The light curve shows a sinusoidal-like modulation with a peak-to-peak amplitude of ∼0.9 mag. The minimum is around phase zero as expected, while the maximum slightly precedes phase 0.5. Two cycles are shown for clarity.
In the text

	Fig. 6. Main metallic absorption lines as observed at the inferior conjunction (ϕ ∼ 0.02) and at the superior conjunction (ϕ ∼ 0.55).
In the text

	Fig. 7. Radial velocities of the companion star in PSR J1048+2339 computed using a K8 V template and folded on the ephemeris derived using SRT. The best-fitting sinusoidal radial velocity curve is superposed.
In the text

	Fig. 8. Trailed spectra averaged in 16 orbital phase bins revealing the orbital evolution of several absorption features and Hα emission lines. Three orbital phases are shown for clarify.
In the text

Fig. 9. Doppler map of the Hα emission lines. We plot the gravitational free-fall gas stream trajectory (lower curve), the Keplerian disk velocity along the stream (upper curve) and the Roche lobe of the donor star using K1 = 72.7639 km s−1 (Deneva et al. 2016), i = 76° (Strader et al. 2019; Yap et al. 2019), MNS = 1.76 M⊙ and q = 0.209. The plus signs represent, from top to bottom, the companion and the NS, while the ‘×’ sign represents the center of mass. The map was computed for a systemic velocity γ = −10.5 km s−1.

In the text

	Fig. 10. Intensity as a function of pulse phase (x-axis) and frequency (y-axis) for the SRT observation performed in the L band on 19 August, 2019. Signal strength is not constant at all frequencies. At ∼1.4 GHz there is a peak of intensity, while below 1.3 GHz and over 1.8 GHz the pulsed signal is almost absent.
In the text