© The Authors 2025

1 Introduction

Millisecond pulsars (MSPs) are formed in binary systems containing a slowly rotating neutron star (NS) that is eventually spun up to millisecond periods by heavy mass and angular momentum accretion from an evolving companion. At the end of the accretion process, during which these systems are observed as low mass X-ray binaries (LMXBs), these NSs become detectable as radio MSPs. In turn, the companion star is expected to become a white dwarf (WD) with a helium (He) core (Alpar et al. 1982; Stairs 2004; Tauris 2011; Ferraro et al. 2015).

Globular clusters (GCs) provide ideal environments for the formation of these objects. In fact, the dynamical interactions occurring in their very dense cores (10³–10⁶ pc⁻³) can promote the formation of peculiar stellar populations (i.e., exotic objects), including LMXBs, cataclysmic variables, and binaries suitable for recycling NSs into MSPs (Paresce et al. 1992; Bailyn 1995; Bellazzini et al. 1995; Ferraro et al. 2001; Pooley & Hut 2006; Freire et al. 2008). Consequently, the number of LMXBs, and therefore MSPs, per unit mass in the Galactic GC population is approximately 10³ times higher than in the Galactic field (Clark 1975; Hui et al. 2010; Turk & Lorimer 2013; Zhao & Heinke 2022). Currently, 350 pulsars are known in 46 GCs¹.

Although the standard evolutionary scenario has been widely supported, several exceptions have emerged. These include MSPs with massive CO WD companions (Tauris 2011; Pallanca et al. 2013) or double NS systems (Jacoby et al. 2006; Tauris et al. 2017; Lian et al. 2025b). Moreover, the synergy of mul-tiwavelength observations has been found to be crucial in the identification of candidate supermassive pulsars (see for example the M13F system found by Cadelano et al. 2020). Of particular interest are the so-called holy grail systems, expected to be formed exclusively in the extreme environments of GC cores through dynamical interactions. One such candidate is NGC 1851E, recently identified by Barr et al. (2024), which represents the first known binary system hosting a pulsar in orbit around a mass-gap companion, either an exceptionally massive NS or the lowest-mass black hole discovered to date.

Another class of MSPs deviating from the canonical evolutionary path includes the eclipsing systems with low-mass, nondegenerate companions (Pallanca et al. 2010; Roberts 2013; Mucciarelli et al. 2013; Cadelano et al. 2015; Kirichenko et al. 2024; Smirnov et al. 2025; Turchetta et al. 2025). These systems, collectively referred to as “spiders”, are typically subdivided into two main categories: black widows (BWs) and redbacks (RBs). The distinction between the two types of spiders lies in the mass of their companion stars. RBs have companions with masses ranging from 0.1–0.5 M_⊙, whereas BWs have companions with masses smaller than 0.05 M_⊙.

These systems are characterized by their small orbital eccentricities, tight orbits (orbital periods P ≲ 1 day), and small mass functions, indicating the presence of a low-mass companion. Moreover, they exhibit periodic eclipses of the radio signal for a significant portion of their orbit, typically during the superior conjunction of the pulsar, due to absorption or scattering of radio photons by ionized material between the emitter and the observer (Fruchter et al. 1988; Ferraro et al. 2001; Ray et al. 2012; Cadelano et al. 2015; Roberts 2013).

In this context, multiwavelength studies of these systems provide a powerful tool for probing their evolution and the physics of accretion processes and for understanding their role within their surrounding environments (Kaluzny et al. 2003; Edmonds et al. 2002; Pallanca et al. 2010, 2014; Li et al. 2014; Cadelano et al. 2015, 2017, 2020; Turchetta et al. 2025; Ettorre et al. 2025). In particular, the identification of the optical counterpart to MSPs allows us to explore how the dynamical interactions shape the formation and evolution of MSPs, and to study stellar evolution under extreme conditions (Ferraro et al. 2003b,c; Sabbi et al. 2003a,b; Mucciarelli et al. 2013; Ferraro et al. 2015).

In systems hosting WD companions, the comparison of observed properties with binary evolution models enables estimates of the companion masses and their cooling ages (Ferraro et al. 2003a; Cadelano et al. 2019, 2020). The latter, provide a more robust estimate of the system’s true age compared to the characteristic spin-down ages inferred from radio observations (see Tauris 2012). In particular, precise age determinations are also fundamental for understanding the spin evolution of MSPs and for constraining the physical mechanisms at play during the recycling process (van Kerkwijk et al. 2005). Finally, when combined with radio timing, optical studies can yield precise mass measurements for both components of the binary system (Ferraro et al. 2003a; Cocozza et al. 2006; Antoniadis et al. 2013; Mucciarelli et al. 2013). These measurements benefit significantly from the well-determined distances and extinction values typical of GCs, leading to smaller uncertainties compared to MSPs in the Galactic field. A very recent example of how the synergy between radio and optical observations enables a full characterization of these systems is presented in Lian et al. (2025a) for the MSP M53A.

In this work we search for the companion stars to five MSPs in the GC M3 (NGC 5272), whose timing solutions were recently updated by Li et al. (2024) using the Five-hundred-meter Aperture Spherical radio Telescope (FAST; Nan et al. 2011). M3 is an old, bright, and extensively studied GC located at a distance of 10 kpc from the Sun (e.g., Carretta et al. 2009; Dotter et al. 2010; Dalessandro et al. 2013). It is characterized by an intermediate metallicity ([Fe/H] = −1.5) and negligible reddening (E(B − V) = 0.01, Harris 2010), given its high galactic latitude (b = 78.71°; Vasiliev & Baumgardt 2021).

Currently, five pulsars are known to be present in M3, all of them in binary systems. Pulsars M3A to D were found by Hessels et al. (2007) using the Arecibo telescope, while M3E and M3F were found with FAST (Pan et al. 2021), and their timing solutions were obtained for the first time by Li et al. (2024). A sixth pulsar candidate, M3C, has not been confirmed. Table 1 and Table 2 summarize the main radio timing parameters for the five pulsars in M3, as reported by Li et al. (2024), along with the physical properties of the binary systems that we were able to characterize, following the methodology described in Sect. 3. The full list of radio timing parameters is available in Table 1 of Li et al. (2024).

The companion star of M3B was identified by Cadelano et al. (2019) as an extremely low-mass WD with a mass of only 0.19 M_⊙. They also proposed a He WD as the likely counterpart to M3D. However, for this system, they adopted the source coordinates from the timing solution of Hessels et al. (2007), which are affected by significantly larger uncertainties (0^′′.6 in RA and a much larger uncertainty of 14^′′ in Dec), compared to the more recent solution by Li et al. (2024), possibly hindering the identification of the true companion.

Building on the recent and precise phase-connected timing solutions provided by Li et al. (2024) for all the confirmed pulsars in M3, we created an extended HST dataset, spanning from the near-UV to the optical I band to search for their optical companions and infer their properties using updated stellar models. While the optical companion to M3B was previously character-ized by Cadelano et al. (2019), we revisited this system using new photometry and improved models to derive its parameters. For the remaining MSPs, this work presents the first systematic search and, where possible, characterization of their optical counterparts.

The paper is organized as follows: In Sect. 2, we present the HST dataset and describe the photometric analysis, including the calibration and astrometric procedures adopted during the data reduction. Sect. 3 outlines our approach for identifying the optical counterparts and deriving their physical properties and provides a detailed discussion of the results for the five MSPs. Finally, in Sect. 4 we summarize our main findings and conclusions.

2 HST observations and data reduction

To search for the optical counterparts, we employed an extended photometric dataset consisting of high-resolution HST images ranging from the near-UV (F275W) to the optical I band (F814W), obtained using the ultraviolet-visible (UVIS) channel of the Wide Field Camera 3 (WFC3) and Wide Field Channel (WFC) of the Advanced Camera for Surveys (ACS). Table 3 provides an overview of the images used in this work to build our photometric catalog, while Fig. 1 shows the footprints of the adopted HST observations.

The WFC3/UVIS dataset comprises 14 images acquired on 2012 May 15 with the F275W, F336W, and F438W filters, as part of proposal 12605 (PI: Piotto). The ACS/WFC data consist of observations taken across three distinct epochs. The first epoch includes four images obtained on 2004 May 9 in the F435W and F555W filters (proposal 10008, PI: Grindlay). The second set comprises ten images taken on 2006 February 20 with the F606W and F814W filters (proposal 10755, PI: Sarajedini). Finally, the third ACS/WFC dataset, from proposal 16298 (PI: Libralato), consists of five images in the F814W filter acquired on 2021 June 3. In total, we used 33 images spanning seven different filters and four observational epochs.

Our photometric reduction process was tailored to achieve high accuracy and completeness at faint magnitudes in crowded fields, an essential requirement for identifying the optical counterparts of MSPs within GC cores. The photometric analysis was performed independently for each of the four datasets following the general procedure outlined below.

In detail, the data reduction was carried out on the calibrated exposures that include the charge transfer efficiency (CTE) correction (i.e., _flc images). For the F275W, F336W, and F438W exposures, cosmic ray removal was performed using the lacosmic Python package (van Dokkum 2001). After correcting for the “pixel-area map (PAM)”, the photometric analysis was conducted using the DAOPHOT package (Stetson 1987). As a first step, we generated a source list using the FIND routine, selecting only those detections with peak counts exceeding a specified threshold above the local background noise. The detection threshold was set empirically for each dataset to effectively minimize contamination from background fluctuations and spurious detections. Specifically, we chose a 3σ threshold for the WFC3/UVIS dataset and the ACS/WFC dataset from proposals 16298 and 10755, while 5σ for the F435W and F555W images, where σ denotes the standard deviation of the local background. Using PHOTO we obtained the concentric aperture photometry of each source. For each image we selected a large number (∼ 200) of bright, unsaturated and isolated stars.

These stars were subsequently used to model the point-spread function (PSF) using the PSF routine. The resulting best-fit PSF models were then applied to all detected sources with ALLSTAR. To refine the photometry, we employed the ALLFRAME package (Stetson 1994), which enables simultaneous PSF-fitting across all available images within each dataset. This step required the construction of a master list containing all objects to be fitted. We used DAOMATCH and DAOMASTER to align the catalogs and produce, for each filter, a list of stars with their average magnitudes and associated uncertainties, computed from the dispersion of individual measurements. Only stars detected in at least half of the exposures +1 were retained. The resulting single-filter lists were then merged into a comprehensive master list using again DAOMATCH and DAOMASTER.

Both the WFC3/UVIS and the ACS/WFC images are char-acterized by significant geometrical distortions across the field of view (FoV). Therefore, we corrected the instrumental positions of the WFC3 dataset with the equations by Bellini & Bedin (2009) for the filter F275W, while we used equations by Meurer et al. (2003) for the ACS/WFC images.

We assigned a photometric quality flag to each star in the catalog by selecting sources based on the CHI and SHARP photometric parameters, following the procedure described in Ettorre et al. (2025). For the WFC3/UVIS dataset, we calibrated the instrumental magnitudes to the VEGAMAG photometric system by crossmatching our catalog with the HST UV Globular Cluster Survey (HUGS) catalog (Piotto et al. 2015; Nardiello et al. 2018), using the CataXcorr software (Montegriffo et al. 1995). We selected only stars with good photometric quality. We then computed the average magnitude difference between stars in common to both catalogs and adopted this value as a zeropoint correction to convert our instrumental magnitudes to the VEGA-MAG system. The same procedure was applied to images from HST proposal ID 10775, using the ACS Survey of Galactic Globular Clusters catalog (Sarajedini et al. 2007) as the photometric reference. For all remaining ACS/WFC images, we applied the appropriate transformation equations and zeropoints provided by Sirianni et al. (2005) to calibrate the instrumental magnitudes.

Using CataXcorr we transformed the instrumental coordinates to the absolute coordinate systems (α, δ) by crossmatching each catalog with the Gaia DR3 dataset (Gaia Collaboration 2023). The coordinate system used in the Gaia dataset is based on the International Celestial Reference System, enabling accurate comparison with MSP positions obtained from timing analyzes that employ solar system ephemerides, which are also referenced to the same celestial framework.

The catalogs derived from the photometric analysis of the four datasets were crossmatched and combined, producing the final catalog. This comprehensive catalog includes all stars detected in at least one filter and contains over 350 000 sources.

Finally, we took advantage of the long temporal baseline offered by the available multiepoch dataset (see Table 3) to perform a relative proper motion (PM) analysis. For this purpose, we selected F814W images from HST programs 10775 (PI: Sarajedini) and 16298 (PI: Libralato), which are separated by approximately 15 years. This combination was chosen as it maxi-mizes the number of stars, particularly the candidate companions to the five MSPs, for which PMs can be reliably measured. For each star detected in both F814W datasets, we determined its position at each epoch by averaging the coordinates from individual exposures. The uncertainties were estimated from the standard deviation of these measurements. We then calculated the displacements in the x and y directions as the difference between the mean positions at the two epochs. Dividing these displacements by the time baseline yielded the proper motions in pixels per year. Finally, we converted these values into mil-liarcseconds per year using the ACS/WFC pixel scale, thereby obtaining the relative PMs.

Fig. 1 Footprints of the FoVs covered by the HST observations listed in Table 3. The underlying stars come from the Gaia DR3 dataset (Gaia Collaboration 2023) centerd on M3, obtained using the Gaia G-band magnitude to scale the size of the data points. Each FoV is color-coded according to the GO proposal number, listed in Table 3. The pink dots represent the position of each MSP and their name is represented. The black cross marks the center of the cluster (RA=205°.548416, Dec=28°.377277) and the yellow circle is the core radius rc = 0.98 pc. Both values are taken from the 2023 version of the GC catalog provided by Baumgardt et al. (2020, but see also Vasiliev & Baumgardt (2021))2.

Fig. 2 Left panel: color-magnitude diagram in the (mF275W, mF275W–mF336W) plane. Black points represent stars from our catalog with good photometric quality. The optical counterparts to M3B, M3D, and M3F are marked with orange, yellow and green squares, respectively. Right panel: Same CMD as shown on the left, but zoomed in to focus on the regions around the counterparts. Error bars for each counterpart reflect the corresponding photometric uncertainties. Overplotted are the He WDs cooling tracks from Cadelano et al. (2020). These correspond to stellar masses, ordered from right to left, between 0.17 M⊙ and 0.46 M⊙, as indicated by the colorbar.

3 Searching and characterizing the optical counterparts

As it can be noticed from Table 1 and Table 2, the positions of the five MSPs obtained from radio timing analysis have extremely small uncertainties, with a maximum combined uncertainty of ∼0.^′′009 for M3E. Consequently, the dominant source of positional uncertainty likely arises from our astrometric calibration, specifically from the residuals associated with the polynomial transformations used to align our instrumental positions with the Gaia reference frame.

Therefore, for each binary MSPs in the cluster, we searched for optical counterparts by carefully inspecting all stars located within a circle of radius 0.^′′1 centered on the radio position. For each candidate star, we analyzed its location in the color–magnitude diagrams (CMDs) and investigated possible photometric variability that could be associated with the binary orbital period.

In the following sections, we present a detailed description of the results from our search and analysis of the five MSPs found in M3. We successfully identified and characterized the optical companions of three MSPs, M3B, M3D, and M3F, as described in details in Sect. 3.1. In all three cases, the optical source coincides precisely with the pulsar’s radio position, demonstrating excellent agreement between the radio timing positions and the astrometry of our catalog, and confirming its high astrometric accuracy. These counterparts align with the WD cooling sequence in the (m_F275W, m_F275W − m_F336W) CMD shown in Fig. 2. For these three counterparts, it was not possible to derive PMs, as they are only detected in WFC3/UVIS images, acquired at the same epoch. On the contrary, no plausible optical counterparts were found for M3A and M3E, as detailed in Sect. 3.2. While the photometric uncertainties of the three identified optical companions, M3D and M3F in particular, make them compatible with being either He or CO WDs, in the following analysis (Sect. 3.1) we limit our study to He WD models as the vast majority of binary MSPs are expected to form through the canonical evolution scenario, which involves a phase of stable mass-transfer from a Roche-Lobe filling red giant secondary, eventually producing a He-WD companion. We stress here that this choice is supported by the observed properties of both M3D and M3F (see Table 2). In detail, both objects are relatively fast-spinning, with spin periods of ∼5.4 and 4.4 ms, respectively. Such short spin periods are expected for MSPs with He WD companions, which typically evolve through a phase with stable Roche-lobe overflow over long timescales. This allows the NS to efficiently accrete material and be spun up to short periods. In contrast, MSPs with CO WDs generally originate from intermediate-mass X-ray binaries, in which the secondary star, in a post-red giant branch phase (likely an asymptotic giant branch star), transfers mass over a brief episode, leading to slower spin periods (>10 ms; Tauris et al. 2012). In addition, in the case of M3D, the long orbital period of ∼129 days further supports a He WD companion, as CO WD MSPs are expected to have shorter orbital periods (<40 days) (Tauris et al. 2012).

Fig. 3 Finding charts for the MSP M3B. The panels display _drc images in the following HST filters, from left to right: WFC3/UVIS F275W and F336W (proposal ID 12605), ACS/WFC F555W (proposal ID 10008), and F814W (proposal ID 10775). Each panel shows a 1′′ × 1′′ FoV centered on the radio position of the pulsar. The magenta cross marks the nominal position of the MSP, while the magenta circle (radius 0.′′1) indicates the search radius. North is up, and east is to the left.

3.1 Detected companions

3.1.1 M3B

The radio timing study performed by Li et al. (2024) suggests that M3B has a minimum companion mass of 0.21 M_⊙ and an orbital period of 34.0 hours. As mentioned in the Introduction, the optical companion to this MSP was previously identified and characterized by Cadelano et al. (2019) as an extremely low-mass He WD with a mass of M_COM = 0.19 ± 0.02 M_⊙, a cooling age of ${1.0}_{-0.3}^{+0.2}$. Gyr, and an effective temperature of (12 ± 1) × 10³ K.

In our search for the optical counterpart using our photometric catalog, we confirm the identification of the same source: a very blue object located precisely at the pulsar’s radio position, as illustrated in the finding chart in Fig. 3, and consistent with the WD cooling sequence in the CMD (see Fig. 2). We were able to measure its magnitude in the F275W, F336W, and F438W filters, which were subsequently used for comparison with theoretical models. No evidence of photometric variability was found for this object. Although the lack of variability could stem from the limited and uneven coverage of the orbital period in the dataset, it is worth mentioning that such variability is uncommon for degenerate companions. Moreover, in cases where it is observed, it is generally not driven by irradiation from the MSP, but rather by intrinsic processes within the WD, such as global stellar pulsations (Cocozza et al. 2006; Maxted et al. 2013; Kilic et al. 2015).

At this point, we derive the physical parameters of the promising optical counterpart candidate by comparing its magnitudes with He WD stellar evolution models. The models adopted here are the updated version of the models presented in Cadelano et al. (2019) and already used by Cadelano et al. (2020), covering a mass range between 0.17 and 0.47 M_⊙. While we recall the reader to that paper for a full description of the latter models, here we briefly report the main ingredients. Binary evolutionary tracks were computed using MESA (version 12 115; Paxton et al. 2011, 2013, 2015, 2018), following a procedure similar to that of Istrate et al. (2014, 2016). The initial binary configuration assumes a 1.4 M_⊙ NS, treated as a point mass, and a 1.1 M_⊙ donor star. All the available models were interpolated to create a fine grid in the mass range 0.17–0.46 M_⊙, effective temperature between ∼3000 and 130 000 K, and cooling ages extending up to the cluster age.

To quantitatively constrain the companion’s properties, such as mass, age, and effective temperature, we defined a multivariate Gaussian logarithmic likelihood. The general equation in the d-dimensional case is

$$\ln ℒ=-\frac{1}{2}\left[\Delta y^{\top }{\Sigma }^{-1}\Delta y+\ln \left|\Sigma \right|+d\ln (2\pi )\right],$$(1)

where ∆y is the vector of residuals between the observed and model quantities, defined in the (m₂₇₅, − m₂₇₅ m₃₃₆, m₃₃₆ − m₄₃₈) space for M3B, and in the (m₂₇₅, m₂₇₅ − m₃₃₆) space for M3D and M3F. The matrix Σ is the corresponding covariance matrix including photometric, distance modulus, and reddening uncertainties. This method allowed us to derive constraints on the key parameters of the companion. For each WD parameter, we determined the best estimates and associated uncertainties from the 0.16, 0.50, and 0.84 percentiles of the corresponding distributions. The values obtained for M3B, listed in Table 1, are consistent with those previously reported by Cadelano et al. (2019), within the associated uncertainties.

The likelihood distributions of the companion’s parameters for M3B are shown in the left-hand panel of Fig. 4. The age distribution peaks at 0.98 Gyr, and the mass distribution is narrow and symmetric around M_COM = 0.18 M_⊙, indicating that this parameter is well constrained and confirming the extremely low mass of the MSP’s companion star. The log T_eff distribution is also symmetric, with a clearly defined peak around log T_eff = 4.08 K.

By determining the companion mass and using the binary orbital parameters derived through radio timing analysis, we can also place constraints on both the NS mass and the orbital inclination of the binary system. This is possible because the masses of the binary components are related to the orbital parameters via the NS mass function:

$$\frac{{\left({\text{M}}_{\text{COM}}\sin i\right)}^3}{{\left({\text{M}}_{\text{COM}}+{\text{M}}_{\text{NS}}\right)}^2}=\frac{4{\pi }^2{\text{x}}_{\text{p}}^3}{{\text{GP}}_{\text{ORB}}^2}.$$(2)

Here, M_COM and M_NS are the masses of the companion star and the NS, respectively, i is the orbital inclination angle, G the gravitational constant, x_p is the projected semimajor axis, and P_ORB is the orbital period. One should notice that the right-hand side of the equation depends only on the binary orbital parameters derived from radio data analysis, making its value well constrained. In contrast, the left-hand side involves the measured companion mass alongside two completely unknown quantities: the NS mass and the orbital inclination angle.

In order to explore the combination of NS masses and orbital inclination angles able to reproduce the observed mass function, we employed the affine-invariant Markov Chain Monte Carlo (MCMC) ensemble sampler, emcee³ (Foreman-Mackey et al. 2013). We used a uniform prior for the distribution of cos i between 0 and 1, and a flat prior for M_NS, within the range of 1.2–2.0 M_⊙, consistent with the NS mass distribution found by Abbott et al. (2023). Additionally, we tested a bimodal Gaussian prior on M_NS in the range 1.0–2.5 M_⊙, as suggested by Antoniadis et al. (2016), with the corresponding results presented in Appendix A.

The best-fit value was determined as the 50th percentile of the inferred distribution, while the lower and upper uncertainties were defined as the 16th and 84th percentiles, respectively. A similar method was used by Cadelano et al. (2019) to constrain the NS mass and orbital inclination of M3B, yielding ${\text{M}}_{\text{NS}}={1.1}_{-0.3}^{+0.3}$ and $i={89}_{-26}^{+1}$. We show the updated posterior distributions shown in the right-hand panel of Fig. 4. For the NS mass we find ${\text{M}}_{\text{NS}}={1.22}_{-0.01}^{+0.02}{\text{M}}_{\odot }$, a slightly larger value than previously reported in Cadelano et al. (2019), but consistent within the uncertainties. The inclination distribution retains the truncated Gaussian shape described by Cadelano et al. (2019), clearly indicating an almost edge-on orientation. As such, we adopt the peak of the distribution as the best estimate and define the lower uncertainty as the 16th percentile of a symmetrized distribution, yielding $i={88.84}_{-6.77}^{+1.16}\text{deg}$.

Fig. 4 Left-hand corner plots: Constraints on the mass, cooling age, and effective temperature of the companion star to M3B are presented. The 1D histograms show the marginalized likelihood distributions for each parameter, with the solid orange and dashed black lines indicating the 50th, 16th, and 84th percentiles, respectively. These percentiles represent the best-fit values and associated uncertainties. In the 2D histograms, contours correspond to the 1σ, 2σ, and 3σ levels, with the best-fit values marked by orange points and lines. The derived values for mass, cooling age, and temperature are reported at the top of each 1D distribution panel. Right-hand corner plot: Constraints on the NS mass and the orbital inclination angle of M3B. The one-dimensional histograms represent the marginalized probability distributions for these two parameters, with the solid green and dashed black lines indicating the best-fit values and their associated uncertainties. The bottom-left panel shows the joint two-dimensional posterior probability distribution, with contours representing the 1σ, 2σ, and 3σ confidence levels.

3.1.2 M3D

This binary system appears to have experienced an unusual evolutionary path, as suggested by the 129 day orbital period derived from radio timing analysis by Li et al. (2024). Using the position derived from the timing solution by Hessels et al. (2007), Cadelano et al. (2019) proposed a He WD as the most plausible companion to M3D. However, the radio position of the system was previously affected by significant uncertainties. The improved timing analysis by Li et al. (2024) markedly enhanced the precision of the system’s coordinates, considerably reducing the errors in both RA and Dec relative to earlier measurements.

Using our photometric catalog along with the updated, more precise position, we were able to identify the real companion: a blue object lying along the WD cooling sequence in the (m_F275W, m_F275W − m_F336W) CMD, which is significantly fainter (m_F275W ∼ 24.1 mag) than the optical counterpart found for M3B (see Figs. 2 and 5). This object is different from the candidate companion identified by Cadelano et al. (2019), which is situated approximately 20^′′ away from our object. For this source, photometric measurements were obtained only in the F275W and F336W filters, with no evidence of variability detected in either band.

The left panel of Fig. 6 shows the corner plot displaying the 1D and 2D marginalized likelihood distributions for the companion star’s parameters, derived from comparison with theoretical models. A close inspection of the 1D and 2D likelihood distributions for the cooling age reveals a pronounced, narrow peak around 0.3 Gyr, followed by another peak at ∼4 Gyr, similar to the case of M13F in Cadelano et al. (2020). The bimodality in the cooling age distribution arises from a dichotomy in the cooling timescales due to the occurrence of diffusion-induced hydrogen-shell flashes in the envelopes of proto-WDs with masses M ≳ 0.20 M_⊙ (Althaus et al. 2001a,b; Istrate et al. 2014, 2016).

In particular, WDs having thinner hydrogen layers are able to cool much faster than those with thicker ones. The models adopted here predict that hydrogen-shell flashes occur only for WDs with masses above approximately 0.24 M_⊙ (see Cadelano et al. 2020). Accordingly, the first peak corresponds to masses exceeding 0.24 M_⊙, which undergo flashes and cool rapidly. In contrast, the secondary peak is linked to lower-mass WDs that avoid flashes and cool more slowly through stable hydrogen burning. In this regime, even small photometric uncertainties can make an observation consistent with both pre- and post-flash models. Hence, the bimodal age posterior arises because of the large uncertainties in the companion’s measured magnitudes (see Fig. 2). When selecting only the young peak of the age likelihood distribution (M_COM > 0.24 M_⊙), the corresponding companion mass is ${\text{M}}_{\text{COM}}={0.33}_{-0.07}^{+0.06}{\text{M}}_{\odot }$, the effective temperature log ${\text{T}}_{\text{eff}}={4.14}_{-0.06}^{+0.05}$ and $\text{age}={0.31}_{-0.11}^{+0.15}\text{Gyr}$, while the older solution (M_COM ≤ 0.24 M_⊙) is associated with a significantly lower mass companion ${\text{M}}_{\text{COM}}={0.20}_{-0.02}^{+0.02}{\text{M}}_{\odot }$, with an effective temperature log ${\text{T}}_{\text{eff}}={4.00}_{-0.03}^{+0.04}$ and $\text{age}={3.93}_{-1.18}^{+1.19}\text{Gyr}$. Due to the large photometric uncertainties, the range of possible solutions for M3D span almost the entire range of masses of the adopted model grid. Consequently, it is worth stressing here that the uncertainties on the derived mass, age, and temperature are mostly driven by the prior limits of the model grid rather than by the observational constraints.

Unfortunately, the pulsar spin-down age generally cannot provide a more accurate constraint on the system’s age, as it is known to be highly unreliable (Tauris 2012). This is due to its dependence on the pulsar’s intrinsic spin-down rate, which is difficult to determine for MSPs in GCs because of contamination from acceleration within the cluster potential (Prager et al. 2017). In this case, the situation is further complicated by the lack of a spin-down age measurement for M3D.

The moderate eccentricity (e ∼ 0.075) and long orbital period (P_orb ∼ 129 d) of M3D suggest that the system has undergone significant dynamical perturbations. These properties might be useful to understand whether M3D is an old system hosting a low-mass WD companion or a relatively younger binary with a more massive WD. Indeed, according to the WD mass–orbital period correlation found by Tauris & Savonije (1999) for binary radio pulsars, M3D’s long orbital period would favour a system with a more massive WD, hence supporting the young scenario. Moreover, as discussed by Hessels et al. (2007) and more recently by Li et al. (2024), binaries with orbital periods between ∼10 and 1000 days are often efficiently formed through exchange interactions between an isolated NS and a primordial binary. Such systems tend to survive preferentially in GCs with relatively low stellar densities, such as M3 (Hut et al. 1992; Sigurdsson & Phinney 1993). In this framework, M3D likely originated from such an exchange encounter, after which the subsequent mass-transfer phase circularized the orbit, as expected for most GC MSPs. A secondary exchange occurring after recycling would typically produce a system with much higher eccentricity and a more massive companion, which is not the case for M3D. Instead, its current mild eccentricity may be explained by cumulative perturbations caused by flyby encounters occurring after the exchange event, as expected in a relatively low-density cluster such as M3.

In this context, we can apply Equation (5) of Rasio & Heggie (1995), which estimates the age of a system given its eccentricity, for e ≳ 0.01. We assume a central 1D velocity dispersion of 8.2 km s⁻¹ and a central mass density of log ρ_c = 3.78 M_⊙ pc⁻³, both taken from the 2023 version of the Baumgardt & Hilker (2018) GC catalog. We estimate the stellar number density by assuming an average stellar mass of ∼0.3 M_⊙. This value was obtained by using the synthetic stellar population provided by the Padova CMD 3.8 web interface⁴, assuming an age of 12.5 Gyr (Dotter et al. 2010) and a metallicity of [Fe/H] = −1.5 (Harris 2010). This approach does not account for the effects of mass segregation within the cluster. As a result, the derived average mass likely underestimates the true value in the cluster core, where M3D is located. The resulting number density of stars is n = 2.01 × 10⁴ pc⁻³ and the age of M3D computed using the equation (5) of Rasio & Heggie (1995) is t_e=0.075 = 1.2 Gyr. It is important to note that this estimate refers to the lifetime of the binary system. Consequently, this value should be interpreted as an upper limit for the age of the WD, favouring the interpretation that the companion of M3D is younger than the age suggested by the second old peak in the probability distribution in Fig. 6, in agreement with what expected from the WD mass– orbital period relation by Tauris & Savonije (1999). However, although the estimated timescale offers useful insight, current uncertainties regarding the formation and evolutionary history of M3D still prevent a definitive determination of its age.

For both solutions corresponding to the two distinct peaks of the likelihood distributions shown in Fig. 6, we derived the NS mass and the binary inclination angle by following the same analysis described in Sect. 3.1.1. When using ${\text{M}}_{\text{COM}}={0.33}_{-0.07}^{+0.06}{\text{M}}_{\odot }$, we find ${\text{M}}_{\text{NS}}={1.72}_{-0.32}^{+0.24}{\text{M}}_{\odot }$ and $i={47.33}_{-7.13}^{+7.91}\text{deg}$. On the other side, in the case of a very low-mass and old companion the NS is less massive than in the previous case, with ${\text{M}}_{\text{NS}}={1.28}_{-0.05}^{+0.09}{\text{M}}_{\odot }$ corresponding to an edge-on orbit with $i={81.43}_{-7.88}^{+5.89}\text{deg}$. The posterior distributions or both the young and old solutions are shown in blue and magenta, respectively, in the right-hand panel of Fig. 6. All values derived for the physical properties of M3D are summarized in Table 2, while the results obtained assuming a bimodal Gaussian prior for M_NS are described in Appendix A.

Fig. 5 Same as Fig. 3 but for M3D.

Fig. 6 Left-hand corner plots: Same as Fig. 4 but for M3D. Here both the young and old solutions are shown, and they are marked with cyan and magenta lines, respectively. Right-hand corner plot: Constraints on the NS mass and orbital inclination angle of M3D, with posterior distributions shown in cyan for the young solution and magenta for the old solution. The derived values for the companion mass, cooling age, effective temperature, NS mass, and orbital inclination for the two solutions are listed in Table 2.

3.1.3 M3F

The timing solution derived by Li et al. (2024) for M3F indicates the presence of a companion star with a minimum mass of 0.17 M_⊙, in a 3 day orbit around the pulsar. Searching for the candidate counterpart in our catalog, we identified a very faint star (m_F275W ∼ 24.7 mag), lying along the WD cooling sequence in the (m_F275W, m_F275W − m_F336W) CMD (marked by the green square in Fig. 2). This source is detected exclusively in the F275W and F336W filters, and it is shown in the finding charts of Fig. 7. No sign of photometric variability was detected in the available filters.

We constrained the physical properties of the companion star through comparison with the updated He WD cooling tracks, and, as illustrated by the likelihood distributions of age, log T_eff, and M_COM shown on the left side of Fig. 8, we found a scenario similar to that of M3D. In particular, the 1D histogram of the cooling age shows a clear bimodal distribution, with a narrow peak around 6.2 Gyr, and a secondary smaller peak at ∼0.5 Gyr. As already explained in Sect. 3.1.2, this bimodality has a clear physical origin, linked to the onset of diffusion-induced hydrogen-shell flashes in proto-WD around M ∼ 0.20 M_⊙. However, the significant photometric uncertainties associated with this faint star affect the mass estimates, broadening the age distribution and preventing the distinction between the two possible age solutions. In this case, selecting only the young peak of the likelihood distribution in Fig. 8 (M_COM > 0.24 M_⊙), we find a companion with mass ${\text{M}}_{\text{COM}}={0.32}_{-0.06}^{+0.07}{\text{M}}_{\odot }$, effective temperature log ${\text{T}}_{\text{eff}}={4.07}_{-0.04}^{+0.06}$ and $\text{age}={0.51}_{-0.18}^{+0.29}\text{Gyr}$. The older solution (M_COM ≤ 0.24 M_⊙) is instead associated with a significantly lower mass companion ${\text{M}}_{\text{COM}}={0.21}_{-0.02}^{+0.02}{\text{M}}_{\odot }$, with an effective temperature log ${\text{T}}_{\text{eff}}={3.98}_{-0.02}^{+0.02}$ and $\text{age}={5.15}_{-1.19}^{+1.20}\text{Gyr}$. According to the Tauris & Savonije (1999) WD mass–orbital period relation, the ∼3 days orbital period of M3F would favour a system with a less massive WD, supporting the old scenario.

Similar to M3D, the large photometric uncertainties associated to M3F result in poorly constrained parameters, with the likelihood distributions mostly reflecting the model boundaries.

We performed the same MCMC analysis as described in Sect. 3.1.1, to simultaneously estimate the NS mass and orbital inclination. For the old solution with a companion mass of ${\text{M}}_{\text{COM}}={0.21}_{-0.02}^{+0.02}{\text{M}}_{\odot }$ this yielded ${\text{M}}_{\text{NS}}={1.72}_{-0.31}^{+0.22}{\text{M}}_{\odot }$ and $i={54.08}_{-8.35}^{+7.06}$ degrees. While the younger solution, corresponding to ${\text{M}}_{\text{COM}}={0.32}_{-0.06}^{+0.07}{\text{M}}_{\odot }$ resulted in a less inclined system, with $i={32.46}_{-4.34}^{+4.61}$ degrees and a slightly less massive NS, of mass ${\text{M}}_{\text{NS}}={1.70}_{-0.32}^{+0.25}{\text{M}}_{\odot }$.

The posterior distributions corresponding to the two distinct solutions are presented in the right-hand panel of Fig. 8, while the derived physical parameters for M3F are summarized in Table 2. As for M3D, it is clear that the significant photometric uncertainties impacting the companion mass estimate are preventing our ability to tightly constrain both the orbital inclination and the NS mass, in contrast to the well-determined solution obtained for M3B. Appendix A shows the results obtained for M_NS and i under the assumption of a bimodal Gaussian prior for the NS mass.

Fig. 7 Same as Fig. 3 but for M3F.

Fig. 8 Same as for Fig 6, but for M3F.

3.2 Nondetections

3.2.1 M3A

According to Li et al. (2024), this system is a BW with a minimum companion mass of 0.01 M_⊙ and an orbital period of 0.14 days. Using our catalog, we searched for potential companions within the defined search radius but found no stellar source at the pulsar’s position. As shown in the finding chart in Fig. 9, the nearest object is a main sequence star located ∼0.^′′2 away.

Given the excellent agreement between the radio positions and the astrometry of our catalog, we conclude that this star is unlikely to be the true companion of the MSP. It is more plausible that the actual optical counterpart is either too faint to be detected or obscured by the nearby bright star.

Therefore, we inspected the PSF-subtracted images obtained after running ALLSTAR, specifically focusing on those taken near the pulsar’s inferior conjunction. In fact, BW companions are known for their pronounced optical variability, with light curves typically displaying a distinct minimum and maximum driven by the heating of the companion’s side directly exposed to the pulsar’s flux. This variability generally peaks at the pulsar’s inferior conjunction and reaches a minimum at the superior conjunction (Pallanca et al. 2014; Cadelano et al. 2015). However, we did not detect any indication of the true counterpart in the residual images.

At this point we aimed to set an upper limit on the companion’s mass by determining an upper magnitude limit achievable at the pulsar position. Since radio timing identifies this system as a BW, we expect the companion to be extremely faint, with most of its emission predominantly in the red filters. Consequently, we estimated a magnitude upper limit from a PSF-subtracted F814W image from HST proposal ID 10775, taken near inferior conjunction. We placed a circular aperture with a radius of 0.5^′′ centered on the radio position of M3A and measured the local number counts. The 3σ upper limit on the total counts within the aperture was computed from the mean and standard deviation of the pixel values inside the aperture. This value was normalized by the exposure time and converted into an instrumental magnitude, which was then calibrated using the filter’s zero point. The resulting 3σ upper limit on the calibrated magnitude is m_F814W,up = 26.06 ± 0.02 mag. This represents the maximum brightness of any undetected optical source at the MSP position.

To constrain the mass of a potential companion, we compared this photometric upper limit with theoretical stellar models. Specifically, we adopted a PARSEC isochrone (Bressan et al. 2012) with an age of 12.5 Gyr (Dotter et al. 2010) and a metallicity of [Fe/H] = −1.5 (Harris 2010), consistent with the properties of the cluster. The isochrone was transformed into the observational plane assuming a color excess of E(B − V) = 0.01, a distance of 10.08 kpc (Baumgardt & Vasiliev 2021), and an extinction coefficient of A_F814W/AV = 0.59696 (Cardelli et al. 1989; O’Donnell 1994). By comparing the calibrated limiting magnitude at the pulsar position with the interpolated isochrone, we derived an upper limit on the companion mass of M_COM,up = 0.108 ± 0.001 M_⊙.

Fig. 9 Same as Fig. 3 but for M3A.

3.2.2 M3E

MSP M3E is in an orbit of 7.1 days with a minimum companion mass of 0.2 M_⊙. The timing solution for this system was presented for the first time by Li et al. (2024). Starting from the radio position, we found no stellar source within 0.^′′1 in any filter, with the exception of the images from proposal IDs 10775 (PI: Sarajedini) and 16298 (PI: Libralato).

In both datasets, a source is located very close to the radio position, at a distance of ∼0.^′′03. This object is detected only in the F814W band images, suggesting it is a red source, predominantly emitting at longer wavelengths. However, as shown in the finding chart in Fig. 10, its position clearly shifts between the two epochs, separated by approximately 15 years.

As anticipated in Sect. 2, these two datasets were used to measure relative stellar PMs. This allowed us to determine the PM of the object located at the position of M3E. The resulting vector point diagram is shown in Fig.11, where gray points correspond to stars detected in both epochs, and the red dot indicates the motion of the object at the M3E radio position, including its associated uncertainty. Its measured motion (∆x = 0.091 ± 0.060 mas/yr; ∆y = 4.567 ± 0.044 mas/yr) lies well outside the cluster’s proper motion locus, strongly suggesting that this object is clearly not a cluster member.

According to Li et al. (2024), the proper motions and distances of the pulsars in M3 have not been individually measured. However, their values are expected to be consistent with those of the cluster itself, as suggested by the remarkably similar dispersion measures of the pulsars. Consequently, if M3E is a cluster member, the detected source near its radio position cannot be its true optical counterpart. Instead, its significant relative proper motion and red color strongly indicate that it is a foreground object unrelated to the cluster.

As in the case of M3A, we attempted to estimate an upper limit on the magnitude of the undetected companion. Given the uncertain nature of the object, we translated this constraint into an upper limit on the companion’s mass by comparing the magnitude limit with both a PARSEC isochrone, assuming the companion is a nondegenerate star, and the binary evolution models from Cadelano et al. (2020), under the alternative hypothesis that the companion is a WD. To do so, we used WFC3/UVIS F275W, F336W and F438W PSF-subtracted images from HST program GO-12605, adopting the same approach described in Sect. 3.2.1.

Comparing the magnitude upper limit obtained for the three filters (m_F275W,up = 25.55 ± 0.01 mag, m_F336W,up = 24.26 ± 0.01 mag and m_F438W,up = 24.12 ± 0.01 mag) with the PARSEC isochrone described in Sect. 3.2.1, we inferred a corresponding upper limit on the companion mass of M_COM,up = 0.484 ± 0.001 M_⊙. Under the WD hypothesis, we selected from the Cadelano et al. (2020) model grid only those models that predict magnitudes fainter than the observed upper limits in all three WFC3 filters. Each consistent model was assigned equal probability, and the posterior distribution for the companion mass was derived using weighted percentiles, adopting the 84th percentile as the mass limit. In the WD scenario, because more massive WDs are intrinsically fainter at a given cooling age, an undetected companion must be fainter, and therefore more massive, than our detection limit, resulting in a lower limit on its mass. Therefore we obtained M_COM,lower = 0.390 ± 0.001 M_⊙.

Fig. 10 Finding charts for the MSP M3E. The panels show the _drc images in the ACS/WFC F814W filter from HST proposal ID 10775 (left) and ID 16298 (right). Each panel spans 1′′ × 1′′ and is centered on the pulsar position.

Fig. 11 Vector point diagram showing the proper motion displacements in x and y for all stars common to the HST datasets from proposals 10775 and 16298 (gray points). The red dot marks the position and uncertainty of the object at the M3E radio position. The 1σ, 2σ, and 3σ contours, representing the modeled cluster PM distribution, are shown for reference. The distribution was modeled with a two-dimensional Gaussian mixture model (GMM) using a single component, implemented via the XDGMM Python framework (Holoien et al. 2017).

4 Conclusions

In this work, we carried out a detailed search for the optical counterparts to five binary MSPs located in the GC M3 (NGC 5272). Our analysis builds upon the recent high-precision, phase-connected timing solutions presented by Li et al. (2024), which provide accurate pulsar positions and orbital parameters. We built an extensive HST photometric dataset spanning from near-ultraviolet (F275W) to the optical I-band (F814W) to carry out an in-depth multiband investigation of each MSP. For each system, we searched for optical sources within a 0.^′′1 radius of the timing position, a conservative threshold that accounts for HST astrometric uncertainties. Candidate counterparts were evaluated based on: (i) positional coincidence with the pulsar, (ii) their location in CMDs, and (iii) comparison with updated binary evolution models (see Cadelano et al. 2020).

For M3A and M3E, no optical counterparts were identified, likely due to crowding by nearby bright stars or intrinsic faintness of the companions. To set upper limits on the companion mass, we estimated a magnitude upper limit at the pulsar positions using PSF-subtracted images. For M3A, a BW system, we derived an upper companion mass limit of M_COM,up = 0.108 ± 0.001 M_⊙ by comparing the limiting magnitude to a suitable PARSEC isochrone. In the case of M3E, a red source detected in F814W images at the radio position was ruled out as a foreground contaminant based on PM measurements. Given the uncertain nature of the companion, we used both PARSEC isochrones and the binary evolution models of Cadelano et al. (2020) to translate the magnitude limits into companion mass limits. We found M_COM,lower = 0.390 ± 0.001 M_⊙ under the WD hypothesis and M_COM,up = 0 484 ± 0.001 M_⊙ in the case of a nondegenerate companion.

On the other side, we identified and characterized the optical counterparts to three MSPs: M3B, M3D, and M3F. While the optical companion to M3B was previously studied by Cadelano et al. (2019), we reanalyzed this system using updated photometry and improved evolutionary models to refine its parameters. For M3D, the optical counterpart we identify is significantly offset (by ∼20^′′) from the candidate proposed by Cadelano et al. (2019). This discrepancy originates from the large positional uncertainties in both RA and Dec associated with the timing solution of M3D from Hessels et al. (2007), which was adopted in the earlier study. Finally, the timing solution for M3F was only recently obtained by Li et al. (2024), making this work the first attempt to identify its companion star.

In all three cases, the counterpart lies precisely at the radio position, demonstrating excellent agreement between optical and radio astrometry, and they exhibit photometric properties consistent with low-mass He WDs in NS binary systems. Based on model comparison, we confirm that M3B is a low-mass He WD with a relatively short cooling age of ${0.98}_{-0.19}^{+0.18}$ Gyr and a mass of ${\text{M}}_{\text{COM}}={0.18}_{-0.01}^{+0.01}{\text{M}}_{\odot }$. All the properties inferred for M3B are in agreement with the results presented by Cadelano et al. (2019).

For M3D and M3F, we report the first identification of their optical companions. In both cases, the counterparts are very faint blue stars whose magnitudes and colors are fully compatible with those of cooling He WDs. Both objects exhibit a bimodal cooling age distribution, with prominent peaks at ages ≲1 Gyr and ∼4–5 Gyr. This bimodality arises from a dichotomy in cooling timescales caused by diffusion-induced hydrogen-shell flashes. However, the significant photometric uncertainties associated with these two faint companions propagate into the physical parameter estimates, broadening the age distributions and hindering a clear distinction between the old and young solutions. When selecting only the young peak of the age likelihood distribution of M3D, the corresponding companion mass is ${\text{M}}_{\text{COM}}={0.33}_{-0.07}^{+0.06}{\text{M}}_{\odot }$, while the older solution is associated with a significantly lower mass companion ${\text{M}}_{\text{COM}}={0.20}_{-0.02}^{+0.02}{\text{M}}_{\odot }$. A similar result is obtained for M3F, with companion masses of ${\text{M}}_{\text{COM}}={0.32}_{-0.06}^{+0.07}{\text{M}}_{\odot }$ and ${\text{M}}_{\text{COM}}={0.21}_{-0.02}^{+0.02}{\text{M}}_{\odot }$ for the young and old solutions, respectively.

It is important to note here that there is a substantial difference in the meaning of the physical parameters derived for the three systems and more specifically of their uncertainties. For M3B, the tight constraints from both magnitude and color measurements locate the companion within a well-defined region of the model grid. Consequently, the inferred mass, cooling age, and effective temperature represent robust physical determinations with meaningful uncertainties. In contrast, the substantially larger photometric errors for M3D and M3F make these companions virtually consistent with nearly the full range of masses in the adopted model grid, as clearly illustrated in Fig. 2. The resulting posterior distributions for mass, age, and temperature and relative uncertaintes, derived in the two separate young and old solutions, are therefore mostly driven by the prior boundaries of our grid rather than by the constraining power of the data. This caveat should be kept in mind when interpreting the fitted values for these two systems.

Although our analysis does not fully break the degeneracy between the orbital inclination and the NS mass, we are able to reduce the range of plausible solutions by combining the derived companion mass with the binary’s orbital parameters. We employ an affine-invariant MCMC ensemble sampler, adopting a flat prior on the orbital inclination and a flat prior for M_NS, within the range of 1.2–2.0 M_⊙, consistent with the NS mass distribution found by Abbott et al. (2023). For M3B we found an almost edge-on configuration $\left(i={88.84}_{-6.77}^{+1.16}\text{deg}\right)$ with a NS mass of ${\text{M}}_{\text{NS}}={1.22}_{-0.01}^{+0.02}{\text{M}}_{\odot }$. For M3D and M3F, we inferred these two parameters separately for the young and old solutions, finding NS masses around ∼1.7 M_⊙ and low to moderately inclined systems with i between ∼30deg and ∼50 deg. The only exception is the solution found for M3D, under the old hypothesis, showing a low-mass NS $\left({\text{M}}_{\text{NS}}={1.28}_{-0.05}^{+0.09}{\text{M}}_{\odot }\right)$ and an edge-on inclination $\left(i={81.43}_{-7.88}^{+5.89}\text{deg}\right)$.

To resolve the degeneracy between NS mass and inclination angle, independent measurements of the companion mass, and ideally the orbital inclination, are required. A possible approach would be to determine the mass ratio of the binary components through spectroscopic measurements, but the companion’s very low luminosity places this beyond the reach of current optical telescopes. At the same time, the radio detection of the Shapiro time delay could yield precise masses for both the MSP and its companion star (Corongiu et al. 2012; Fonseca et al. 2021; Corongiu et al. 2023), especially for highly inclined systems such as M3B; however, precise timing and dense radio observations are required.

This study highlights the effectiveness of combining high-precision radio timing with deep, multiband HST imaging. The synergy between radio and optical observations, supported by comparisons with updated evolutionary models, enables the identification and characterization of MSP companions even under challenging observational conditions. However, as clearly demonstrated by the cases of M3D and M3F, future efforts to reach fainter magnitudes with significantly improved photometric precision are essential. Only by reducing uncertainties we can achieve a secure and detailed characterization of the faintest counterparts, as photometric errors play a crucial role in limiting our ability to constrain the nature of the entire system. This goal will become increasingly achievable with the advent of next-generation facilities such as the Extremely Large Telescope, whose MORFEO adaptive optics system and MICADO camera will deliver the sensitivity and angular resolution needed to overcome current limitations in the most crowded and faint observational regimes.

Acknowledgements

The research activities described in this paper were carried out with contribution of the Next Generation EU funds within the National Recovery and Resilience Plan (PNRR), Mission 4 – Education and Research, Component 2 – From Research to Business (M4C2), Investment Line 3.1 – Strengthening and creation of Research Infrastructures, Project IR0000034 – “STILES – Strengthening the Italian Leadership in ELT and SKA”. GE and ED also acknowledge financial support from the INAF Data analysis Research Grant (PI: E. Dalessandro) of the “Bando Astrofisica Fondamentale 2024”. CP acknowledges support from the INAF Large Grant 2022 “GCjewels” (PI: Andrea Possenti) approved with the Presidential Decree 30/2022. We thank the referee for their comments and suggestions, which have helped us improve the quality and clarity of the manuscript. G.E. is grateful to A. Della Croce for the useful discussions.

References

Appendix A Bimodal Gaussian prior on the NS mass

Figure A.1 shows the posterior distributions for the NS mass and the orbital inclination angle for for the three systems with identified optical counterparts: M3B, M3D, and M3F (from left to right). These results were obtained assuming a flat prior on the cosine of the inclination angle and a bimodal Gaussian prior on the NS mass, peaking at ∼1.4 M_⊙ and 1.8 M_⊙ and uniformly sampled between 1 and 2.5 M_⊙ (Antoniadis et al. 2016). This choice of prior differs from the analysis presented in the main text, where a flat prior between 1.2 and 2 M_⊙ was adopted, following the NS mass distribution inferred by Abbott et al. (2023) from gravitational wave observations (see Sect. 3.1.1).

On one side, the flat-prior approach removes any prior assumptions about the expected value of M_NS, allowing the data to explore a broader, physically plausible mass range, and yielding results that reflect a deliberately non-informative perspective. On the other hand, adopting a bimodal Gaussian prior enables us to test whether the observed properties of each system are consistent with a canonical NS mass distribution, while also providing an estimate of the orbital inclination under this assumption, restricting the parameter space to values consistent with the typical mass range observed in recycled MSPs. The results for the three systems are summarized in Table A.1.

Under the latter assumption, M3B shows a very narrow posterior on both the inclination angle and the NS mass, with the best-fit values consistent with the ones inferred under the assumption of a flat prior for M_NS, yielding a NS mass of ${1.22}_{-0.04}^{+0.04}{\text{M}}_{\odot }$ and a moderately high inclination angle $\left({88.75}_{-7.44}^{+1.25}\text{deg}\right)$. This confirms that the system is observed at high inclination angle and hosts a NS with a mass consistent with the typical values measured for recycled NSs in compact binary systems.

For M3D and M3F, we show the posterior distributions of M_NS and i corresponding to the two companion masses derived from the young and old solutions in cyan and magenta, respectively, consistently with Figs. 6 and 8 (see Sections 3.1.2 and 3.1.3).

For both systems, the posterior distributions are much broader with respect to the case of M3B, with the NS mass posterior distribution that still reflects the shape of the adopted bimodal prior. The only exception is the posterior distribution of M3D in the old companion scenario, constraining a nearly edge-on system with $i={83.27}_{-6.70}^{+4.67}\text{deg}$ and a NS mass of ${1.34}_{-0.05}^{+0.06}{\text{M}}_{\odot }$, similarly to to what is obtained assuming a flat prior on the NS mass. On the other side, the old solution of M3F corresponds to a relatively inclined system $\left({48.97}_{-4.00}^{+12.43}\text{deg}\right)$ and a NS with ${\text{M}}_{\text{NS}}={1.48}_{-0.12}^{+0.46}{\text{M}}_{\odot }$. Regarding the young solutions, the best-fit values indicate relatively low-inclination systems (${44.44}_{-5.40}^{+8.86}\text{deg}$ and ${30.70}_{-3.31}^{+5.14}\text{deg}$, for M3D and M3F respectively), hosting a less massive NS (${1.47}_{-0.11}^{+0.48}{\text{M}}_{\odot }$ and ${1.46}_{-0.10}^{+0.47}{\text{M}}_{\odot }$, respectively) compared to the results obtained assuming a flat prior on MNS (see Sect. 3.1.2 and 3.1.3).

These results once more confirm that the degeneracy in the physical parameters observed for both M3D and M3F significantly restricts our ability to constrain M_NS and the system’s orbital geometry. Hence, achieving higher photometric precision, and thereby obtaining more accurate measurements of the companion star’s properties, is crucial for improving the determination of the system’s properties.

Fig. A.1 Constraints on the NS mass and the orbital inclination angle of M3B, M3D and M3F, from left to right, assuming a flat prior for the inclination angle and a bimodal Gaussian prior for the NS mass, following Antoniadis et al. (2016). The one-dimensional histograms represent the marginalized probability distributions for these two parameters, with the solid blue and dashed black lines indicating the best-fit values and their associated uncertainties. The bottom-left panel shows the joint two-dimensional posterior probability distribution, with contours representing the 1σ, 2σ, and 3σ confidence levels. For M3D and M3F, the posterior distributions shown in cyan for the young solution and magenta for the old solution. The derived values for the companion mass, cooling age, effective temperature, NS mass and orbital inclination for the two solutions are listed in Table A.1.

All Tables

All Figures

Fig. 1 Footprints of the FoVs covered by the HST observations listed in Table 3. The underlying stars come from the Gaia DR3 dataset (Gaia Collaboration 2023) centerd on M3, obtained using the Gaia G-band magnitude to scale the size of the data points. Each FoV is color-coded according to the GO proposal number, listed in Table 3. The pink dots represent the position of each MSP and their name is represented. The black cross marks the center of the cluster (RA=205°.548416, Dec=28°.377277) and the yellow circle is the core radius rc = 0.98 pc. Both values are taken from the 2023 version of the GC catalog provided by Baumgardt et al. (2020, but see also Vasiliev & Baumgardt (2021))2.

In the text

Fig. 2 Left panel: color-magnitude diagram in the (mF275W, mF275W–mF336W) plane. Black points represent stars from our catalog with good photometric quality. The optical counterparts to M3B, M3D, and M3F are marked with orange, yellow and green squares, respectively. Right panel: Same CMD as shown on the left, but zoomed in to focus on the regions around the counterparts. Error bars for each counterpart reflect the corresponding photometric uncertainties. Overplotted are the He WDs cooling tracks from Cadelano et al. (2020). These correspond to stellar masses, ordered from right to left, between 0.17 M⊙ and 0.46 M⊙, as indicated by the colorbar.

In the text

Fig. 3 Finding charts for the MSP M3B. The panels display _drc images in the following HST filters, from left to right: WFC3/UVIS F275W and F336W (proposal ID 12605), ACS/WFC F555W (proposal ID 10008), and F814W (proposal ID 10775). Each panel shows a 1′′ × 1′′ FoV centered on the radio position of the pulsar. The magenta cross marks the nominal position of the MSP, while the magenta circle (radius 0.′′1) indicates the search radius. North is up, and east is to the left.

In the text

Fig. 4 Left-hand corner plots: Constraints on the mass, cooling age, and effective temperature of the companion star to M3B are presented. The 1D histograms show the marginalized likelihood distributions for each parameter, with the solid orange and dashed black lines indicating the 50th, 16th, and 84th percentiles, respectively. These percentiles represent the best-fit values and associated uncertainties. In the 2D histograms, contours correspond to the 1σ, 2σ, and 3σ levels, with the best-fit values marked by orange points and lines. The derived values for mass, cooling age, and temperature are reported at the top of each 1D distribution panel. Right-hand corner plot: Constraints on the NS mass and the orbital inclination angle of M3B. The one-dimensional histograms represent the marginalized probability distributions for these two parameters, with the solid green and dashed black lines indicating the best-fit values and their associated uncertainties. The bottom-left panel shows the joint two-dimensional posterior probability distribution, with contours representing the 1σ, 2σ, and 3σ confidence levels.

In the text

	Fig. 5 Same as Fig. 3 but for M3D.
In the text

Fig. 6 Left-hand corner plots: Same as Fig. 4 but for M3D. Here both the young and old solutions are shown, and they are marked with cyan and magenta lines, respectively. Right-hand corner plot: Constraints on the NS mass and orbital inclination angle of M3D, with posterior distributions shown in cyan for the young solution and magenta for the old solution. The derived values for the companion mass, cooling age, effective temperature, NS mass, and orbital inclination for the two solutions are listed in Table 2.

In the text

	Fig. 7 Same as Fig. 3 but for M3F.
In the text

	Fig. 8 Same as for Fig 6, but for M3F.
In the text

	Fig. 9 Same as Fig. 3 but for M3A.
In the text

	Fig. 10 Finding charts for the MSP M3E. The panels show the _drc images in the ACS/WFC F814W filter from HST proposal ID 10775 (left) and ID 16298 (right). Each panel spans 1′′ × 1′′ and is centered on the pulsar position.
In the text

Fig. 11 Vector point diagram showing the proper motion displacements in x and y for all stars common to the HST datasets from proposals 10775 and 16298 (gray points). The red dot marks the position and uncertainty of the object at the M3E radio position. The 1σ, 2σ, and 3σ contours, representing the modeled cluster PM distribution, are shown for reference. The distribution was modeled with a two-dimensional Gaussian mixture model (GMM) using a single component, implemented via the XDGMM Python framework (Holoien et al. 2017).

In the text

Fig. A.1 Constraints on the NS mass and the orbital inclination angle of M3B, M3D and M3F, from left to right, assuming a flat prior for the inclination angle and a bimodal Gaussian prior for the NS mass, following Antoniadis et al. (2016). The one-dimensional histograms represent the marginalized probability distributions for these two parameters, with the solid blue and dashed black lines indicating the best-fit values and their associated uncertainties. The bottom-left panel shows the joint two-dimensional posterior probability distribution, with contours representing the 1σ, 2σ, and 3σ confidence levels. For M3D and M3F, the posterior distributions shown in cyan for the young solution and magenta for the old solution. The derived values for the companion mass, cooling age, effective temperature, NS mass and orbital inclination for the two solutions are listed in Table A.1.

In the text