© The Authors 2024

1. Introduction

To date, there are over 550 pulsars¹ with spin periods ≲30 ms. Such rapidly rotating neutron stars (NSs) are named millisecond pulsars (MSPs). The commonly accepted ‘recycling’ scenario of their formation (Bisnovatyi-Kogan & Komberg 1974; Alpar et al. 1982) suggests that these NSs are spun up in binary systems through accretion from a secondary star, which eventually becomes a low-mass white dwarf (Tauris 2011). This also explains why most MSPs are found in binary systems. However, not all binary companions to MSPs are white dwarfs: some of them tend to have a more exotic nature. Those binary MSPs that contain low-mass non- or semi-degenerate companions in tight orbits with periods of P_b ≲ 1 d are called ‘spider’ pulsars. This class includes redbacks (RBs) and black widows (BWs) (see, e.g. Roberts 2013). In these systems, one side of a spider companion is irradiated by the energetic pulsar wind, which leads to heating and evaporation of its stellar material. The RBs have more massive companions with masses of 0.1–1 M_⊙, while masses of BW companions are ≲0.05 M_⊙. The formation of such systems as well as the link between the two subclasses are matters of debate (Chen et al. 2013; Benvenuto et al. 2014; Ablimit 2019; Guo et al. 2022).

The number of known spiders has increased significantly in the last few years. Many of them were discovered thanks to the multi-wavelength studies of unidentified γ-ray sources detected with the Fermi observatory (e.g. Salvetti et al. 2017; Swihart et al. 2022). In addition to a broad range of spiders’ fundamental parameters commonly derived from radio follow-up searches, optical observations can also have a significant impact. In particular, they allow one to derive the temperature distribution over a companion surface, its Roche-lobe filling factor, irradiation efficiency, distance to the binary, and masses of its components (Strader et al. 2019; Draghis et al. 2019; Swihart et al. 2022). With a larger sample of newly detected spider systems, further investigations can help to clarify their origin and evolution.

Two BW pulsars, PSRs J2017−1614 and J1513−2550 (hereafter J2017 and J1513), were discovered in the radio with the Green Bank Telescope during the search for counterparts of unassociated Fermi sources (Sanpa-arsa 2016). Their parameters are listed in Table 1. The radio emission of both binary MSPs was eclipsed, but only in some observational runs. After the discovery in the radio, their periodic pulsations were also found in γ-rays (Sanpa-arsa 2016), but no studies in X-rays have been reported.

The optical counterpart of J2017 was first observed in the R band with the 2.4-m telescope of the Michigan-Dartmouth-MIT (MDM) observatory (Sanpa-arsa 2016). It shows significant magnitude variations, from R ≈ 21.8 mag at the maximum brightness to R ≳ 24 mag at the minimum, with a period of ≈2.3 h, which is compatible with the binary period derived from the radio data (Sanpa-arsa 2016).

A possible optical companion of J1513 could be preliminarily identified by positional coincidence with the star PSO J228.3472−25.8420 from the Panoramic Survey Telescope and Rapid Response System DR 2 catalogue (Pan-STARRS, Flewelling et al. 2020). Its mean brightness in the i-band is about 22 mag.

No detailed optical studies of J1513 and J2017 have been reported so far. Here, we present the results of the first multi-band time-series optical photometry of J2017 and J1513 and analyse X-ray data for both BWs obtained with the XMM-Newton observatory. Optical observations and data reduction are described in Sec. 2 and analysed in Sec. 3, X-ray data are presented in Sec. 4, and the results are discussed and concluded in Secs. 5 and 6.

2. Optical observations, data reduction, and calibration

2.1. J2017

The phase-resolved photometric observations of the J2017 field² were performed on September 7 2021 with the HiPERCAM instrument³ (Dhillon et al. 2016, 2018, 2021) at the 10.4-metre Gran Telescopio Canarias (GTC). The HiPERCAM allows for simultaneous imaging in five (u_s, g_s, r_s, i_s, and z_s) high-throughput ‘Super’ SDSS filters.

In each band, we obtained 79 images with individual exposure times of 140-s. The field of view (FoV) is 2.​​′8 × 1.​​′4, with the 0.​​″081 image scale in the 1 × 1 binning mode. The weather conditions during the observations were photometric. The log of observations is presented in Table 2.

We reduced the data following the HiPERCAM pipeline manual⁴, including bias subtraction and flat field correction (Dhillon et al. 2016). We also produced bad pixel correction with the defect files obtained from the HiPERCAM website. In addition, we applied the defringing procedure to the z_s-band images. An example of an individual image in the r_s-band is presented in Fig. 1, left panel. The pulsar counterpart position is marked by the blue box. Its variability is demonstrated in the insets.

Fig. 1. Optical images of the fields of the two pulsars. Left panel: 2.​​′8 × 1.​​′4 FoV of the HiPERCAM instrument in the rs band containing J2017. Right panel: 2.​​′7 × 2.​​′7 image of the J1513 field in the r′ band obtained with the Magellan-1 telescope. The pulsars’ vicinities within the small boxes are enlarged in the insets to demonstrate the maximum (bottom) and minimum (top) brightness phases of their companions, the positions of which are marked by blue circles.

Using the optimal extraction algorithm (Naylor 1998), we measured instrumental magnitudes of the counterpart and a dozen of bright stars in the FoV whose magnitudes are available in the Pan-STARRS catalogue. To avoid centroiding problems of the companion during its faint brightness stages, its position was fixed relative to a nearby star when the source was in its maximum brightness phase. We performed photometric calibration using the GTC atmospheric extinction coefficients k_us = 0.48, k_gs = 0.17, k_rs = 0.1, k_is = 0.05, and k_zs = 0.05 and the spectro-photometric standards WD1606+422 and WD2047+372⁵ (Dhillon et al. 2021) observed during the same night as the target. To eliminate possible systematic errors and to account for photometric zero-point variations during the observation, we used the technique described in Honeycutt (1992) and the Pan-STARRS stars mentioned above. The measured magnitudes of the latter were compared with their catalogue values applying the photometric transformation equations, which establish relations between the magnitudes in the HiPERCAM and Pan-STARRS bands (Tonry et al. 2012). This allowed us to verify the calibration, which resulted in the final zero-points z_us = 27.30(5) mag, z_gs = 27.80(5) mag, z_rs = 27.65(5) mag, z_is = 27.41(5) mag, and z_zs = 27.06(5) mag. The detection limits varied in the ranges of u_s = 24.9–25.2 mag, g_s = 23.5–25.0 mag, r_s = 23.3–24.5 mag, i_s = 22.7–24.2 mag, and z_s = 23.0–24.3 mag.

In the u_s-band, at given individual exposures the J2017 companion brightness was below the HiPERCAM sensitivity limit. As a result, we obtained J2017 light curves only in the g_s, r_s, i_s, and z_s bands.

2.2. J1513

Time-series optical photometry of J1513 was carried out using the 2.1-m telescope of Observatorio Astronómico Nacional San Pedro Mártir (OAN-SPM) and the 6.5-metre Magellan-1 telescope (Baade).

OAN-SPM observations were performed in the B, V, R, and I bands with the ‘Rueda Italiana’ instrument during three observing runs in May 2019. The detector FoV is 6′×6′ with a pixel scale of 0.​​″34 in the 2 × 2 CCD pixel binning mode. Magellan-1 observations were carried out in the Sloan r′ band with the Inamori-Magellan Areal Camera and Spectrograph (IMACS) during one observing run in April 2022 (PI R. E. Mennickent). The pixel scale was 0.​​″111 in the 1 × 1 CCD pixel binning mode. The target was exposed on chip 3 with the FoV of 3.​​′8 × 7.​​′5. The weather conditions were photometric during both observing runs. The log of all observations is given in Table 2.

Standard data reduction, including bias subtraction, flat-fielding, and cosmic ray rejection, was performed for each dataset utilising the Image Reduction and Analysis Facility (IRAF) package. The defringing procedure was applied to the OAN-SPM images in the I-band.

Astrometric referencing was performed using 12 stars from the Gaia DR 2 astrometric catalogue (Gaia Collaboration 2016, 2018) Formal rms uncertainties of the resulting astrometric fit were Δα ≲ 0.​​″08 and Δδ ≲ 0.​​″08 for the Magellan-1 and Δα ≲ 0.​​″04 and Δδ ≲ 0.​​″04 for the OAN-SPM.

An example of an individual 50 s r′-band image of the J1513 field obtained with the Magellan-1 telescope is shown in the right panel of Fig. 1. The co-ordinates of the J1513 detected companion candidate (marked by the blue circle in the inserts) are fully consistent with the pulsar radio position given in Table 1, confirming its preliminary identification based on the Pan-STARRS catalogue.

Photometric calibrations were performed using photometric standards observed during the same nights as the target. For the Magellan-1 data, we used three standards from the PG 1525-071 field (Landolt 1992) and the atmospheric extinction coefficient k_r′ = 0.1 provided by the IMACS team⁶. As a result, the zero point in the r′-band was 28.15(1) and the 3σ sensitivity level was ≈24.8 in the AB system. For the OAN-SPM data, we used the PG1323 spectrophotometric standard and the extinction coefficients k_B = 0.25 mag, k_V = 0.14 mag, k_R = 0.07 mag, and k_I = 0.06 mag relevant to the telescope site. The resulting zero points were z_R = 24.96(2) mag, z_B = 24.91(3) mag, z_V = 25.16(2) mag, and z_I = 24.15(2) mag, and the 3σ sensitivity levels were ≈23.6 mag, ≈24.1 mag, ≈23.8 mag, and ≈22.0 mag, respectively⁷.

We performed standard aperture photometry of the source and a number of stars from the Pan-STARRS catalogue in each image. Aperture radii were chosen according to the optimal extraction algorithm (Naylor 1998). The measured magnitudes of the stars from the catalogue remained consistent with their catalogue values during the observations.

3. Light curve analysis

The J2017 and J1513 optical companions’ barycentre-corrected light curves folded with their known orbital periods (see Table 1) are presented in Fig. 2. The times of the ascending node were chosen as the reference times for each band. The curves show one broad peak per period and peak-to-peak amplitudes of ≳2 mag, as was expected for BW systems in which the heating of the companion by the pulsar is dominated over the ellipsoidal modulation (e.g. Draghis et al. 2019; Mata Sánchez et al. 2023). This firmly establishes both optical sources as the true companions of the considered binary MSPs.

Fig. 2. Light curves of the two pulsars and results of their modelling. Top panel: gs, rs, is, and zs band light curves of the J2017 system folded with the orbital period and the best-fitting model (solid lines). Two periods are shown for clarity. Bottom panel: J1513 light curves folded with the orbital period. Data obtained in different filters and different nights are shown with various colours and symbols, indicated in the legend. All magnitudes are in the AB system. The best-fitting model is shown by the solid line. Shaded areas mark the part of the light curves where small-scale wavering was detected (see text for details). In both panels, points with dashed error bars lie below the 3σ detection limits. The phases ϕ = 0.0 correspond to the minima of the systems’ model brightness. Sub-panels show residuals calculated as the difference between the observed (O) and calculated (C) magnitudes for each data point in terms of the magnitude error, σ.

The light curves of J2017 are symmetric relative to the brightness maxima where the colour indices are g_s − r_s = 0.7, r_s − i_s = 0.4, and i_s − z_s = 0.2. The colours slowly become redder as the system rotates towards the minimum, reflecting the difference between the ‘day-side’ and ‘night-side’ temperatures of the pulsar companion. When the system approaches the minimum, its brightness drops below the detection limit.

The J1513 light curves are more complex. They demonstrate significant asymmetry with a steeper initial brightness increase and a more gradual fading to the minimum. We also note apparent regular brightness variations (wavering) on a short timescale with a peak-to-peak amplitude of ∼0.1 mag, which are clearly seen in the r′-band at the rising part of the light curve within the orbital phase range of 0.1–0.6. They are also present in the R-band curve, though with smaller significance due to the worse time resolution and sensitivity. The variations are not visible in the phase range of 0.6–1.0 in either of the two bands. To check whether they are caused by some atmospheric effects, we investigated light curves of several isolated stars with magnitudes within the range of r = 17.8–20 mag in the J1513 field. We did not find similar patterns there. For this reason, we concluded that the observed features are real and can be attributed to some physical process in the J1513 system.

To derive parameters of the systems, we fitted the obtained light curves applying a binary model consisting of an NS and a low-mass companion heated by the pulsar wind. The model and fit method are described in detail in Zharikov et al. (2013, 2019). It is assumed that each surface element of the companion emits a black-body spectrum whose effective temperature varies from element to element. This direct heating model produces symmetrical light curves in respect to their brightness maxima.

In the case of J2017, we used only g_s-, r_s-, and i_s-band light curves because the black-body approximation does not describe well the observed spectral energy distribution in the wide spectral range from the g_s- to z_s band. As for J1513, since the utilised model does not include components that can describe the asymmetry and small-scale variations in the light curves, we fitted only the smooth part of the r′-band light curve within the orbital phase range of ϕ ∈ [0.6, 1.0]. We did not include the other bands in the fitting procedure as their respective light curves contain too few data points, and some of them cover phases with variations. For both objects, we excluded points that lie below the 3σ detection limits from the modelling.

The fitted parameters are the interstellar reddening, E(B − V), the distance, D, the pulsar mass, M_p, the orbit inclination angle, i, the effective irradiation factor, K_irr, which defines heating of the companion, the companion Roche lobe filling factor, f, and the companion ‘night-side’ temperature, T_n. The mass of the companion, M_c, or the mass ratio, q = M_c/M_p, was defined by the mass function (Table 1), which links M_p, M_c, and i. The allowed range of E(B − V) was limited by the total galactic interstellar absorption in the directions of the pulsars (≈0.12 mag, Green et al. 2019). The masses of the pulsars started from 1.4 M_⊙, and it was assumed that they cannot exceed 3 M_⊙. Other parameters did not have any special restrictions. They were started from random values inside the natural physical limitations. The gradient descent method was used to find the minimum of the χ² function. The resulting parameters for both pulsars are presented in Table 3 and the best-fitting models are shown by solid lines in Fig. 2.

The light curves of the J2017 are perfectly reproduced by the model with reduced χ² = 59.6 for 99 degrees of freedom (d.o.f.). The fit residuals are within the 2σ data errors (Fig. 2, top) and the fit parameters are consistent with those expected for BWs. In contrast, J1513 shows light curves different from most spider companions (see discussion below). Formally, the fit of a part of the light curve is acceptable, with χ² = 69.4 for 44 d.o.f.. However, as can be seen from the bottom panel of Fig. 2, the observed light curve demonstrates a strong emission excess over the symmetric model within the range of ϕ ∼ 0.1–0.35, indicating the presence of some additional source of radiation at these orbital phases. Moreover, there are small-scale brightness variations mentioned above that deserve a separate analysis.

To study the variations, it is natural to filter out the large-scale variations related with the orbital motion. To do that, we fitted the light curve at ϕ ∈ [0.1, 0.6] in the r′-band with a low order Chebyshev polynomial (Chebyshev 1853) and subtracted this approximation from the light curve. The result is presented in the top panel of Fig. 3. The filtered light curve demonstrates an apparent periodic variability for about 75 min between orbital phases of ∼0.15 and 0.45. We used the data points from this interval and the Lomb-Scargle periodogram method (Lomb 1976; Scargle 1982) provided by the astropy package⁸ to search for the brightness periodicity. The highest peak in the power spectrum corresponds to the period P_V ≈ 15.3 min (Fig. 3, middle panel). To estimate the power spectrum uncertainty, we used an equivalent of the bootstrap method. We simulated 10 000 light curves randomly scattered within the brightness intervals defined by ±1σ errors of the data, assuming a normal distribution. We obtained periodograms for each of the simulated curves. Then the periodogram distribution for each frequency bin was fitted with the Gaussian whose maximum position and width values were attributed to the most probable power spectrum value and its error in the bin, respectively.

Fig. 3. Result of the small-scale brightness variations analysis. Top panel: J1513 filtered light curve in the r′ band (black points with error bars). Points within shaded areas are excluded from the periodicity search. The red line is the best fit of the data by a sine function with the period of 15.3 min found from the periodicity search. The sub-panel shows fit residuals. Middle panel: power spectrum of the filtered light curve (the solid line) with its uncertainties (vertical error bars). The largest peak corresponds to the 15.3 min period. Bottom panel: light curve folded with this period and fitted with the sine function of the same period. The sub-panel shows fit residuals.

The calculated false alarm probability of the 15.3 min peak is ≈0.46, making it questionable. However, we constructed the power spectra of a non-variable reference star and the window function of the data for the same time range as for J1513. We did not find any signature of the 15.3-min period there. This indicates that the 15.3-min peak in J1513 is not an artefact.

In addition, we fitted the data with a constant and then with a sine function of the 15.3-min period plus a constant. The latter resulted in χ² = 24.1 for 30 d.o.f. and the average peak-to-peak amplitude of ≈0.08 mag (see Fig. 3, top panel). With the constant model, we obtained a slightly worse but also acceptable fit with χ² = 33.3 for 32 d.o.f..

The F-test applied to compare the two models gave a probability of ≈8 × 10⁻³ of improvement by chance of the sine fit with respect to the constant one, implying strong evidence against the constant model. The light curve folded with the obtained period shown in the bottom panel of Fig. 3, together with the best sine fit, appear to favour the presence of the periodicity. The remaining part of the filtered light curve shows only stochastic fluctuations of the same amplitude and no signs of the periodicity. Nevertheless, we finally conclude that while the small-scale variability is clearly visible in the considered phase interval, its apparent periodicity has a low statistical significance of ≈2.5σ and new observations are needed to confirm it.

4. X-ray data

We found X-ray counterparts to J2017 and J1513, 4XMM J201746.0−161416 and 4XMM J151323.3−255029, in the XMM-Newton Serendipitous Source Catalogue Data Release 13 (4XMM-DR13; Webb et al. 2020). According to the catalogue, 4XMM J151323.3−255029 is very weak, with an absorbed flux f_X = (5.8 ± 2.2)×10⁻¹⁵ erg s⁻¹ cm⁻² in the 0.2–12 keV band. 4XMM J201746.0−161416 is brighter, with f_X = (1.8 ± 0.3)×10⁻¹⁴ erg s⁻¹ cm⁻². Below, we present an independent analysis of the observations.

For both pulsars, the European Photon Imaging Camera-Metal Oxide Semiconductor (EPIC-MOS) detectors and EPIC-pn (PN hereafter) were operated in the full frame mode with the medium filter. The XMM-Newton Science Analysis Software (XMM-SAS) v.21.0.0 was used for data reduction. We reprocessed the data, applying the emproc and epproc tasks.

4.1. J2017

The XMM-Newton observation of the J2017 field was performed on October 15 2016 (ObsID 0784770601, PI M. Roberts, duration = 25 ks). Examination of the high-energy light curves extracted from the FoVs of EPIC detectors revealed some strong background flares. We used the espfilt task to filter event lists, which resulted in effective exposure times of 17, 18, and 13 ks for the MOS1, MOS2, and PN instruments, respectively.

The PN image of the pulsar field is presented on the left of Fig. 4. The X-ray source is clearly seen in the pulsar radio position. We extracted its spectra using the 14″-radius aperture. The background spectra were extracted from the source-free circle region with a radius of 28″. The resulting net counts in the 0.2–10 keV bands are 27, 29, and 72 for MOS1, MOS2, and PN data, respectively. The redistribution matrix and ancillary response files were generated by the rmfgen and arfgen routines.

Fig. 4. 3.​​′5 × 3.​​′5 XMM-Newton/PN images of the J2017 (left) and J1513 (right) fields in the 0.2–2 keV band. The pulsars’ radio positions (Table 1) are marked by the ‘X’ symbols. Dashed circles show regions used for the background extraction.

We grouped the spectra to ensure at least 1 count per energy bin and fitted them simultaneously in the 0.2–10 keV band in the X-Ray Spectral Fitting Package (XSPEC) v.12.13.1 (Arnaud et al. 2023). We used the W-statistics appropriate for Poisson data with Poisson background and the power law (PL) model that is usually utilised to describe X-ray spectra of spider pulsars (Swihart et al. 2022; Strader et al. 2019). To account for the interstellar absorption, we chose the tbabs model with the wilm abundances (Wilms et al. 2000). The equivalent hydrogen column density, N_H, was fixed at 8 × 10²⁰ cm⁻². The latter was calculated using the reddening value obtained from the optical light curves modelling (Table 3) and the relation from Foight et al. (2016). As a result, we obtained the photon index Γ = 2.7 ± 0.2, the unabsorbed flux in the 0.5–10 keV band F_X = (1.5 ± 0.2)×10⁻¹⁴ erg s⁻¹ cm⁻², and W = 80 for 81 degrees of freedom.

4.2. J1513

The J1513 field was observed with XMM-Newton on September 1 2016 (ObsID 0784770701, PI M. Roberts, duration = 25 ks). We detected only one short flare, which is clearly seen only in the PN light curve. Thus, we filtered the PN data using the count rate threshold of 0.5 counts s⁻¹, which resulted in the effective exposure time of 15 ks.

The J1513 X-ray counterpart is very faint (Fig. 4, right) and most prominent in the images from the PN camera, which is more efficient in the soft (≲2 keV) band than the MOS detectors. We extracted the source spectra using a 10″-radius circle. For the background spectrum, we chose a source-free region with a radius of 20″. As a result, we obtained 18 source counts in the 0.2–10 keV band from the PN instrument, while the MOS cameras provided 11 counts.

The spectra were grouped to ensure at least 1 count per energy bin and fitted simultaneously in the 0.2–10 keV band. The reddening for J1513 derived from the optical light curves analysis (Table 3) was converted to the equivalent hydrogen column density, N_H = 6 × 10²⁰ cm⁻². This value was fixed during the fitting procedure. The PL model resulted in the photon index $\mathrm{\Gamma }=2.0_{-0.4}^{+0.3}$, the unabsorbed flux in the 0.5–10 keV band $F_X=6.2_{-1.6}^{+2.2}\times {10}^{-15}$ erg s⁻¹ cm⁻², and W = 29 for 44 degrees of freedom.

5. Discussion

We have presented results of optical and X-ray observations of J2017 and J1513. For the first time, both pulsars were detected in X-rays and J1513 was firmly identified in the optical. Below, we discuss their characteristics.

5.1. J2017

J2017 belongs to the class of short orbital period systems and shows a typical optical light curve observed in spider pulsars. Modelling the optical data suggests a small companion mass of ≈0.04 M_⊙, which confirms the BW nature of J2017. The system NS mass is apparently larger than the canonical one of 1.4 M_⊙. The companion is cool, with a night-side temperature of $T_{\mathrm{n}}={3000}_{-100}^{+200}$ K. The difference between this and the temperature of a hotter day-side is about a factor of two. These values are typical for BWs (e.g. Mata Sánchez et al. 2023). The distance, $D=2.{40}_{-0.05}^{+0.10}$ kpc, is larger than the DM distance estimates (Table 1). This is not surprising, as according to Koljonen & Linares (2023), DM distances to spider pulsars are systematically lower than those derived by other methods. The E(B − V) value is consistent with the expected value of $0.{11}_{-0.02}^{+0.01}$ mag derived for the line-sight of the pulsar and its distance using the 3D Galactic dust map by Green et al. (2019).

The X-ray spectrum of the system can be described by the PL model with a photon index of 2.7 ± 0.2. The X-ray luminosity in the 0.5–10 keV band is ≈10³¹ erg s⁻¹ for the distance D = 2.4 kpc. The flux and corresponding luminosity in the 0.1–100 GeV range is F_γ^J2017 = 6.5(6)×10⁻¹² erg s⁻¹ cm⁻² (Abdollahi et al. 2022) and L_γ^J2017 = 4.5 × 10³³ erg s⁻¹. These values are in agreement with those obtained for other BWs (Swihart et al. 2022).

The irradiation efficiency, η = 0.6, calculated for the observed $\dot{E}$ overshoots values typical for BWs (e.g. Draghis et al. 2019; Mata Sánchez et al. 2023). The γ-ray efficiency, ${\eta }_{\gamma }=L_{\gamma }/\dot{E}=0.58$, is also sufficiently high. The correction of the pulsar spin-down luminosity for acceleration due to differential Galactic rotation (Nice & Taylor 1995; Lynch et al. 2018) results in a slightly lower intrinsic value, ${\dot{E}}_i=7.7\times {10}^{33}$ erg s⁻¹, implying higher efficiencies. However, as was noted above, J2017 can be a rather massive NS with a larger than canonical moment of inertia, and therefore larger ${\dot{E}}_i$. Assuming a plausible NS radius, R_p = 13 km (e.g. Salmi et al. 2022; Vinciguerra et al. 2024), and using the formula from Ravenhall & Pethick (1994), we obtained ${\dot{E}}_i=1.7\times {10}^{34}$ erg s⁻¹ for the lower bound of the NS mass, M_p = 1.8 M_⊙ (Table 3). The corresponding η = 0.28 and η_γ = 0.26 are more plausible for BW systems, and they can be even lower if J2017 is more massive and/or if the beaming factor correction is taken into account (e.g. Draghis et al. 2019; Romani et al. 2021).

5.2. J1513

J1513 is an example of a spider pulsar with an asymmetric optical light curve. This feature appears in both OAN-SPM and Magellan observations performed about three years apart, indicating that this property is persisting. The simple model of direct heating of the companion cannot reproduce the asymmetry. There are a number of other spider systems demonstrating similar behavior; for example, BW J1311−3430 (Romani et al. 2012), BW J1653−0158 (Nieder et al. 2020), BW J1810+1744 (Romani et al. 2021), RB J2039−5617 (Clark et al. 2021), and RB J2339−0533 (Kandel et al. 2020). Different models have been proposed to explain the asymmetries: cold spots on a companion star surface caused by its magnetic activity (Clark et al. 2021), asymmetric heating from the intrabinary shock (IBS; Romani & Sanchez 2016), IBS particles channelling to the magnetic poles of the companion star (Sanchez & Romani 2017), and heat redistribution over the companion surface via convection and diffusion (Kandel & Romani 2020; Voisin et al. 2020). It is not clear which of those ideas is truly applicable to describe J1513 due to the possible presence of puzzling periodic brightness wavering with a peak-to-peak amplitude of ∼0.1 mag and a periodicity of ∼15 min in the light curve in the ∼0.15–0.45 orbital phases (Fig. 3). No such variations have been reported so far for other spider systems.

Possible wavering could be attributed to the flaring activity of J1513 , which is observed in some other spider systems (e.g. Romani et al. 2012; Halpern 2022; Swihart et al. 2022; Zyuzin et al. 2024). However, the observed periodicity calls this suggestion into question as flaring is usually a stochastic process. In the case of our data, it is unclear whether these variations always appear in the same orbital phases.

On the other hand, if the variations are proved in the future, then they could be explained by the sunspot-like oscillations (e.g. Kolobov et al. 2016; Sych 2016) of the companion’s hot spot, produced by IBS particles being channelled to the off-axis magnetic pole of the companion star (see Fig. 3 in Sanchez & Romani 2017).

Though the used symmetric model does not describe the J1513 light curves for the whole phase range, it provides preliminary constraints on the system’s parameters. The night- and day-side temperatures are close to those obtained for J2017 (Table 3). The derived distance to J1513 is about 2 kpc, which is close to the value provided by the NE2001 model. Taking into account the pulsar proper motion, μ = 7.2(2) mas yr⁻¹, the transverse velocity is v_t = 68 km s⁻¹, which is consistent with the velocity distribution for binary pulsars (Hobbs et al. 2005). The best-fitting interstellar absorption for the pulsar is below the maximum value in this direction, $0.{12}_{-0.01}^{+0.02}$ mag, derived from the dust map of Green et al. (2019).

The J1513 system has a weak counterpart in X-rays. Its spectrum can be described by the PL model with a photon index of $2.0_{-0.4}^{+0.3}$. The X-ray and γ-ray luminosities of J1513 for a distance of 2 kpc are $L_{\mathrm{X}}^{\mathrm{J}1513}=3.0\times {10}^{30}$ erg s⁻¹ and L_γ^J1513 = 3.5 × 10³³ erg s⁻¹. The latter was calculated from the γ-ray flux F_γ^J1513 = 7.4(6)×10⁻¹² erg s⁻¹ cm⁻² (Abdollahi et al. 2022). These values are typical of BW systems (Swihart et al. 2022).

The intrinsic spin-down luminosity of J1513 corrected for the Shklovskii effect (Shklovskii 1970) and acceleration due to the Galactic differential rotation is ${\dot{E}}_i=1.7\times {10}^{35}$ erg s⁻¹ for D = 2 kpc, M_p = 1.7 M_⊙, and R_p = 13 km. This implies the efficiencies η_γ ≈ 0.02 and η ≈ 0.03. While the γ-ray efficiency is typical of BWs, the irradiation efficiency is very low. The lowest values provided by Mata Sánchez et al. (2023) for a sample of BWs are 0.068 ± 0.010 for PSR J2256−1024 and 0.09 ± 0.05 for PSR J0251+2606. J1513 has a rather high spin-down luminosity in comparison with other BW pulsars, including that of J2017. However, its orbital separation is also larger, which could imply a similar heating rate of J2017 and J1513 companions.

6. Conclusions

We performed first multi-band time-series optical photometry of two BW pulsars, J2017 and J1513. This allowed us to firmly identify the latter in the optical. Both objects show similar peak-to-peak amplitudes of orbital modulation, ≳2 mag. However, while the J2017 light curves are symmetric, the J1513 demonstrates strong asymmetry and periodic wavering of unknown origin.

The light curve modelling provides parameters for the J1513 and J2017 companions (masses, radii, temperatures, and Roche-lobe filling factors) comparable with values derived for other BWs (e.g. Zharikov et al. 2019; Swihart et al. 2022; Mata Sánchez et al. 2023). Inclinations of the systems are about 70°, which is compatible with observed irregular eclipses of the pulsars’ radio emission.

For the first time, J2017 and J1513 are identified in X-rays. Their X-ray spectra can be described by the absorbed PL models providing photon indices and luminosities consistent with the values derived for the BW sample.

New optical observations of J1513 with a high time resolution are necessary to confirm the wavering and to clarify its nature. Detailed multi-band observations would provide information about source colours, which could help us to understand the mechanism causing the light curves asymmetry and to better constrain the system parameters.

Acknowledgments

We thank the anonymous referee for useful comments. The work is based on observations carried out at the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile, the Gran Telescopio Canarias (GTC), installed at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma, the Observatorio Astronómico Nacional on the Sierra San Pedro Mártir (OAN-SPM), Baja California, México and on observations obtained with XMM-Newton, a ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. We thank the daytime and night support staff at the OAN-SPM for facilitating and helping obtain our observations. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. AYK acknowledges the DGAPA-PAPIIT grant IA105024. SVZ acknowledges PAPIIT grant IN119323. DAZ thanks Pirinem School of Theoretical Physics for hospitality. The analysis of the X-ray data by AVK and DAZ was supported by the Russian Science Foundation project 22-12-00048, https://rscf.ru/project/22-12-00048/. Optical data reduction produced by AVB and YAS was supported by the baseline project FFUG-2024-0002 of the Ioffe Institute. REM acknowledges support by BASAL Centro de Astrofísica y Tecnologías Afines (CATA) FB210003.

References

All Tables

All Figures

Fig. 1. Optical images of the fields of the two pulsars. Left panel: 2.​​′8 × 1.​​′4 FoV of the HiPERCAM instrument in the rs band containing J2017. Right panel: 2.​​′7 × 2.​​′7 image of the J1513 field in the r′ band obtained with the Magellan-1 telescope. The pulsars’ vicinities within the small boxes are enlarged in the insets to demonstrate the maximum (bottom) and minimum (top) brightness phases of their companions, the positions of which are marked by blue circles.

In the text

Fig. 2. Light curves of the two pulsars and results of their modelling. Top panel: gs, rs, is, and zs band light curves of the J2017 system folded with the orbital period and the best-fitting model (solid lines). Two periods are shown for clarity. Bottom panel: J1513 light curves folded with the orbital period. Data obtained in different filters and different nights are shown with various colours and symbols, indicated in the legend. All magnitudes are in the AB system. The best-fitting model is shown by the solid line. Shaded areas mark the part of the light curves where small-scale wavering was detected (see text for details). In both panels, points with dashed error bars lie below the 3σ detection limits. The phases ϕ = 0.0 correspond to the minima of the systems’ model brightness. Sub-panels show residuals calculated as the difference between the observed (O) and calculated (C) magnitudes for each data point in terms of the magnitude error, σ.

In the text

Fig. 3. Result of the small-scale brightness variations analysis. Top panel: J1513 filtered light curve in the r′ band (black points with error bars). Points within shaded areas are excluded from the periodicity search. The red line is the best fit of the data by a sine function with the period of 15.3 min found from the periodicity search. The sub-panel shows fit residuals. Middle panel: power spectrum of the filtered light curve (the solid line) with its uncertainties (vertical error bars). The largest peak corresponds to the 15.3 min period. Bottom panel: light curve folded with this period and fitted with the sine function of the same period. The sub-panel shows fit residuals.

In the text

	Fig. 4. 3.​​′5 × 3.​​′5 XMM-Newton/PN images of the J2017 (left) and J1513 (right) fields in the 0.2–2 keV band. The pulsars’ radio positions (Table 1) are marked by the ‘X’ symbols. Dashed circles show regions used for the background extraction.
In the text