#  The black widow pulsar J1641+8049 in the optical, radio and X-rays

A. Yu. Kirichenko1,2, S. V. Zharikov1, A. V. Karpova2, E. Fonseca3,4, D. A.
Zyuzin2, Yu. A. Shibanov2, E. A. López5,6, M. R. Gilfanov7,8, A. Cabrera-
Lavers9,10, S. Geier9,10, F. A. Dong11, D. C. Good12, J. W. McKee13,14, B. W.
Meyers15,11, I. H. Stairs11, M. A. McLaughlin3,4, J. K. Swiggum16
1Instituto de Astronomía, Universidad Nacional Autónoma de México, Apdo.
Postal 877, Baja California, México, 22800
2Ioffe Institute, 26 Politekhnicheskaya, St. Petersburg, 194021, Russia
3Department of Physics and Astronomy, West Virginia University, PO Box 6315,
Morgantown, WV 26506, USA
4Center for Gravitational Waves and Cosmology, West Virginia University,
Chestnut Ridge Research Building, Morgantown, WV 26505, USA
5Instituto de Física, Universidad Nacional Autónoma de México, POB 20-364,
Cd.Mx. 01000, México
6Instituto de Investigación en Ciencias Físicas y Matemáticas, USAC, Ciudad
Universitaria, 01012, Zona 12, Guatemala
7Space Research Institute, Russian Academy of Sciences, Profsoyuznaya 84/32,
117997 Moscow, Russia
8Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85741
Garching, Germany
9Instituto de Astrofísica de Canarias, Vía Láctea s/n, E38200, La Laguna,
Tenerife, Spain
10GRANTECAN, Cuesta de San José s/n, E-38712, Breña Baja, La Palma, Spain
11Department of Physics and Astronomy, University of British Columbia, 6224
Agricultural Road, Vancouver, BC V6T 1Z1 Canada
12Department of Physics and Astronomy, University of Montana, 32 Campus Drive,
Missoula, MT
13E. A. Milne Centre for Astrophysics, University of Hull, Cottingham Road,
Kingston-upon-Hull, HU6 7RX, UK
14Centre of Excellence for Data Science, Artificial Intelligence and Modelling
(DAIM), University of Hull, Cottingham Road, Kingston-upon-Hull, HU6 7RX, UK
15International Centre for Radio Astronomy Research (ICRAR), Curtin
University, Bentley, WA 6102, Australia
16Department of Physics, Lafayette College, Easton, PA 18042, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

PSR J1641+8049 is a 2 ms black widow pulsar with the 2.2 h orbital period
detected in the radio and $\gamma$-rays. We performed new phase-resolved
multi-band photometry of PSR J1641+8049 using the OSIRIS instrument at the
Gran Telescopio Canarias. The obtained data were analysed together with the
new radio-timing observations from the Canadian Hydrogen Intensity Mapping
Experiment (CHIME), the X-ray data from the Spectrum-RG/eROSITA all-sky
survey, and all available optical photometric observations. An updated timing
solution based on CHIME data is presented, which accounts for secular and
periodic modulations in pulse dispersion. The system parameters obtained
through the light curve analysis, including the distance to the source 4.6–4.8
kpc and the orbital inclination 56–59 deg, are found to be consistent with
previous studies. However, the optical flux of the source at the maximum
brightness phase faded by a factor of $\sim$2 as compared to previous
observations. Nevertheless, the face of the J1641+8049 companion remains one
of the most heated (8000–9500 K) by a pulsar among the known black widow
pulsars. We also report a new estimation on the pulsar proper motion of
$\approx$2 mas yr-1, which yields a spin down luminosity of
$\approx$4.87$\times 10^{34}$ ergs s-1 and a corresponding heating efficiency
of the companion by the pulsar of 0.3–0.7. The pulsar was not detected in
X-rays implying its X-ray-luminosity was $\la$3 $\times$ 1031 erg s-1 at the
date of observations.

###### keywords:

stars: neutron – binaries: general – pulsars: individual: PSR J1641+8049

††pubyear: 2023††pagerange:  The black widow pulsar J1641+8049 in the optical,
radio and X-rays– The black widow pulsar J1641+8049 in the optical, radio and
X-rays

##  1 Introduction

Among about 3400 pulsars discovered to date, more than 550 belong to the class
of millisecond pulsars (MSPs; Manchester et al.
(2005))111https://www.atnf.csiro.au/people/pulsar/psrcat/. These objects have
short spin periods ($P<30$ ms) and low spin-down rates ($\dot{P}\sim
10^{-20}$–$10^{-18}$ s s-1). The most generally accepted scenario implies that
MSPs are old neutron stars (NSs) which were spun-up (or ‘recycled’) through
angular momentum transfer by accretion from their main-sequence companions
during a low-mass/intermediate-mass X-ray binary stage (Bisnovatyi-Kogan &
Komberg, 1974; Alpar et al., 1982). The binary MSP population comprises
several classes depending on the type of the companion. In the so-called
‘spider’ systems with tight orbits ($P_{b}\lesssim 1$ d), low-mass companions
are heated and ablated by the pulsar wind of relativistic particles and
radiation (e.g. Manchester, 2017). Evaporated material often causes eclipses
of the pulsar radio emission. Black widows (BWs), which represent a subclass
of such binaries, have very low-mass ($M_{c}\lesssim 0.05$M⊙) almost ablated
degenerate companions. The origin and formation of such systems is not well
understood, but it is actively discussed (Chen et al., 2013; Benvenuto et al.,
2014, 2015; Ablimit, 2019; Ginzburg & Quataert, 2021; Guo et al., 2022).

About 70 BWs have been discovered so far thanks to radio and $\gamma$-ray
observations. Half of them reside in the Galactic disk (Swihart et al., 2022),
while others are associated with globular clusters222See https://www3.mpifr-
bonn.mpg.de/staff/pfreire/GCpsr.html. According to the recent census (Swihart
et al., 2022), only about 20 BWs have been detected in the optical. However,
optical studies allow one to determine fundamental parameters of BWs such as
spectral type and temperature of a companion, irradiation efficiency,
distance, as well as masses of its components when they cannot be derived from
radio timing observations alone.

Figure 1: The $3.05\times 3.05$ arcmin2 J1641 FoV imaged with the GTC/OSIRIS in the $r^{\prime}$ band. The J1641 vicinity is shown by the box in the centre of the image, and the optical companion is indicated by the arrow. The left and right inserts correspond to the enlarged J1641 vicinity imaged near the companion maximum and minimum brightness phases. The stars from the Pan-STARRS catalogue listed in Table 2 and used for the photometric calibration, are marked by the capital letters.  Table 1: Log of the J1641 observations with GTC/OSIRIS. Date | Filter | Exposure time, | Airmass | Seeing,
---|---|---|---|---
|  | s |  | arcsec
12/07/2020 | $r^{\prime}$ | 150$\times$43 | 1.62$-$1.70 | 0.8$-$1.0
13/08/2020 | $g^{\prime}$ | 200$\times$21 | 1.64$-$1.76 | 0.7$-$1.1
13/08/2020 | $i^{\prime}$ | 120$\times$21 | 1.64$-$1.76 | 0.7$-$0.9
Table 2: Stars from the Pan-STARRS Catalogue detected in the J1641 field and their magnitudes. Star | $g^{\prime}$ | $r^{\prime}$ | $i^{\prime}$
---|---|---|---
A | 21.619(068) | 20.575(025) | 19.598(017)
B | - | 21.295(137) | 20.717(019)
C | - | 21.316(044) | 20.636(032)
D | 21.142(054) | 19.970(021) | 18.980(011)
E | 19.617(014) | 18.559(009) | 18.054(008)
F | 20.416(030) | 19.253(013) | 18.683(008)
G | 19.297(012) | 18.744(006) | 18.505(008)
H | 21.403(095) | 21.058(054) | 19.889(012)
Table 3: The J1641 parameters derived from observations with the CHIME/Pulsar
backend.

Logistics of CHIME/Pulsar radio observations
---
MJD Range of Observations, MJD | 59030–60063
Frequency range, MHz | 400–800
Number of channelized TOAs per epoch | 32 (maximum)
Number of observing epochs | 581
Number of TOAs, total | 7793
tempo goodness-of-fit ($\chi^{2}$) statistic | 9970
tempo reduced $\chi^{2}$ statistic | 1.04
tempo root-mean-square residual, $\mu$s | 1.88
tempo scaling factor for raw uncertainties | 1.7
Best-fit timing parameters from tempo
Right ascension (J2000), hh:mm:ss | 16:41:20.83311(3)
Declination (J2000), dd:mm:ss | 80:49:52.92335(8)
Proper motion in R. A., mas yr-1 |  $-$1.74(10)
Proper motion in Dec., mas yr-1 |  $-$1.03(11)
Timing parallax, mas | $-4(4)$
Spin frequency $f$, Hz | 494.760636473171(4)
Time derivative in spin frequency, Hz2 | $-2.3928(2)\times 10^{-19}$
Epoch of reference for pulsar spin, MJD | 59315.1255
DM used for observations, pc cm-3 | 31.091883
Best-fit DMs, binned | see Section 3
Binary model used by tempo | ELL1
Projected semi-major axis, lt-s | 0.06407531(17)
Epoch of passage through asc.-node longitude, MJD | 59315.04093360(5)
ELL1 eccentricity parameter #1 | 0.000019(4)
ELL1 eccentricity parameter #2 | $-0.000010(4)$
Orbital frequency $n_{b}$, Hz | $1.2736404201(2)\times 10^{-4}$
First time-derivative in $n_{b}$, Hz2 | $-9.9(3)\times 10^{-20}$
Second time-derivative in $n_{b}$, Hz3 | $-1.89(16)\times 10^{-27}$
Derived parameters from CHIME/Pulsar radio observations
Proper motion $\mu$, mas yr-1 | $2.02(10)$
Observed spin-period derivative $\dot{P}$, s s-1 | $9.7748(9)\times 10^{-21}$
Intrinsic spin-period derivative $\dot{P}_{i}$, s s-1 | $10.169(12)\times 10^{-21}$
Mass function $f_{M}$, M⊙ | $3.42\times 10^{-5}$
Characteristic age $\tau_{c}\equiv P/2\dot{P}$, Gyr | 3.28
Observed spin-down luminosity $\dot{E}$, erg s-1 | $4.67\times 10^{34}$
Intrinsic spin-down luminosity $\dot{E}_{i}$, erg s-1 | $4.86\times 10^{34}$
Magnetic field at the equator $B_{\rm eq}$, G | $1.42\times 10^{8}$

  * •

Numbers in parentheses denote 1$\sigma$ uncertainties relating to the last
significant digit quoted.

  * •

Spin-down luminosities are calculated assuming the canonical moment of inertia
of 1045 g cm2. $\dot{E}_{i}$ is calculated based on the distance $D=4.6(2)$
kpc derived from our modelling (see text).

The binary MSP PSR J1641+8049 (hereafter J1641) was discovered in the Green
Bank North Celestial Cap (GBNCC) pulsar survey (Stovall et al., 2014; Lynch et
al., 2018). This is an eclipsing radio pulsar which is also presented in the
list of Fermi-detected pulsars (Ray, 2023). Its flux in the 0.1–100 GeV range
is $2.0(3)\times 10^{-12}$ erg s-1 cm-2 (Abdollahi et al., 2022). J1641 was
not observed in X-rays. It has an orbital period of 2.18 h which is one of the
shortest among the known BWs. The companion’s minimum mass was estimated to be
0.04 M⊙ (Lynch et al., 2018). The dispersion measure (DM) distances to the
pulsar are $D_{\rm YMW16}$ = 3.0 kpc and $D_{\rm NE2001}$ = 1.7 kpc based on
the YMW16 (Yao et al., 2017) and NE2001 (Cordes & Lazio, 2002) models for the
distribution of free electrons in the Galaxy, respectively.

Lynch et al. (2018) found a faint optical counterpart ($r=24.0(3)$) to J1641
in the archival data, and performed photometric observations of the pulsar
with the 4.3-m Lowell Discovery Telescope333Formerly known as the Discovery
Channel Telescope (LDT) in the $g$, $r$, $i$ and $z$ filters and with the
McDonald Observatory 1 m telescope in the $r^{\prime}$ and $i^{\prime}$
filters. The counterpart revealed strong brightness variations tied to the
orbital period, confirming the optical identification of the pulsar companion.
Further optical studies of J1641 were recently reported by Mata Sánchez et al.
(2023) who performed its phase-resolved multiband photometry with the HiPERCAM
instrument at the 10.4-m Gran Telescopio Canarias (GTC) in 2019. They analysed
the multiband light curves and obtained the fundamental parameters of the
system, including the inclination, the Roche-lobe filling factor, the
companion mass and the temperature gradient over the companion surface, as
well as the distance to the system.

In this paper, we report the results of our independent phase-resolved multi-
band optical observations of J1641 obtained with the OSIRIS instrument at the
GTC in 2020. We analyse the OSIRIS, HiPERCAM, LDT, and the 1 m telescope data
of J1641 all together. In addition, we present the updated parameters of the
system derived from ongoing observations with the Canadian Hydrogen Intensity
Mapping Experiment (CHIME) telescope, and report an upper limit on the pulsar
X-ray flux based on observations with eROSITA (Predehl et al., 2021) aboard
the Spectrum-RG (SRG) orbital observatory (Sunyaev et al., 2021). The paper is
organised as follows: observations and data reduction are described in Sec. 2,
the radio timing analysis is presented in Sec. 3, while the modelling of the
light curves is described in Sec. 4. Discussion and conclusions are given in
Sec. 5.

##  2 Observations and data reduction

###  2.1 Optical data

The phase-resolved photometric observations444Proposal GTC11-20AMEX, PI A.
Kirichenko of the J1641 field were carried out during two observing runs in
the Sloan $g^{\prime}$, $r^{\prime}$ and $i^{\prime}$ bands with the Optical
System for Imaging and low Resolution Integrated Spectroscopy (OSIRIS)
instrument at the GTC. In order to reduce the CCD readout time and increase
the efficiency of the phase-resolved observations, we used windowing. The
target was exposed on CCD1, and the windowed FoV was $3.05\times 3.05$
arcmin2. To avoid effects from CCD defects, 5 arcsec dithering between the
individual exposures was used in both observing runs. The observations roughly
covered two orbital periods in total, i.e. one orbital period per each
observing run. The first period was observed in the $r^{\prime}$ band only,
whereas the second one was covered one month later using the alternating
$g^{\prime}$ and $i^{\prime}$ bands. The log of observations is given in Table
1. The $r^{\prime}$-band image of the pulsar field is presented in Fig. 1,
where the inserts demonstrate the variability of the pulsar companion.

Using the Image Reduction and Analysis Facility (iraf) package, we performed a
standard data reduction, including bias subtraction and flat-fielding. The
cosmic rays were removed from all images with the L.A.Cosmic algorithm (van
Dokkum, 2001). Astrometric referencing was performed using a single 180-s
$r^{\prime}$-band image and a set of stars from the Gaia DR3 Catalogue (Gaia
Collaboration et al., 2023). Given the reduced image size in the windowing
mode, only five field stars detected by Gaia appeared to be suitable for the
astrometric purposes. Using these stars, we computed the astrometric solution
with the formal rms uncertainties $\Delta$RA $\la$ $0.13$ arcsec and
$\Delta$Dec $\la$ $0.14$ arcsec.

Photometric calibration was performed using the Sloan photometric standards
SA111-1925 for the $r^{\prime}$ band and PG1528+062B for the $g^{\prime}$ and
$i^{\prime}$ bands (Smith et al., 2002) observed during the same nights as the
target. Using their instrumental magnitudes and the site extinction
coefficients $k_{g^{\prime}}$ = 0.15(2), $k_{r^{\prime}}$ = 0.07(1), and
$k_{i^{\prime}}$ = 0.04(1) (Cabrera-Lavers et al., 2014), we calculated the
zero points $Z_{g^{\prime}}=28.53(2)$, $Z_{r^{\prime}}=28.94(1)$, and
$Z_{i^{\prime}}=28.42(1)$. To verify these zero points, we checked the sky
transparency variability by comparing the instrumental magnitudes of a field
star during the observations of the target. The variation did not exceed
errors of the flux measurements for the used star. Nevertheless, since the
standards and the target were observed at different sky positions, we also
compared the magnitudes of several field stars with those from the Sloan
Digital Sky Survey (SDSS) Release 14 Catalogue (Abolfathi et al., 2018) and
the Pan-STARRS Catalogue (Flewelling et al., 2020). To convert the Pan-STARRS
magnitudes to the SDSS photometric system, we used the equation (6) from Tonry
et al. (2012). The stars used in this analysis are shown in Fig. 1 and listed
in Table 2. Their catalogue $g^{\prime}$-band magnitudes were found to be
consistent within uncertainties with those calibrated using the photometric
standard. However, in case of the $r^{\prime}$ and $i^{\prime}$ bands, we
found a slight discrepancy between the respective magnitudes. In addition, we
calculated the colour term corrections and found that they are negligible in
the $r^{\prime}$ and $i^{\prime}$ bands, and only slightly affect the
$g^{\prime}$-band measurements. Accounting for all of the mentioned
corrections, the resulting zero points are $Z_{g^{\prime}}=28.51(2)$,
$Z_{r^{\prime}}=29.03(3)$, and $Z_{i^{\prime}}=28.34(3)$. The point source
$3\sigma$ upper limits in individual $g^{\prime}r^{\prime}i^{\prime}$
exposures are $g^{\prime}=26.4$, $r^{\prime}=25.8$, and $i^{\prime}=24.8$. We
note that all magnitudes presented in this paper are in the AB system.

###  2.2 Radio data

J1641 is being observed by the CHIME telescope (CHIME; CHIME Collaboration et
al., 2022). For our present work, we processed and analysed high-cadence
timing data recorded with the pulsar-timing backend built for CHIME
(CHIME/Pulsar Collaboration et al., 2021). The CHIME/Pulsar backend generates
folded profiles evaluated over 10-s integrations and 1,024 frequency channels
that span the 400–800 MHz range. These data are coherently dedispersed in real
time and prior to folding, using the value of DM listed in Table 3.

The CHIME/Pulsar data set on J1641 currently spans $\sim$3 years of
observations between early 2020 and 2023, and is being collected in support of
ongoing GBNCC analyses (McEwen & et al., prep). All data were processed –
through statistical cleaning of radio-frequency interference and downsampling
to 32 channels – using the psrchive (van Straten et al., 2012) and clfd
(Morello et al., 2019) analysis suites. Once initially cleaned, a subset of
data was co-added to form a high-significance pulse profile that turned into a
‘standard template’ after de-noising; this standard template was then used to
compute times of arrival (TOAs) for the entire CHIME/Pulsar data set, yielding
32 TOAs per epoch that are each evaluated across 32 downsampled frequency
channels. Any TOAs with S/N values less than 8.0, as determined with the pat
utility in psrchive, were excluded from analysis based on sub-optimal
detection statistics. For the remaining data set, a small, additional amount
of TOAs ($\ll$1 per cent) were excised during the timing analysis due to
corruption of the pulse profile from sub-threshold interference that was not
detected during the data-preparation process.

###  2.3 X-ray data

The J1641 field was observed in the course of SRG/eROSITA all-sky survey in
four visits spanning between April 9, 2020 and Oct. 13, 2021 with the total
exposure time of 2.8 ks. No source was statistically significantly detected at
the pulsar position. The derived upper limit on the unabsorbed flux is
$\approx 1.3\times 10^{-14}$ erg s-1 cm-2 in the 0.5–10 keV range (90 per cent
confidence), assuming a power law (PL) model with the photon index $\Gamma$ =
2.5 (the average value for BWs, Swihart et al., 2022) and the absorbing column
density $N_{\rm H}=7\times 10^{20}$ cm-2. The latter was derived using the
reddening $E(B-V)=0.08$ mag obtained for J1641 from the extinction map of
Green et al. (2019) and the empirical relation from Foight et al. (2016). Note
that the reddening is equal to that obtained from the optical light curve
modelling (see below).

##  3 Updated Radio Timing of PSR J1641+8049

We refined the timing model of J1641, based on the solution developed by Lynch
et al. (2018), using the tempo pulsar-timing package to obtain updated
estimates of the spin, astrometric, and orbital parameters based on the
CHIME/Pulsar data set alone. We also incorporated 355 additional degrees of
freedom to fit for DM in contiguous time bins of 1.5-day extent that span the
timespan of the data set. Due to the turbulent environments often observed in
BW systems, we also explored the fitting of parameters that quantify
variations in the orbital elements.

Figure 2: Best-fit timing residuals ($\mathcal{R}$, top) and DM values
determined for the CHIME/Pulsar data set described in Section 2.2. The DM
values are plotted as changes relative to the value reported in Table 3.

A summary of the best-fitting CHIME/Pulsar timing residuals for J1641 and DM
timeseries is shown in Fig. 2, and timing model parameters are reported in
Table 3. One important result from our updated modelling is that the
CHIME/Pulsar data yielded statistically different estimates of the proper
motion than those obtained by Lynch et al. (2018), with the new magnitude of
proper motion $\mu=2.02(10)$ mas yr-1 being lower by a factor of $\sim$20\. We
attribute this difference to the high-cadence nature of the CHIME/Pulsar TOAs,
which allow for better estimates of short-term variations that are typically
observed in other BW systems. This change in proper motion mainly impacts the
derived estimate of ‘intrinsic’ spin-down of the pulsar, i.e., the time rate
of change in spin frequency corrected for biases induced by proper motion and
acceleration in the Galactic potential (e.g., Nice & Taylor, 1995). The
corrected spin-down we derive for J1641 is discussed in Sec. 5 as it depends
on results obtained from the optical analysis presented below.

Our timing model of orbital motion uses the ELL1 formalism to describe low-
eccentricity orbits (Lange et al., 2001). The best-fit ELL1 model indicates
that deviations from purely-Keplerian motion are detectable in the
CHIME/Pulsar data set. The inclusion of one time-derivative in orbital
frequency $n_{b}=2\pi/P_{b}$ as a degree of freedom improves the fit of the
timing model by an amount $\Delta\chi^{2}\approx 17,000$. We found that
fitting for two time-derivatives in $n_{b}$, along with the five Keplerian
elements, yielded an optimal fit to the timing data. The interpretation of
these parameters in terms of macroscopic quantities (e.g., mass and/or
geometry) is nontrivial due to the complex environment in BW systems that
produce stochastic orbital variations (e.g., Shaifullah et al., 2016). Future
analysis of orbital evolution in the J1641 system will be performed once
several years of additional CHIME/Pulsar data are obtained.

The reliable detection of J1641 indicates minimal eclipsing of radio signal at
superior conjunction. However, the DM variations in Fig. 2 show secular and
quasi-periodic trends over time. The estimate of $\mu$ for J1641 nonetheless
remains comparable whether we use a many-bin DM model or a polynomial
expansion in DM, which differ in functional form and by hundreds of fit
parameters. This circumstance indicates that proper motion is robustly
measured with the CHIME/Pulsar data set despite only spanning $\sim$3 years in
time and regardless of the choice in DM model.

In order to better assess these variations, we generated a separate timing
solution based on data acquired in the MJD 59000–59400 range555While
arbitrary, this range was chosen since the default version of tempo is unable
to accommodate the total number of per-epoch DM parameters needed to model the
entire CHIME data set. and setting maximum DM-bin extents to be 0.5 days,
i.e., to ensure a single DM is estimated for each observing epoch. Despite the
shortened data set, we fit for the same parameters reported in Table 3 in
order to separately measure variations in celestial position, DM and orbital
motion.

Figure 3: Per-epoch measurements of DM for J1641 plotted as a function of time
(top) and orbital phase (bottom), estimated by fitting to all CHIME/Pulsar
data within the MJD 59000–59400 range. In the top panel, the grey line is a
best-fit estimate of the secular variation over time that is presumed to be
linear. In the bottom panel, the DM values have been corrected for the secular
variation.

The results of this data-subset modelling is shown in Fig. 3, which shows this
residual subset as function of time and orbital phase. In the ELL1 binary
model, the orbital phase $\Phi=n_{b}(t-T_{\rm asc})/(2\pi)$, where $T_{\rm
asc}$ is the epoch of passage through the longitude of ascending node
specified in Table 3; superior conjunction corresponds to $\Phi=0.25$ in
Figure 3. The fitting of per-epoch DMs allows for better resolution of
periodic variations, though requires at least 270 DM-bin fit parameters for
the MJD 59000–59400 portion of the data alone. Nonetheless, the per-epoch DMs
exhibit clear periodic variations that occur on a timescale equal to the
orbital period for J1641; the $\sim$monthly variation in the DM timeseries is
a manifestation of aliasing due to CHIME/Pulsar observations occuring once
every sidereal day. While an excess in DM coincides with superior conjunction,
the DM appears to modulate over the whole orbit and thus indicates a
structured circumbinary medium. Further analysis of the DM variations over the
entirety of the J1641 data set, as well as other BW systems discovered and
monitored by GBNCC, will be presented in future work.

##  4 PSR J1641+8049 optical light curves and the system parameters

Figure 4:  Top-left: Multi-colour light curves of J1641 obtained with the
OSIRIS/GTC ($g^{\prime}$, $r^{\prime}$, $i^{\prime}$) folded with the orbital
period. The phase zero corresponds to the orbital phase when the companion
day-side is facing an observer. Two orbital cycles are shown for clarity.
Solid lines represent the result of the best fit of the data by the model
where the companion is heated by the pulsar. Horizontal dashed lines show
3$\sigma$ detection limits of the observations. Photometric bands are marked
by different colours. Top-right: The HiPERCAM/GTC multi-band data of J1641 and
the best fit of the $g^{\prime}_{s}$, $r^{\prime}_{s}$, $i^{\prime}_{s}$ light
curves by the model. Bottom panels: Fit residuals calculated as the difference
between the observed ($O$) and the calculated ($C$) magnitudes for each data
point in terms of the magnitude error $\sigma$.

Figure 5: The broad-band optical spectra of the J1641 optical component at the
photometric maximum for different epochs of observations.

Figure 6: The r′-band light curves folded with the orbital period obtained
with the GTC, LDT, and McD. The orbital phase $\phi=0.0$ corresponds to the
maximum of the GTC/HiPERCAM light curve and it is shifted to -0.76 from the
radio convention orbital phase, where $\phi=0.0$ corresponds to the ascending
node. LDT and McDonald 1-m telescope data were obtained in March 2017.

The resulting $g^{\prime}$, $r^{\prime}$, $i^{\prime}$-band light curves
folded with the orbital period are presented in Fig. 4, left. In the right
panel of Fig. 4 we show the HiPERCAM data obtained in the $u^{\prime}_{s}$-,
$g^{\prime}_{s}$-, $r^{\prime}_{s}$-, $i^{\prime}_{s}$-, and
$z^{\prime}_{s}$\- bands about one year before our observations (Mata Sánchez
et al., 2023). The shapes of the light curves are found to be consistent,
nevertheless the object appears to be slightly brighter and bluer at the
maximum of the light curves in the HiPERCAM data. To demonstrate this, in Fig.
5 we present the source broad-band optical spectra in the maximum of the light
curves, corresponding to the OSIRIS (blue) and the HiPERCAM (green) data. The
slopes666Without the $u^{\prime}_{s}$ band for the HiPERCAM data. of the
spectra $F\sim\lambda^{\alpha}$ are $\alpha_{\rm OSIRIS}=-1.74(21)$ and
$\alpha_{\rm HiPERCAM}=-2.01(18)$, respectively. As it can be seen, there is a
brightness-decreasing tendency together with a relative reddening of the
object spectrum on the time scale of these observations. We also note that the
flux measurements close to the photometric maximum obtained with the McDonald
Observatory 1 m telescope were brighter than those reported for the
GTC/HiPERCAM observations, whereas the LDT measurements were mostly close to
them (Lynch et al. (2018) and D. Kaplan, private communication). For
comparison, in Fig.6 we show the r′-band light curves with all available data
points.

Figure 7:  The Roche lobe of the system components and a magnification of the secondary. The colours mark the effective temperature distribution on the secondary surface. The image corresponds to the OSIRIS data modelling, when all parameters are free.  Table 4: The light-curve fitting results for J1641. Mass function $\times 10^{-5}$, M⊙ (fixed) | - | - | 3.42 | 3.42
---|---|---|---|---
Fitted parameters | OSIRIS | HiPERCAM | OSIRIS | HiPERCAM
Pulsar mass $M_{\rm NS}$, M⊙ | 2.0(6) | 1.3(4) | 1.3(3) | 1.3(1)
Mass ratio $q$ = $M_{\rm c}/M_{\rm NS}$ | 0.018(4) | 0.019(2) | 0.035 | 0.036
Distance $D$, kpc | 4.60(20) | 4.70(20) | 4.64(12) | 4.83(13)
Reddening $E(B-V)$, mag | 0.072(23) | 0.056(30) | 0.09(2) | 0.07(3)
‘Night-side’ temperature $T_{\rm n}$, K | 3300(100) | 3400${}_{-300}^{+400}$ | 3380(120) | 3500${}_{-180}^{+300}$
Inclination $i$, deg | 57(7) | 56(10) | 59(8) | 58(10)
Roche lobe filling factor $f_{x}$ | 0.95${}_{-0.10}^{+0.05}$ | 0.90${}_{-0.09}^{+0.10}$ | 0.99${}_{-0.08}^{+0.01}$ | 0.87${}_{-0.07}^{+0.13}$
Irradiation factor $K_{\rm irr}$, |  |  |  |
$\times 10^{20}$ erg cm-2 s-1 sr-1 | $4.3(3)$ | 5.8(5) | 1.9(1) | 3.6(2)
$\chi^{2}$/d.o.f. | 138/73 | 767/956 | 167/74 | 2525/957
Derived parameters
Companion mass $M_{\rm c}$, M⊙ | 0.038 | 0.025 | 0.046 | 0.047
Companion radius $R_{\rm c,x}$, R⊙ | 0.194 | 0.165 | 0.220 | 0.198
Companion radius $R_{\rm c,y}$, R⊙ | 0.143 | 0.121 | 0.167 | 0.144
Lowest ‘day-side’ temperature $T_{\rm d}^{\rm min}$, K | 3800 | 4250 | 3940 | 4410
Highest ‘day-side’ temperature $T_{\rm d}^{\rm max}$, K | 8200 | 9500 | 8000 | 8700
Irradiation efficiency $\eta$ | 0.45 | 0.7 | 0.3 | 0.5

Figure 8: The masses of the pulsar and its companion from the fits of the
optical light curves. The lines show the relation between masses at different
orbit inclinations based on the mass function found from the pulsar radio
timing. The black triangles correspond to the fits when the masses and
inclinations are free parameters. The red triangle marks the fits which take
into account the mass function. The green triangle shows the result from Mata
Sánchez et al. (2023).

To estimate the system parameters, we fitted the light curves using the
emission model of a binary system described in Zharikov et al. (2013);
Zharikov et al. (2019). The model consists of an NS as the primary which heats
a low-mass companion as the secondary. The spectrum of each surface element of
the companion is approximated by a blackbody with an effective temperature
which is distributed non-uniformly over the star surface accounting for its
heating by the pulsar. The contribution of the pulsar into the observed
optical flux is negligible for any expected distance and NS brightness values.
Following Zharikov et al. (2019), the effective irradiation factor of the
secondary is related to the heating efficiency $\eta$ and the spin-down
luminosity of the pulsar as

$K_{\rm irr}=\frac{\eta\dot{E_{i}}}{4\pi^{2}R_{\rm NS}^{2}},$ (1)

and it defines the effective radiate flux $F_{\rm in}$ transferred from the
pulsar to the secondary:

$F_{\rm in}=\mathrm{cos}(\alpha_{\rm norm})\Omega\Delta SK_{\rm irr},$ (2)

where $\alpha_{\rm norm}$ is the angle between the incoming flux and the
normal to the surface, $\Omega=\pi R^{2}_{\rm NS}/a^{2}$ is the solid angle
from which the pulsar is visible from the surface element $\Delta S$ of the
companion, $R_{\rm NS}=13$ km is the NS radius and $a$ is the orbit
separation. The corresponding ‘day-side’ temperature of the companion star
surface element is

$T_{\rm d}=T_{\rm n}\left[1+\frac{F_{\rm in}}{\sigma\Delta S(T_{\rm
n})^{4}}\right]^{1/4},$ (3)

where $\sigma$ is the Stefan-Boltzmann constant. The phase-resolved light
curves were calculated by integrating the flux from all visible elements of
the secondary in the corresponding band. The gradient descent method was used
to find the minimum of $\chi^{2}$ defined as

$\chi^{2}=\sum_{j}^{g^{\prime},r^{\prime},i^{\prime}}\sum_{k}^{N_{j}}\left(\frac{O_{k}-C_{k}}{\sigma_{k}}\right)^{2},$
(4)

where $N_{j}$ is the number of observations in a given filter, $O_{k}$,
$C_{k}$, and $\sigma_{k}$ are the observed and the calculated magnitudes, and
the error of the observed magnitude, respectively. The free fitted parameters
were the distance $D$, the reddening $E(B-V)$, the binary system inclination
$i$, the Roche lobe filling factor $f_{x}$ defined as a ratio of distances
from the centre of mass of the secondary to the star surface and to the
Lagrange point $L_{1}$, the ‘night-side’ temperature $T_{\rm n}$ of the
secondary, the effective irradiation factor $K_{\rm irr}$ [ergs s-1 cm-2
sr-1], the pulsar mass $M_{\rm NS}$, and the component mass ratio.

The best-fitting model light curves for the OSIRIS data are shown in the left
panel of Fig. 4 by solid lines. The model parameters are given in the second
column of Table 4. The uncertainty of each fitted parameter was calculated
following the method proposed by Lampton et al. (1976). The geometry of the
system and the distribution of the effective temperature at the companion
surface for the OSIRIS data are shown in Fig 7. For comparison, we also used
the model for the HiPERCAM $g^{\prime}_{s}$-, $r^{\prime}_{s}$-,
$i^{\prime}_{s}$\- band light curves. We limited our analysis of the HiPERCAM
data to the three close optical bands because a simple blackbody spectrum
approximation for radiation from a star surface element used in the model
cannot describe the companion spectrum in a wide spectral range from the
$u^{\prime}_{s}$\- up to the $z^{\prime}_{s}$-band. On the other hand, this
approach was applied to simplify the comparison of the fit results achieved
using the same model setup in both cases. Since the source was slightly
brighter and bluer in the HiPERCAM data compared to the OSIRIS data, we first
fixed all model parameters excluding the secondary heating at the values given
in Table 4 for the OSIRIS data, and then fitted the HiPERCAM data. The fit
provided a hotter day-side part of the secondary with a maximum temperature of
about 9200 K. After that, we thawed all parameters as in the case of the
OSIRIS data. The best-fitting parameters of the last fit are given in the
third column of Table 4. In general, they are close to or inside the 1$\sigma$
error range of the OSIRIS data fit results, except for the effective
irradiation factor $K_{\rm irr}$, which, in turn, gives a higher day-side
temperature of the secondary. Another difference is related to the mass of the
pulsar. The fitting of the two data sets results in a significantly lower mass
in the HiPERCAM data compared to that in the OSIRIS data. Nevertheless, their
1$\sigma$ uncertainties overlap. We note that the light curves shapes and
fluxes are mainly defined by the size of the companion Roche lobe, its filling
factor, temperature distribution and system inclination. The Roche lobe size
is changed by 15 per cent at the variation of the pulsar mass within its
reasonable limits. To get the observed light curves, these changes can be
compensated by variations of the temperature and filling factor. Thus, these
parameters are correlated and additional information is needed to better
constrain them.

For instance, we can utilise the mass function from the new radio timing
measurements (see Table 3). It decreases the numbers of free parameters,
because it links the pulsar mass, the inclination and the mass of the
companion. Taking into account the mass function we repeated the fitting of
the data from both instruments. The results are presented in the last two
columns of Table 4. They are very close to those obtained by Mata Sánchez et
al. (2023) who used the icarus code (Breton et al., 2012) to model the
HiPERCAM data (see table 3 in the respective paper). However, we note that the
formal $\chi^{2}$ values are higher as compared to the cases when the masses
of components and inclination are free parameters. All results are summarised
in Fig. 8 and discussed below.

##  5 Discussion and conclusions

As it was mentioned before, the J1641 orbital period of 2.18 h is one of the
shortest among the known BWs (Swihart et al., 2022), and its highly modulated
optical light curves are typical for such tight binary systems with MSPs. The
shapes of the light curves do not demonstrate significant changes on a one-
year time scale between the HiPERCAM and OSIRIS observations. However, the
fluxes at the photometric maximum exhibit a significant decrease between the
observations by up to a factor of two. This can indicate variability of the
pulsar wind on a short (days) and long (years) time scales. We note that
variable heating (increases and decreases of the companion irradiation by the
pulsar) is observed for some other MSP binary systems, such as, e.g., the BW
candidate 4FGL J0935.3+0901 (Halpern, 2022) or the redback PSRs J1048+2339
(Yap et al., 2019) and J2129-0429 (Bellm et al., 2016). Follow-up observations
are needed to confirm and study this effect for J1641.

The previously reported proper motion $\mu$ = $39(3)$ mas yr-1 for J1641 led
to a negative value of the pulsar intrinsic spin-down luminosity, implying a
spin-up scenario (Lynch et al., 2018; Mata Sánchez et al., 2023), which could
indicate the presence of accretion. However, the latter would complicate
detection of radio pulsations from the pulsar and affect the optical light
curves of the system. Neither the HiPERCAM nor the OSIRIS observations support
this scenario. In addition, using the distance obtained through the HiPERCAM
light-curve modelling, Mata Sánchez et al. (2023) concluded that the maximum
proper motion for this system, which would allow to avoid the spin-up
scenario, is $\mu\leq 19$ mas yr-1. Indeed, the updated proper motion $\mu$ =
$2.02(10)$ mas yr-1 derived from the new radio observations (see Sec. 2.2) is
significantly lower than the one provided by Lynch et al. (2018). It yields a
new estimation on the intrinsic period derivative, $\dot{P_{i}}$ and the
intrinsic spin-down luminosity $\dot{E_{i}}$. Considering the distance D =
4.6(2) kpc obtained from our optical light curve modelling, we estimated the
Shklovskii correction (Shklovskii, 1970), $\dot{P}_{\rm Shk}\approx 0.08\times
10^{-21}$ s s-1. Taking into account the corrections due to the differential
Galactic rotation and the pulsar acceleration, the corresponding intrinsic
period derivative is $\dot{P_{i}}=10.2\times 10^{-21}$ s s-1, and the spin-
down luminosity is $\dot{E_{i}}=4.87\times 10^{34}$ erg s-1, rejecting the
spin-up and accretion scenario. We note that in the case of this particular BW
system, the Shklovskii correction is subdominant leading to the fact that the
intrinsic spin-down luminosity is $\approx$ 4 per cent higher than the
observed one. In addition, taking into account the newly determined proper
motion, we can now estimate the J1641 transverse velocity. Considering the
distance 4.6 kpc derived from the light-curve modelling and following Verbunt
et al. (2017) to estimate the contribution of the Galactic rotation and solar
peculiar velocity, we obtain $V_{T}\approx 36$ km s-1 in the local standard of
rest of the pulsar. This velocity is typical for pulsar binary systems (Hobbs
et al., 2005).

According to Mata Sánchez et al. (2023), J1641 has one of the heaviest
companions ($M_{\rm c}$=0.055${}^{+0.016}_{-0.014}$ M⊙) among the known BWs.
More massive secondaries were found for, e.g., PSR J1555$-$2908
($0.060^{+0.005}_{-0.003}$ M⊙; Kennedy et al. 2022) and PSR J1810+1744
(0.065(1) M⊙; Romani et al. 2021). However, our fit for all free parameters
resulted in a lower companion mass of 0.035(13) M⊙ and 0.025(4) M⊙ for the
OSIRIS and HiPERCAM data, respectively. These values are close to the average
mass $\bar{m}_{c}=0.026(15)$ of pulsar companions in BW systems (Swihart et
al., 2022, see table 5 therein). However, the value for the HiPERCAM fit is
inconsistent with the mass function derived from the radio timing (see Fig.
8), while the value for the OSIRIS fit agrees with it. Taking the mass
function into account we got the mass of the companion of 0.046(10) M⊙ for
both data sets. This is in agreement with the values reported by Mata Sánchez
et al. (2023). Thus, the companion indeed can be rather heavy in comparison
with other BW systems.

Another distinct feature of J1641 is the high day-side temperature of the
companion, 8200–9500 K, making it one of the five most heated sources among
the known BWs (see Fig. 9). The other four sources are PSR J1311$-$3430 with
$T_{\rm d}\gtrsim 12000$ K (Romani et al., 2012, 2015), the BW candidate ZTF
J1406+1222 with $T_{\rm d}\approx 10500$ K (Burdge et al., 2022), PSR
J1810+1744 with $T_{\rm d}\approx 9400$ K (Romani et al., 2021) and PSR
J1555$-$2908 with $T_{\rm d}\approx 9400$ K (Kennedy et al., 2022). The J1641
day-side temperature is $\gtrsim 2$ times larger than the respective
temperatures of, e.g., PSR J0251+2606 ($\approx$3400 K) or PSR J0636+5129
($\approx$4600 K) (Mata Sánchez et al., 2023) near the low end of the source
temperature distribution. It seems natural that the number of companions with
a lower day-side temperature is larger than that with a hotter one, as the
latter have to be evaporated faster. However, it remains unclear how the
companion heating is related to the ‘spin-down flux’ defined as
$\dot{E}P_{b}^{-4/3}$ (Zharikov et al., 2019, see figure 5 and table 5
therein). It is possible that the absence of a clear dependence on the spin-
down luminosity and/or system separation indicates the importance of the
pulsar wind instability and/or asymmetry in the companion heating.

Using the intrinsic spin-down luminosity $\dot{E_{i}}$ and the irradiation
factor $K_{\rm irr}$, we estimated the irradiation efficiency range
$\eta\approx$ 0.3–0.7. These values seem to overshoot those typically observed
for BWs (e.g., Draghis et al., 2019). An irradiation luminosity even larger
than $\dot{E}$ ($\eta\ga 1$), calculated assuming a canonical momentum of
inertia value of 1045 g cm2, was derived, e.g., for BW PSR J1810+1744 by
Romani et al. (2021). It is also an additional argument that the true one can
be lower due to the beaming factor (see Draghis et al., 2019; Romani et al.,
2021).

The derived companion mass of J1641 (Table 4) is close to typical masses of
brown dwarfs (0.01–0.07 M⊙) implying the possible origin of the companion.
However, a field brown dwarf with an age of $\ga 1$ Gyr, similar to the J1641
characteristic age (Table 3), would be twice more compact ($\la$0.1 R⊙) and
colder ($\la 1500$ K, Marley et al. 2021) as compared to the derived night-
side temperature of J1641. This disfavours the brown dwarf nature of the
companion. Nevertheless, such a discrepancy appears to be not a unique
property of J1641. As was noted before, most BW companions, especially the
strongly heated, are bloated up by a factor of two in comparison with the
Galactic field brown dwarfs (Kandel & Romani, 2023, see table 6 therein). The
strong irradiation by the pulsar can affect the global structure of the
presumed brown dwarf leading to the increase of its night side temperature as
well.

Finally, J1641 is not detected in X-rays in the 2.8 ks data accumulated during
the SRG/eROSITA all-sky survey. The upper limit (90 per cent) on its X-ray
luminosity is about 3.3 $\times$ 1031 erg s-1 for a distance of 4.6 kpc. This
and the $\gamma$-ray luminosity $5\times 10^{33}$ erg s-1 (Swihart et al.,
2022) are typical for BWs.

In addition to follow-up optical broad-band observations to confirm the
optical variability of the companion, optical spectroscopy would be useful to
establish its spectral type and to measure its radial velocity curve, to
better constrain the parameters of the system.

Figure 9: Distribution of the day-side temperature $T_{\rm d}$ of BW
companions. 18 sources are included. For PSRs J0023+0923, J0251+2606,
J0636+5129, J0952$-$0607, J1124$-$3653, J1301+0833, J1544+4937, J1555$-$2908,
J1653$-$0158, J1810+1744, J1959+2048, J2051$-$0827, J2052+1219, J2241$-$5236
and J2256$-$1024, the day-side temperatures are calculated using the base
(night-side) and the irradiation temperatures ($T_{\rm n}$ and $T_{\rm irr}$)
from Mata Sánchez et al. (2023) (see their table A1 and references therein).
For PSR J1311$-$3430, Kandel & Romani (2023) obtained very low $T_{\rm n}$,
negligible in comparison with $T_{\rm irr}$, which was derived from the
irradiation luminosity and parameters of the orbit (i.e. $T_{\rm d}\approx
T_{\rm irr}$). The position of J1641 according to our fit of the OSIRIS data
is marked by the arrow. We also include the BW candidate ZTF J1406+1222
(Burdge et al., 2022). The best-fitting value of the base temperature of PSR
J0610$-$2100 is below the lowest temperature covered by the spectral models
(van der Wateren et al., 2022); thus, we excluded it from the sample.

## Acknowledgements

We thank the referee for their useful comments and suggestions which allowed
us to improve the manuscript. We also thank Daniel Mata Sanchez and David
Kaplan for providing the optical data and for useful discussions, and Andrei
Igoshev for his helpful comments. The 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. This work used data obtained with eROSITA telescope
onboard SRG observatory. The SRG observatory was built by Roskosmos in the
interests of the Russian Academy of Sciences represented by its Space Research
Institute (IKI) in the framework of the Russian Federal Space Program, with
the participation of the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The
SRG/eROSITA X-ray telescope was built by a consortium of German Institutes led
by MPE, and supported by DLR. The SRG spacecraft was designed, built, launched
and is operated by the Lavochkin Association and its subcontractors. The
science data are downlinked via the Deep Space Network Antennae in Bear Lakes,
Ussurijsk, and Baykonur, funded by Roskosmos. The eROSITA data used in this
work were processed using the eSASS software system developed by the German
eROSITA consortium and proprietary data reduction and analysis software
developed by the Russian eROSITA Consortium. We acknowledge that CHIME is
located on the traditional, ancestral, and unceded territory of the
Syilx/Okanagan people. We are grateful to the staff of the Dominion Radio
Astrophysical Observatory, which is operated by the National Research Council
of Canada. CHIME is funded by a grant from the Canada Foundation for
Innovation (CFI) 2012 Leading Edge Fund (Project 31170) and by contributions
from the provinces of British Columbia, Québec and Ontario. The CHIME/FRB
project is funded by a grant from the CFI 2015 Innovation Fund (Project 33213)
and by contributions from the provinces of British Columbia and Québec, and by
the Dunlap Institute for Astronomy and Astrophysics at the University of
Toronto. Additional support is provided by the Canadian Institute for Advanced
Research (CIFAR), McGill University and the McGill Space Institute thanks to
the Trottier Family Foundation, and the University of British Columbia. The
CHIME/Pulsar instrument hardware is funded by the Natural Sciences and
Engineering Research Council (NSERC) Research Tools and Instruments (RTI-1)
grant EQPEQ 458893-2014. The work of DAZ and AVK was supported by the Russian
Science Foundation, grant number 22-22-00921,
https://rscf.ru/project/22-22-00921/. DAZ thanks Pirinem School of Theoretical
Physics for hospitality. SVZ acknowledges PAPIIT grant IN119323. Pulsar
research at UBC is funded by an NSERC Discovery Grant and by the Canadian
Institute for Advanced Research. F. A. D. is supported by the UBC Four Year
Fellowship. M. A. M. is supported by the NANOGrav NSF Physics Frontiers Center
award numbers 1430284 and 2020265 and NSF award number 2009425. J. K. S. is
supported by the NANOGrav NSF Physics Frontiers Center award numbers 1430284
and 2020265. 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.

## Data Availability

The optical data are available through the GTC data archive:
https://gtc.sdc.cab.inta-csic.es/gtc/, the CHIME data and the eROSITA data –
upon request.

## References

  * Abdollahi et al. (2022) Abdollahi S., et al., 2022, ApJS, 260, 53
  * Ablimit (2019) Ablimit I., 2019, ApJ, 881, 72
  * Abolfathi et al. (2018) Abolfathi B., et al., 2018, ApJS, 235, 42
  * Alpar et al. (1982) Alpar M. A., Cheng A. F., Ruderman M. A., Shaham J., 1982, Nature, 300, 728
  * Bellm et al. (2016) Bellm E. C., et al., 2016, ApJ, 816, 74
  * Benvenuto et al. (2014) Benvenuto O. G., De Vito M. A., Horvath J. E., 2014, ApJ, 786, L7
  * Benvenuto et al. (2015) Benvenuto O. G., De Vito M. A., Horvath J. E., 2015, MNRAS, 449, 4184
  * Bisnovatyi-Kogan & Komberg (1974) Bisnovatyi-Kogan G. S., Komberg B. V., 1974, Soviet Astronomy, 18, 217
  * Breton et al. (2012) Breton R. P., Rappaport S. A., van Kerkwijk M. H., Carter J. A., 2012, ApJ, 748, 115
  * Burdge et al. (2022) Burdge K. B., et al., 2022, Nature, 605, 41
  * CHIME Collaboration et al. (2022) CHIME Collaboration et al., 2022, ApJS, 261, 29
  * CHIME/Pulsar Collaboration et al. (2021) CHIME/Pulsar Collaboration et al., 2021, ApJS, 255, 5
  * Cabrera-Lavers et al. (2014) Cabrera-Lavers A., Pérez-García A., Abril Abril M., Bongiovanni A., Cepa G., 2014, Canarian Observatories Updates CUps, 3, https://www.gtc.iac.es/instruments/osiris/media/CUPS_BBpaper.pdf
  * Chen et al. (2013) Chen H.-L., Chen X., Tauris T. M., Han Z., 2013, ApJ, 775, 27
  * Cordes & Lazio (2002) Cordes J. M., Lazio T. J. W., 2002, preprint,  (arXiv:astro-ph/0207156)
  * Draghis et al. (2019) Draghis P., Romani R. W., Filippenko A. V., Brink T. G., Zheng W., Halpern J. P., Camilo F., 2019, ApJ, 883, 108
  * Flewelling et al. (2020) Flewelling H. A., et al., 2020, ApJS, 251, 7
  * Foight et al. (2016) Foight D. R., Güver T., Özel F., Slane P. O., 2016, ApJ, 826, 66
  * Gaia Collaboration et al. (2023) Gaia Collaboration et al., 2023, A&A, 674, A1
  * Ginzburg & Quataert (2021) Ginzburg S., Quataert E., 2021, MNRAS, 500, 1592
  * Green et al. (2019) Green G. M., Schlafly E., Zucker C., Speagle J. S., Finkbeiner D., 2019, ApJ, 887, 93
  * Guo et al. (2022) Guo Y., Wang B., Han Z., 2022, MNRAS, 515, 2725
  * Halpern (2022) Halpern J. P., 2022, ApJ, 932, L8
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Kandel & Romani (2023) Kandel D., Romani R. W., 2023, ApJ, 942, 6
  * Kennedy et al. (2022) Kennedy M. R., et al., 2022, MNRAS, 512, 3001
  * Lampton et al. (1976) Lampton M., Margon B., Bowyer S., 1976, ApJ, 208, 177
  * Lange et al. (2001) Lange C., Camilo F., Wex N., Kramer M., Backer D. C., Lyne A. G., Doroshenko O., 2001, MNRAS, 326, 274
  * Lynch et al. (2018) Lynch R. S., et al., 2018, ApJ, 859, 93
  * Manchester (2017) Manchester R. N., 2017, Journal of Astrophysics and Astronomy, 38, 42
  * Manchester et al. (2005) Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993
  * Marley et al. (2021) Marley M. S., et al., 2021, ApJ, 920, 85
  * Mata Sánchez et al. (2023) Mata Sánchez D., et al., 2023, MNRAS, 520, 2217
  * McEwen & et al. (prep) McEwen A., et al. in prep., MNRAS
  * Morello et al. (2019) Morello V., et al., 2019, MNRAS, 483, 3673
  * Nice & Taylor (1995) Nice D. J., Taylor J. H., 1995, ApJ, 441, 429
  * Predehl et al. (2021) Predehl P., et al., 2021, A&A, 647, A1
  * Ray (2023) Ray P., 2023, Public List of LAT-Detected Gamma-Ray Pulsars, https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-Detected+Gamma-Ray+Pulsars, (accessed September 28, 2023)
  * Romani et al. (2012) Romani R. W., Filippenko A. V., Silverman J. M., Cenko S. B., Greiner J., Rau A., Elliott J., Pletsch H. J., 2012, ApJ, 760, L36
  * Romani et al. (2015) Romani R. W., Filippenko A. V., Cenko S. B., 2015, ApJ, 804, 115
  * Romani et al. (2021) Romani R. W., Kandel D., Filippenko A. V., Brink T. G., Zheng W., 2021, ApJ, 908, L46
  * Shaifullah et al. (2016) Shaifullah G., et al., 2016, MNRAS, 462, 1029
  * Shklovskii (1970) Shklovskii I. S., 1970, Soviet Ast., 13, 562
  * Smith et al. (2002) Smith J. A., et al., 2002, AJ, 123, 2121
  * Stovall et al. (2014) Stovall K., et al., 2014, ApJ, 791, 67
  * Sunyaev et al. (2021) Sunyaev R., et al., 2021, A&A, 656, A132
  * Swihart et al. (2022) Swihart S. J., Strader J., Chomiuk L., Aydi E., Sokolovsky K. V., Ray P. S., Kerr M., 2022, ApJ, 941, 199
  * Tonry et al. (2012) Tonry J. L., et al., 2012, ApJ, 750, 99
  * Verbunt et al. (2017) Verbunt F., Igoshev A., Cator E., 2017, A&A, 608, A57
  * Yao et al. (2017) Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29
  * Yap et al. (2019) Yap Y. X., Li K. L., Kong A. K. H., Takata J., Lee J., Hui C. Y., 2019, A&A, 621, L9
  * Zharikov et al. (2013) Zharikov S., Tovmassian G., Aviles A., Michel R., Gonzalez-Buitrago D., García-Díaz M. T., 2013, A&A, 549, A77
  * Zharikov et al. (2019) Zharikov S., Kirichenko A., Zyuzin D., Shibanov Y., Deneva J. S., 2019, MNRAS, 489, 5547
  * van Dokkum (2001) van Dokkum P. G., 2001, PASP, 113, 1420
  * van Straten et al. (2012) van Straten W., Demorest P., Oslowski S., 2012, Astronomical Research and Technology, 9, 237
  * van der Wateren et al. (2022) van der Wateren E., et al., 2022, A&A, 661, A57