© ESO 2022

1. Introduction

Amongst the radio pulsars residing in binaries, black widows are a class of binary millisecond pulsars with very low-mass binary companions (M_c ≲ 0.05 M_⊙) in compact orbits (orbital period P_b ≲ 1 day) (Roberts et al. 2012; Chen et al. 2013). These binary systems display a set of observational phenomena not seen in the more typical binary millisecond pulsars with low-mass (0.1 M_⊙ <  M_c <  0.5 M_⊙) white dwarf companions in wider orbits (P_b >  1 day). The most common phenomenon in black widows is the presence of radio eclipses, where the pulsed radio emission from the pulsar disappears over a fraction of the orbit as the binary companion moves close to the line of sight of the pulsar. However, radio eclipses are not seen in all black widows (e.g., Breton et al. 2013; Stovall et al. 2014). For those black widow systems whose binary companions are detected at optical wavelengths, the companions display large sinusoidal variations in their brightness and colour, which is in phase with the binary orbit (Stappers et al. 2001; Breton et al. 2013). High precision pulsar timing of the millisecond pulsars in black widow systems reveal that some show variability of the orbital period on the timescales of months to years (e.g., Arzoumanian et al. 1994; Ng et al. 2014; Shaifullah et al. 2016). Members of the class of redback binary millisecond pulsars display these phenomena as well, typically showing longer radio eclipses and larger orbital variations; these systems can be clearly distinguished from black widow systems due to their higher companion masses (M_c ≃ 0.1−0.4 M_⊙; Roberts et al. 2012; Chen et al. 2013).

All of these observational phenomena are a result of the combination of the energetic wind from the radio pulsar (Ė ∼ 10³⁴⁻³⁵ erg s⁻¹) and the compactness of the binary orbit. This wind irradiates the hemisphere of the binary companion facing the pulsar, leading to the observed optical and near-infrared variations. With the companions close to filling their Roche lobes, they can undergo wind-driven mass loss, leading to the observed radio eclipses. Finally, the orbital period variations are explained through magnetically, rotationally and/or tidally induced quadrupole distortions, leading to coupling between the binary orbit and the pulsar spin (Applegate 1992; Applegate & Shaham 1994; Voisin et al. 2020a).

High precision pulsar timing of radio millisecond pulsars provides the possibility of detecting nano-Hertz gravitational waves with pulsar timing arrays (PTAs, Foster & Backer 1990). The sensitivity of PTAs strongly depends on the number of millisecond pulsars that can be included (Siemens et al. 2013). Hence, black widow systems are observed by the European Pulsar Timing Array (EPTA; Desvignes et al. 2016), the North American Nano-Hertz Observatory for Gravitational Waves (NANOgrav; Arzoumanian et al. 2018), and the Parkes Pulsar Timing Array (PPTA; Kerr et al. 2020) to determine whether their suitability is limited by orbital period variations (e.g., Bak Nielsen et al. 2020).

Many new radio millisecond pulsars in black widow systems have been discovered by targeting γ-ray sources discovered by the Fermi Gamma-Ray Space Telescope since 2009 (Ray et al. 2012). Hence, long pulsar timing baselines are not yet available for these systems. The first two black widow systems discovered before Fermi, PSRs B1957+20 and J2051−0827, were discovered in 1988 (Fruchter et al. 1988) and 1995 (Stappers et al. 1996) and both exhibit orbital period variations (Arzoumanian et al. 1994; Doroshenko et al. 2001), radio eclipses (Fruchter et al. 1988; Stappers et al. 1996) as well as optical variability of their binary companions (Stappers et al. 1999; Reynolds et al. 2007). The third system, PSR J0610−2100 discovered by Burgay et al. (2006) in 2003, is a 3.86 ms pulsar orbiting a M_c ≈ 0.02 M_⊙ companion in a 6.86-h orbit. While the system shows optical variations (Pallanca et al. 2012), high precision pulsar timing obtained to date does not report orbital variations or show evidence of radio eclipses (Espinoza et al. 2013; Desvignes et al. 2016).

In this paper, we present radio observations of PSR J0610−2100 obtained with several radio telescopes. We combine the radio observations with optical observations of the companion of PSR J0610−2100 from the European Southern Observatory (ESO) to constrain the properties of this black widow system. The observations are described in Sect. 2 and the data analysis is presented in Sect. 3. We present our results in Sect. 4 and discuss and conclude in Sect. 5.

2. Observations

2.1. Radio observations

The radio observations presented in this paper were obtained with the 76-m Lovell telescope at the Jodrell Bank Observatory in the United Kingdom, the 94-m diameter equivalent Nançay Radio Telescope in France, the 64-m Parkes radio telescope in Australia, and the 14 × 25-m diameter Westerbork Synthesis Radio Telescope in the Netherlands. Together, these datasets cover a 16-year time span from the discovery of PSR J0610−2100 in 2003 May (Burgay et al. 2006) until 2019 June. The Parkes observations were obtained with an Analogue Filterbank (AFB) at a centre frequency of 1390 MHz and the Parkes Digital Filterbanks (PDFB3 and PDFB4) at 1369 MHz and consist of pulse times-of-arrival (TOAs) previously presented in Burgay et al. (2006) and Espinoza et al. (2013). The Nançay observations consist of TOAs which have been published as part of the European Pulsar Timing Array (EPTA) Data Release 1.0 (DR1; Desvignes et al. 2016). The observations were obtained with the Berkeley-Orleans-Nançay (BON) instrument at centre frequencies 1400, 1600, and 2000 MHz, with bandwidths of 64 or 128 MHz.

These TOAs were supplemented with folded and dedispersed profiles (in the PSRCHIVE format; Hotan et al. 2004) from Parkes, Nançay, and the Lovell telescope. Parkes PDFB3 observations were retrieved from the CSIRO Data Access Portal. Additional Nançay observations were obtained with the newer Nançay Ultimate Pulsar Processing Instrument (NUPPI; Liu et al. 2014) at various L-band frequencies with a bandwidth of 512 MHz. The observations from the Lovell Telescope were obtained with a Digital Filterbank (DFB) since 2010 June and simultaneously with the Reconfigurable Open Architecture Computing Hardware backend (ROACH; Bassa et al. 2016) from 2011 April onwards. These instruments observe at centre frequencies of 1524 and 1532 MHz with bandwidths of 384 and 400 MHz, respectively. The DFB observations have been previously presented in EPTA DR1 (Desvignes et al. 2016). For our analysis, we used the DFB observations for the time period before 2011 April and ROACH observations afterwards, as the latter uses coherent dedispersion to minimise dispersive smearing.

Additionally, six observations from the Westerbork Synthesis Radio Telescope were used to create TOAs and search for radio eclipses at low frequencies. These observations were obtained with the PuMa II backend (Karuppusamy et al. 2008) at 345 MHz with 70 MHz of observing bandwidth. An overview of all observations is presented in Table 1. Figure 1a shows the observing frequency of the observations as a function of time.

Fig. 1. Radio timing observations obtained with the Lovell, Nançay, Parkes, and Westerbork telescopes. (a): timing observations as a function of observing frequency. (b) and (c): the residuals of the timing model as a function of orbital phase and time from observations with Lovell, Nançay, and Parkes. (d): the residuals of the timing model as a function of orbital phase from observations with Westerbork. The orbital phases covered by the WSRT observations is indicated with square brackets. Orbital phase ϕ = 0 is defined at the time of ascending node.

2.2. Optical observations

The faint companion (V ∼ 26.7) to PSR J0610−2100 was discovered by Pallanca et al. (2012) in archival observations obtained with FORS2, the Focal Reducer, and low dispersion spectrograph at the ESO VLT in Chile (Appenzeller et al. 1998), in 2004 December and 2005 January. Though only a few V (9 × 1010 s) and R-band images (20 × 590 s) were obtained, these observations showed that the companion was variable in brightness and that the modulation correlated with the orbital period and phase of PSR J0610−2100.

To improve the light curves, we have obtained follow-up RI¹ observations of the companion of PSR J0610−2100, again with FORS2, in 2010 October, November and December 2010 and 2011 April, with exposure times of 6 × 480 s and 12 × 540 s in R and 41 × 240 s in I.

In all observations, FORS2 was used with 2 × 2 binning, providing a pixel scale of $0\stackrel{″}{.}25$ pix⁻¹, with the target placed on the master chip. Exposure times were 1010 s in V, 590 s, 540 s, or 480 s in R, and 240 s in I. The seeing, as measured from the width of the point spread function (PSF), varied from excellent at $0\stackrel{″}{.}27$ to poor at $1\stackrel{″}{.}2$, though 80 per cent of the observations had a seeing better than $0\stackrel{″}{.}8$. For photometric calibration, short 30 s VRI exposures of the PSR J0610−2100 field and the Markarian A photometric standard field were obtained under photometric conditions during the 2010 October 15th observing run.

3. Analysis

3.1. Timing

We determined TOAs from the folded and dedispersed radio observations of PSR J0610−2100 using tools from the PSRCHIVE software package (Hotan et al. 2004). TOAs for the Westerbork observations were fully averaged in frequency and partially averaged in time to about 6 min per sub-integration. All other observations were fully averaged in time and frequency. After radio frequency interference (RFI) removal, all observations were manually inspected and non-detections and observations still badly affected by RFI were removed. We calculated TOAs by cross correlating the profiles against template profiles. Template profiles were created for each dataset, which we defined as a combination of instrument and centre frequency. Each template consisted of several summed (up to six) von Mises profiles with different phases, widths, and amplitudes to model the averaged profile.

We used the standard timing approach detailed in Edwards et al. (2006) where TEMPO2 (Hobbs et al. 2006) was used to convert the topocentric TOAs to the Solar-system bary-centre with the DE436 solar system ephemeris (Folkner & Park 2016) and placed them on the Terrestrial Time standard (BIPM2011; Petit 2010) using Barycentric Coordinate Time (TCB).

Due to arbitrary phase rotations between templates, as well as cable length differences between instruments, the pulsar timing model includes phase offsets, so-called jumps, between each dataset. We used the recommendation from Verbiest et al. (2016) to select the dataset with the lowest value of ${\overline{\sigma }}_{\mathrm{TOA}}/\sqrt{N_{\mathrm{TOA}}}$, where ${\overline{\sigma }}_{\mathrm{TOA}}$ is the median TOA uncertainty and N_TOA the number of TOAs in that dataset, as the reference dataset to reference jumps to. This is the NUPPI dataset at 1484 MHz. Due to the simultaneity of some of the Lovell DFB and ROACH observations, phase offsets between these datasets were determined directly and not included as free parameters in the timing model. Afterwards, TOAs from Lovell DFB observations that were obtained simultaneously with ROACH observations were removed. Phase offsets between the previously published Parkes PDFB3 TOAs and the additional PDFB3 TOAs from folded profiles were determined in a similar fashion, as there was an overlap between these datasets. After obtaining a phase offset, the TOAs from the folded profiles that coincided with the previously published PDFB3 TOAs were removed.

We started with the EPTA DR1 pulsar timing model of PSR J0610−2100 by Desvignes et al. (2016), which includes proper motion, a polynomial model for variations in dispersion measure (DM), and models the low-eccentricity binary orbit of PSR J0610−2100 with the ELL1 binary model (Lange et al. 2001). The ELL1 model removes the strong covariance between eccentricity e and longitude of periastron ω for small eccentricities by fitting the Laplace-Lagrange parameters η = esinω and κ = ecosω instead. TOAs from the Lovell, Nançay, and Parkes telescopes were included in the fit, while the Westerbork TOAs were only used to assess the presence of a radio eclipse or orbital phase dependent DM variations (see below). After an initial fit of the timing model using all datasets, we applied an iterative procedure where the profiles were refolded with the new model to create improved templates with which TOAs were redetermined. To minimise the impact of outliers on the fit, TOAs for which no profiles were available and which resulted in outlying residuals (>3σ) were removed. Additionally, any TOAs with uncertainties exceeding 20 ms were also removed.

After three such iterations, we obtained the timing model given in Table 2. The residuals for this model are displayed as a function of orbital phase and time in Figs. 1b and c. To allow for comparison with the Desvignes et al. (2016) timing model, we fitted a pulsar timing model using the same reference epoch and time of ascending node. Furthermore, we fitted a timing model where these are referenced to the centre of the timespan covered by the TOAs, to minimise covariances between fitted parameters. For the latter model, we included orbital period derivatives to fit for variations in the orbital period (P_b).

3.2. Orbital period variations

Since high precision pulsar timing of other black widow systems show variations of the orbital parameters (e.g., Shaifullah et al. 2016), we checked for variability by computing pulsar timing models for segments of the full timing baselines of 16 yr to constrain the presence of similar variations in PSR J0610−2100. Three different segment lengths of 180, 360, and 720 days were used, with 15, 30, and 60-day overlaps, respectively. For each segment we used TEMPO2 to propagate the timing model from Table 2 to a reference epoch at the centre of the segment, a step which adjusts the position, spin frequency, and DM for the effects of proper motion and spindown. The time of ascending node T_ASC was also adjusted to the centre of each segment. Next, we refitted the orbital parameters (orbital period P_b, T_ASC, the projected semi-major axis x, η, and κ) for each segment, keeping all other parameters fixed. These fits did not include Ṗ_b, which was set to zero.

Figure 2 shows the resulting variations in the orbital parameters for the different segment lengths compared to the parameters of the timing model in Table 2. These are discussed in Sect. 4.

Fig. 2. Constraints on orbital variations. Variations in the orbital period ΔPb (panel a), the epoch of the ascending node ΔTASC (panel b), and the projected semi-major axis of the orbit Δx (panel c) are plotted against time for 180-, 360-, and 720-day segments. Here ΔPb = 0 refers to the orbital period in Table 2 found for the full span and similarly for the other parameters. The right-hand panels show distributions of the variations for each segment length. The dotted lines denote the mean of the uncertainty on the measurements. The uncertainties are generally of the same order of magnitude as the standard deviations on the variations and the measurements do not show any clear trends.

3.3. Optical photometry

The FORS2 observations were corrected for bias using bias frames obtained during daytime and flat-fielded using twilight sky flats. To align the images all were registered using integer pixel offsets. A $2\stackrel{\prime }{.}1\times 2\stackrel{\prime }{.}1$ subsection of these images was analysed with DAOPHOT II (Stetson 1987), where instrumental magnitudes were determined using PSF fitting. Offsets in instrumental magnitude between the 30 s exposure VRI images and images obtained using the same filter were measured to place the magnitudes on the same zeropoint. The instrumental magnitudes were subsequently calibrated using 70 photometric standards from the Markarian A field using calibrated Johnson-Cousins VRI magnitudes by Stetson (2000). We fitted for zeropoint and colour terms, but used tabulated ESO extinction coefficients of 0.113, 0.109, and 0.087 mag airmass⁻¹ for VRI, respectively. The rms residuals of these calibrations were 0.09 mag in V, 0.05 mag in R and 0.11 mag in I.

3.4. Optical modelling

We modelled the optical photometry using the Icarus binary light curve modelling software (Breton et al. 2012)². Icarus generates a model star with a tessellated surface accounting for the tidal deformation in the binary potential. It then computes photometric fluxes in each filter band, at each orbital phase, for each surface element according to the local surface temperature, gravity, and viewing angle, and sums the fluxes from all surface elements to produce a model light curve. Here, we computed these fluxes from the model spectra of the Göttingen Spectral Library³ (Husser et al. 2013).

We used the Multinest nested sampling algorithm (Feroz et al. 2019) to explore the parameter space and estimate posterior distributions for model parameters. For exposures in which the optical counterpart to PSR J0610−2100 was not detected, we assumed a flux of zero and an uncertainty that was one third of the image’s 3σ detection threshold, and included these points in our fitting to ensure the model predicted fluxes below the detection limits around the minimum. We allow for 0.1 mag offsets in the flux calibrations in each band (guided by the residuals quoted in Sect. 3.3), which accounts for systematic uncertainties in both the photometric calibration, and in the continuum level of the model spectra from which Icarus computes fluxes.

The parameters of interest in our models were: the binary inclination angle, i, for which we assumed a prior probability of sin(i) between 0° (a face-on orbit) and 90° (edge-on); the Roche-lobe filling factor, f_RL, defined as the ratio between the stellar and Roche-lobe radii in the direction from the companion’s centre-of-mass to the L1 Lagrange point; the base temperature, T_base, of the star before irradiation and gravity darkening are applied; the irradiating temperature T_irr (see below); and the distance, d. For the distance, we assumed a prior that combines a model of the millisecond pulsar density within the Galaxy (Levin et al. 2013) as a function of the distance along the line of sight, multiplied by a log-Gaussian distribution centred on the DM distance of 3.26 kpc (see Sect. 4), with 1σ fractional uncertainties of ±45%.

In addition, we fitted for the extinction, E(B − V), and pulsar mass M_psr, neither of which are constrained by the data, but which were included here as ‘nuisance’ parameters and marginalised over to account for their effects on our parameter uncertainties. For the extinction, we adopted a Gaussian prior of E(B − V)=0.05 ± 0.02, based on the estimate from Green et al. (2019) for the extinction towards PSR J0610−2100, truncated such that E(B − V)> 0. For the pulsar mass, we assumed a uniform prior between 1.0 M_⊙ <  M_psr <  2.5 M_⊙, encompassing the observed neutron star mass range (Özel & Freire 2016). To account for gravity darkening, the surface temperature was modified according to the local surface gravity g, by T ∝ g^β, with a fixed index of β = 0.08, which assumes that the companion star has a convective envelope.

We also derived additional parameters using the sampled values of the optical model parameters, along with the parameters obtained from radio timing (see Table 2): the pulsar’s semi-major axis, x; orbital period, P_b; spin period, P; and intrinsic spin-down rate, Ṗ_int. The companion mass M_comp and binary mass ratio, q ≡ M_psr/M_comp were derived from M_psr, i, and the binary mass function $f(M)=M_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}^3{\sin }^{3}i/(M_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}}+M_{\mathrm{p}\mathrm{s}\mathrm{r}}{)}^2=4{\pi }^2{x}^3/(G{P}_{\mathrm{b}}^2)=5.2\times {10}^{-6}{M}_{\odot }$. The heating efficiency, ϵ ≡ L_irr/Ė_int(d) (Breton et al. 2013) is the ratio between the heating luminosity, $L_{\mathrm{i}\mathrm{r}\mathrm{r}}=4\pi A^2\sigma T_{\mathrm{i}\mathrm{r}\mathrm{r}}^4$ where A = x(1 + q)/sini is the orbital separation, and the pulsar’s intrinsic spin-down power budget, Ė_int(d) = 4π²IṖ_int(d)/P³ (where the distance dependence is due to the Shklovskii correction described below, and we assumed that the pulsar’s moment of inertia I = 10⁴⁵ g cm²). In our model, we assumed that the pulsar acts as a point source of an isotropic heating flux and that this flux raises the surface temperature when it impinges on the companion star. The efficiency ϵ therefore incorporates unknown beaming and stellar albedo factors. While ϵ >  1 is possible (e.g., if the heating flux is strongly beamed towards the companion), we excluded models requiring ϵ >  200, effectively ruling out large distances above which Ė_int becomes very low or negative.

4. Results

4.1. Timing

The pulsar timing dataset presented here has a timespan of 16 years, doubling that of Desvignes et al. (2016). In general, the pulsar timing solution derived from our datasets is consistent with that found by Desvignes et al. (2016). The few parameters that differ at the 3σ level are the astrometric parameters and the dispersion model. We attribute these differences to our longer dataset. The uncertainties on the rotational, astrometric, and binary parameters have decreased significantly. We note that the Laplace-Lagrange parameters yield an eccentricity of e = 2.3(4) × 10⁻⁵ with longitude of periastron ω = 79(10)°, also consistent with the values obtained by Desvignes et al. (2016).

Due to the large proper motion of PSR J0610−2100 of μ = 18.81(3) mas yr⁻¹, the pulsar is subject to an apparent acceleration, the Shklovskii effect (Shklovskii 1970) of cṖ_shk/P=μ²d which offsets the observed spin period derivative Ṗ_obs = 1.232661(10) × 10⁻²⁰ s s⁻¹. We can use this apparent acceleration to set an upper limit to the distance of the pulsar of d <  3.73 kpc, above which the intrinsic spin period derivative Ṗ_int = Ṗ_obs−Ṗ_shk would be negative and hence unphysical. This is consistent with the DM derived distance estimates of 3.54 kpc from the NE2001 model (Cordes & Lazio 2002) and 3.26 kpc from the YMW16 model (Yao et al. 2017), yielding Ṗ_shk = 1.17 × 10⁻²⁰ s s⁻¹ and 1.08 × 10⁻²⁰ s s⁻¹, respectively. Other apparent accelerations impacting the spin period derivative are caused by differential Galactic rotation, introducing Ṗ_dgr, and acceleration towards to the Galactic disk Ṗ_kz (Damour & Taylor 1991; Nice & Taylor 1995). For distances up to d = 3.73 kpc, the sum of these effects ranges from −1.8 × 10⁻²² to −1.2 × 10⁻²² s s⁻¹, two orders of magnitude smaller than that caused by the Shklovskii effect.

Similarly, these effects also impact the observed orbital period, with the Shklovskii effect dominating. For the above distance upper limit, we obtained Ṗ_{b, shk} < 7.91 × 10⁻¹⁴ s s⁻¹. The observed value of the orbital period derivative Ṗ_{b, obs} = −7(3) × 10⁻¹⁴ s s⁻¹ in our timing solution is not significant. In none of the models including higher order derivatives of P_b were any of the derivatives significant (>3σ). From refitting the orbital parameters for smaller time segments of the full time span, we found that the scatter in the measurements of P_b, T_ASC, and x are of the same order of magnitude as their uncertainties and, except for T_ASC, do not show any clear trends or periodic structure (Fig. 2). There is a possible weak trend seen in T_ASC during 2017 and 2018, but its significance is low (∼3σ). To quantify the variations, we obtained the standard deviation of the normalised variations on the shortest timescales as σ_ΔPb/P_b = 3.3 × 10⁻⁸, σ_ΔTasc/P_b = 5.1 × 10⁻⁶, and σ_Δx/x = 3.2 × 10⁻⁵. Hence, we conclude that no significant orbital variations of P_b and x are present in our observations of PSR J0610−2100 while slight variations in T_ASC may be present.

Previously, no eclipses of the radio signal have been detected in observations of PSR J0610−2100 between frequencies of 1374 and 2048 MHz. Confirming previous findings, our observations show no evidence of radio eclipses near superior conjunction at orbital phase ϕ = 0.25 (Fig. 1d) down to observing frequencies of 310 MHz. The residuals near superior conjunction do not show additional dispersive delays due to a ionised material evaporating from the pulsar companion. The rms in the residuals between the orbital phases of ϕ = 0.2 and 0.3 sets a limit on additional DM delays of 0.02 pc cm⁻³ at 1532 MHz and 0.005 pc cm⁻³ at 345 MHz.

4.2. Optical modelling

The best-fitting model of the optical light curve for the companion star is shown in Fig. 3 with posterior distributions for the fit and derived parameters shown in Fig. 4. The best-fitting model has a reduced chi-square of χ²/n = 98.3/82 = 1.20. All quoted uncertainties here are 95% confidence intervals from the posterior distributions estimated by MultiNest.

Fig. 3. Optical light curve and Icarus model for the companion star of PSR J0610−2100. Vertical lines show the measured magnitudes and their 1σ uncertainties. Where the counterpart is not detected in an image, we plot 2σ upper limits instead, denoted by downward-pointing arrows. Dashed curves show the flux in each band as predicted by our best-fitting Icarus model, solid curves show the same model after allowing for small (0.1 mag) offsets between the photometric calibration and the model normalisation. An orbital phase of ϕ = 0 corresponds to the pulsar’s ascending node; the companion star’s superior conjunction is at phase ϕ = 0.75.

Fig. 4. Posterior distributions for Icarus model parameters. Contours on the 2-d conditional distributions are at the 1σ, 2σ and 3σ levels. 1-d marginalised distributions are shown on the diagonal. The final parameter, ϵ, is derived from i, Tirr, d (which enters via the Shklovskii corrections to Ėint) and the pulsar mass (not shown), which is not constrained by the data, but which we marginalised over as a nuisance parameter.

From our Icarus models, we find that the base temperature T_base <  1800 K is very low, with the samples converging towards the lower limit allowed in our modelling of 1500 K. This limit is quite far below the lowest temperature (2300 K) covered by the spectral library used by Icarus to compute fluxes. Extrapolating to temperatures far below the range of the spectral libary is problematic, as unmodelled effects such as dust settling start to become important when simulating low-temperature spectra (Husser et al. 2013). While this unreliable extrapolation must be treated with caution, this will not badly affect our fits, as the model light curves are dominated by hotter surface elements on the irradiated side of the star. Nevertheless, the inferred temperature values far below the 2300 K limit of the spectral library cannot be trusted. We therefore conclude only that the base temperature must be well below the lower limit of the spectral library.

We find also that the irradiation temperature T_irr = 2820 ± 190 K is unusually low in this system, with T_irr >  4000 K being more typical for black widow companions (e.g., Breton et al. 2013, and see Fig. 5). However, this weak irradiation seems to be consistent with the pulsar’s low intrinsic spin-down power (1.9 × 10³³ erg s⁻¹ < Ė_int < 4.8 × 10³³ erg s⁻¹), as the heating efficiency 0.15 <  ϵ <  0.77 is typical for black widow systems. Owing to the low base temperature, this weak irradiation is still sufficient to produce a large relative temperature difference between the heated and unheated sides and hence the observed 4-mag variation in the I-band flux across the orbit.

Fig. 5. Companion properties and presence of radio eclipses of black widow systems for which the optical light curves of irradiated companions have been modelled with ICARUS. Higher irradiation luminosities and Roche-lobe filling factors are preferred for eclipsing black widow systems compared to non-eclipsing systems. The dashed line in the bottom panel denotes the companion radius for which physical eclipses would occur. The properties of PSR J1653−0158 are based on two different light curve models (see Nieder et al. 2020 for details).

The mechanism through which pulsars heat their companion stars remains unclear, but one explanation is that the pulsar’s gamma-ray flux is deposited in the upper stellar atmosphere of the companion. The heating efficiency can therefore be compared to the efficiency through which the pulsar converts its spin-down power into gamma-ray emission, η_γ ≡ 4πF_γd²/Ė_int(d), where F_γ is the observed gamma-ray flux and d is the best fitting distance of 2.24 kpc. The DR2 (Ballet et al. 2020) iteration of the Fermi-LAT Fourth Source Catalogue (Abdollahi et al. 2020) gives the gamma-ray flux of PSR J0610−2100 above 100 MeV as F_γ = (6.8 ± 0.5) × 10⁻¹² erg cm⁻² s⁻¹, which corresponds to 0.5 <  η_γ <  3.7 (where efficiencies above unity suggest non-isotropic emission preferentially beamed towards the observer). The ratio between ϵ and η_γ is better constrained than either individual efficiency (as the distance dependence in Ė_int cancels out), which we find to be 0.14 <  ϵ/η_γ <  0.34, suggesting that the pulsar’s gamma-ray flux has sufficient power to explain the irradiation in this system. However, we note that in the overall population of black widows, these efficiencies are not well correlated with one another (Draghis et al. 2019).

We find that the Roche-lobe filling factor is preferred to be below unity (0.72 <  f_RL <  0.94). We also performed the fitting with a fixed f_RL = 1.0, which gave parameter estimates that were consistent with those found when f_RL was free, albeit with slightly higher i and T_irr, but with slightly worse chi-square Δχ² = 3.8 and log-evidence ΔlogZ = 2.6⁴.

The observed light curve in Fig. 3 has hints of asymmetry, with data points after the conjunction at orbital phase ϕ = 0.75 being systematically lower than those at the equivalent pre-conjunction phase. Such asymmetries in black-widow light curves can be caused by various effects: reprocessing of heating flux by an asymmetric intra-binary shock (Romani & Sanchez 2016); channelling of charged particles onto the magnetic poles of the companion star (Sanchez & Romani 2017); or convective winds on the stellar surface (Kandel & Romani 2020; Voisin et al. 2020b). We therefore repeated our modelling with the addition of a hot spot (whose size, temperature, and location were new free parameters) on the trailing face of the companion star to mimic these effects. This model gave a better chi-squared, thanks to the additional degrees of freedom, with a slightly lower range of filling factors ($f_{\mathrm{RL}}=0.{78}_{-0.13}^{+0.1}$) and slightly higher inclination (i >  62°), but all parameters were consistent within the uncertainties of the simpler symmetric-heating model.

5. Discussion and conclusions

The radio timing measurements of PSR J0610−2100 show no evidence of radio eclipses down to low frequencies (>310 MHz) or significant orbital period variations over the 16 yr timing baseline. Modelling of the optical light curve of PSR J0610−2100’s low-mass companion confirms the irradiation of the companion by the pulsar, as previously reported by Pallanca et al. (2012), and shows that the companion is likely not filling its Roche lobe, while having a high inclination (54° < i <  89°). We found a surprisingly low temperature for the irradiated hemisphere of the companion.

While eclipses of the pulsed radio emission are seen in several black widow systems, the absence of radio eclipses in PSR J0610−2100 is not unique. We compared companion properties between all black widow systems modelled with ICARUS to investigate the relation between irradiation, Roche-lobe properties, and the presence of radio eclipses. Draghis et al. (2019) modelled optical light-curves of irradiated black widow companions with ICARUS for nine systems, of which no radio eclipses are observed in four of them (PSRs J0023+0923; Bak Nielsen et al. 2020, J0636+5128; Stovall et al. 2014, J0952−0607; Bassa et al. 2017, and J2241−5236; Keith et al. 2011). Radio eclipses have been observed in PSRs J0251+2606, J1124−3653, J1301+0833, B1957+20, and J2052+1219 (Deneva et al. 2021; Fruchter et al. 1988, Deneva & Ray, priv. comm.). Furthermore, light-curve modelling has also been performed for PSR J1810+1744 (Romani et al. 2021) and PSR J1653−0158 (Nieder et al. 2020), where the PSR J1810+1744 shows radio eclipses (Polzin et al. 2018), while PSR J1653−0158 remains undetected in radio, which we attribute to eclipses of the radio emission over large fractions of the orbit (Nieder et al. 2020).

Radio eclipses or orbitally modulated DM variations detected in the pulsar observations are caused by (ionised) material in the evaporative wind from the companion passing the line of sight. Conversely, the absence of eclipses can be caused by either the absence of ionised material, or by the material missing the line of sight. The mass lost in the evaporative wind is predicted to be higher for systems filling their Roche lobe, compared to detached systems (f_RL <  1) and dependent on the irradiation luminosity (Ginzburg & Quataert 2020). This is qualitatively confirmed in Fig. 5, which collates the companion properties and highlights the presence or absence of radio eclipses; the black widow systems showing radio eclipses tend to have a higher Roche-lobe filling factor and irradiation luminosity compared to the non-eclipsing systems. Material missing the line of sight is predominantly dependent on the orbital inclination and a statistical analysis of a larger sample of black widow pulsars indicates that pulsar mass functions for eclipsing black widow systems are higher than non-eclipsing systems, indicating that the former have higher inclinations (Guillemot et al. 2019). We quantified this in the lower panel in Fig. 5 where the minimum distance R_min between the companion and the line of sight is compared to the Roche-lobe radius R_RL. The distance R_min is a strong function of inclination through the orbital projection, which sets R_min = x(1 + q)/tani, as well as mass ratio q constrained by the inclination and the mass function (we assumed a pulsar mass of 1.5 M_⊙). Furthermore, the Roche-lobe radius depends on inclination through the mass ratio. For R_min/R_RL <  1 and (R_min/R_RL)/f_RL <  1 the Roche lobe and actual companion radius would actually pass the line of sight. We find that for those systems with radio eclipses, the eclipses are caused by material that extends multiple Roche-lobe radii away from the companion. PSR J0610−2100 is one of the few non-eclipsing systems where the companion passes within a few Roche-lobe radii from the line of sight. Hence, it is tempting to speculate that the low spin-down of PSR J0610−2100 and the correspondingly low temperature of the irradiated side of the companion, resulting in less mass lost in the evaporative wind, may explain the absence of radio eclipses observed in this system.

Where the original black widow pulsars show lMarge orbital variations leading to ΔT_ASC/P_b ≈ 6 × 10⁻⁴, ΔP_b/P_b ≈ 1 × 10⁻⁷, and Δx/x ≈ 9 × 10⁻⁴ in the case of PSR J2051−0827 (Shaifullah et al. 2016) and ΔP_b/P_b ≈ 1.6 × 10⁻⁷ for PSR B1957+20 (Arzoumanian et al. 1994), these variations were not detected in P_b and x for PSR J0610−2100, with limits of |ΔP_b|/P_b <  3.3 × 10⁻⁸ and |Δx|/x <  3.2 × 10⁻⁵. There is a possible weak trend (∼3σ) seen in T_ASC with ΔT_ASC/P_b ∼ −7 × 10⁻⁶. Hence, the orbital variations in PSR J0610−2100 are at least an order of magnitude smaller compared to those in PSRs B1957+20 and J2051−0827.

In the spider timing model by Voisin et al. (2020a), orbital variations are modelled by a gravitational quadrupole moment, which is split into spin J_s, tidal J_t, and time-dependent variable J_v(t) components. The spin and tidal components are strong functions of the Roche-lobe radius, set by the mass ratio q and Roche-lobe filling factor f_RL. We compared J_s and J_t for PSRs J0610−2100 and J2051−0827. Assuming a 1.5 M_⊙ pulsar mass and using the orbital inclinations of i = 73° and i ∼ 40° for PSRs J0610−2100 and J2051−0827 (Stappers et al. 2001) respectively, we used the mass function to obtain mass ratios. With f_RL = 0.86 for PSR J0610−2100 and f_RL = 1 for PSR J2051−0827 (Stappers et al. 2001), we found that these components are a factor 3.5–4 higher in PSR J2051−0827 compared to PSR J0610−2100. Unfortunately, we cannot compare the time-dependent variable components J_v(t), as these depend on the observed orbital variations and are not linked to physical quantities in Voisin et al. (2020a). It is unclear if the variable components scale similarly to the spin and tidal components of the quadrupole moment.

While the orbital period derivative obtained from timing itself is not a significant measurement, it is offset by the Shlovskii effect by Ṗ_{b, shk} < 7.91 × 10⁻¹⁴ s s⁻¹. At the best-fitting distance of 2.24 kpc distance from light-curve modelling, we found that Ṗ_{b, shk} = 4.75 × 10⁻¹⁴ s s⁻¹, leading to an intrinsic orbital period derivative of Ṗ_b,int = −12(3) × 10⁻¹⁴ s s⁻¹. This tentative detection is an order of magnitude larger than the orbital period derivative expected due to emission of gravitational radiation, for which general relativity predicts Ṗ_b,GR = −4.6 × 10⁻¹⁵ s s⁻¹ (Peters 1964). The presence of similar long-term orbital period derivatives in other black widow systems without showing detectable shorter term orbital period variations (e.g., PSR J0023+0923, Bak Nielsen et al. 2020 and PSR J0636+5128, Stovall et al. 2014) indicates that these may be the result of spin-orbit coupling as suggested by the models of Applegate & Shaham (1994) and Voisin et al. (2020a).

PSR J0610−2100 is currently timed as part of PTA experiments by the EPTA (Desvignes et al. 2016), NANOgrav Arzoumanian et al. (2018), and the PPTA (Kerr et al. 2020). Since the sensitivity of PTAs to gravitational waves is strongly dependent on the number of millisecond pulsars that can be timed to high precision (e.g. Siemens et al. 2013), continued timing of PSR J0610−2100 for PTA purposes is warranted. The absence of evidence for radio eclipses and significant orbital variations in the timing of PSR J0610−2100 confirms the suitability of this pulsar for PTAs. The long timing baseline of PSR J0610−2100 may provide clues whether the tentative detection of the orbital period derivative is a manifestation of spin-orbit coupling induced orbital variations that operate on time scales longer than those observed in systems such as PSRs B1957+20 and J2051−0827. This will help to assess the suitability of more recently discovered black widow systems for PTA purposes, even if optical observations show that the companion is irradiated by the millisecond pulsar.

Acknowledgments

We thank the anonymous referee for the careful reading of and useful comments on the manuscript. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 386.D-0207. This research was made possible by support from the Dutch National Science Agenda, NWA Startimpuls – 400.17.608. C.J.C. and R.P.B acknowledge support from the ERC under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 715051; Spiders). This work was supported by the Max-Planck-Gesellschaft (MPG). The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). We acknowledge financial support from ‘Programme National de Cosmologie et Galaxies’ (PNCG) of CNRS/INSU, France. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from The Netherlands Foundation for Scientific Research (NWO). The Parkes radio telescope is part of the Australia Telescope National Facility (grid.421683.a) which is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Wiradjuri people as the traditional owners of the Observatory site.

References

All Tables

All Figures

Fig. 1. Radio timing observations obtained with the Lovell, Nançay, Parkes, and Westerbork telescopes. (a): timing observations as a function of observing frequency. (b) and (c): the residuals of the timing model as a function of orbital phase and time from observations with Lovell, Nançay, and Parkes. (d): the residuals of the timing model as a function of orbital phase from observations with Westerbork. The orbital phases covered by the WSRT observations is indicated with square brackets. Orbital phase ϕ = 0 is defined at the time of ascending node.

In the text

Fig. 2. Constraints on orbital variations. Variations in the orbital period ΔPb (panel a), the epoch of the ascending node ΔTASC (panel b), and the projected semi-major axis of the orbit Δx (panel c) are plotted against time for 180-, 360-, and 720-day segments. Here ΔPb = 0 refers to the orbital period in Table 2 found for the full span and similarly for the other parameters. The right-hand panels show distributions of the variations for each segment length. The dotted lines denote the mean of the uncertainty on the measurements. The uncertainties are generally of the same order of magnitude as the standard deviations on the variations and the measurements do not show any clear trends.

In the text

Fig. 3. Optical light curve and Icarus model for the companion star of PSR J0610−2100. Vertical lines show the measured magnitudes and their 1σ uncertainties. Where the counterpart is not detected in an image, we plot 2σ upper limits instead, denoted by downward-pointing arrows. Dashed curves show the flux in each band as predicted by our best-fitting Icarus model, solid curves show the same model after allowing for small (0.1 mag) offsets between the photometric calibration and the model normalisation. An orbital phase of ϕ = 0 corresponds to the pulsar’s ascending node; the companion star’s superior conjunction is at phase ϕ = 0.75.

In the text

Fig. 4. Posterior distributions for Icarus model parameters. Contours on the 2-d conditional distributions are at the 1σ, 2σ and 3σ levels. 1-d marginalised distributions are shown on the diagonal. The final parameter, ϵ, is derived from i, Tirr, d (which enters via the Shklovskii corrections to Ėint) and the pulsar mass (not shown), which is not constrained by the data, but which we marginalised over as a nuisance parameter.

In the text

Fig. 5. Companion properties and presence of radio eclipses of black widow systems for which the optical light curves of irradiated companions have been modelled with ICARUS. Higher irradiation luminosities and Roche-lobe filling factors are preferred for eclipsing black widow systems compared to non-eclipsing systems. The dashed line in the bottom panel denotes the companion radius for which physical eclipses would occur. The properties of PSR J1653−0158 are based on two different light curve models (see Nieder et al. 2020 for details).

In the text