# Mass estimates from optical modelling of the new TRAPUM redback PSR J1910-5320

O.G Dodge, R.P Breton, C.J Clark, M Burgay, J Strader, K.-Y Au, E.D Barr, S Buchner, V.S Dhillon, E.C Ferrara, et al.

► To cite this version:

O.G Dodge, R.P Breton, C.J Clark, M Burgay, J Strader, et al. Mass estimates from optical modelling of the new TRAPUM redback PSR J1910-5320. Monthly Notices of the Royal Astronomical Society, 2024, 528 (3), pp.4337-4353. 10.1093/mnras/stae211. hal-04446890

HAL Id: hal-04446890

https://hal.science/hal-04446890

Submitted on 20 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Mass estimates from optical modelling of the new
TRAPUM redback PSR J1910-5320

O.G Dodge, R.P Breton, C.J Clark, M Burgay, J Strader, K.-Y Au, E.D Barr, S Buchner, V.S Dhillon, E.C Ferrara, et al.

► To cite this version:

O.G Dodge, R.P Breton, C.J Clark, M Burgay, J Strader, et al.. Mass estimates from optical modelling of the new TRAPUM redback PSR J1910-5320. Monthly Notices of the Royal Astronomical Society, 2024, 528 (3), pp.4337-4353. 10.1093/mnras/stae211. hal-04446890

HAL Id: hal-04446890

https://hal.science/hal-04446890v1

Submitted on 20 Apr 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mass estimates from optical modelling of the new TRAPUM redback PSR J1910−5320

O. G. Dodge Ⓓ,<fnref n="1" />★ R. P. Breton Ⓓ,<fnref n="1" /> C. J. Clark Ⓓ,<fnref n="2" />,<fnref n="3" /> M. Burgay Ⓓ,<fnref n="4" /> J. Strader,<fnref n="5" /> K.-Y. Au,<fnref n="6" /> E. D. Barr Ⓓ,<fnref n="7" />

S. Buchner,<fnref n="8" /> V. S. Dhillon Ⓓ,<fnref n="9" />,<fnref n="10" /> E. C. Ferrara,<fnref n="11" />,<fnref n="12" />,<fnref n="13" /> P. C. C. Freire Ⓓ,<fnref n="7" /> J.-M. Griessmeier,<fnref n="14" />,<fnref n="15" />

M. R. Kennedy Ⓓ,<fnref n="1" />,<fnref n="16" /> M. Kramer Ⓓ,<fnref n="1" />,<fnref n="7" /> K.-L. Li,<fnref n="6" /> P. V. Padmanabh Ⓓ,<fnref n="2" />,<fnref n="3" />,<fnref n="7" /> A. Phosrisom,<fnref n="1" /> B. W. Stappers Ⓓ,<fnref n="1" />

S. J. Swihart<fnref n="17" /> and T. Thongmeearkom Ⓓ,<fnref n="1" />,<fnref n="18" />

Accepted 2024 January 17. Received 2024 January 12; in original form 2023 October 17

## ABSTRACT

Spider pulsars continue to provide promising candidates for neutron star mass measurements. Here we present the discovery of PSR J1910−5320, a new millisecond pulsar discovered in a MeerKAT observation of an unidentified Fermi-LAT gamma-ray source. This pulsar is coincident with a recently identified candidate redback binary, independently discovered through its periodic optical flux and radial velocity. New multicolour optical light curves obtained with ULTRACAM/New Technology Telescope in combination with MeerKAT timing and updated SOAR/Goodman spectroscopic radial velocity measurements allow a mass constraint for PSR J1910−5320. ICARUS optical light curve modelling, with streamlined radial velocity fitting, constrains the orbital inclination and companion velocity, unlocking the binary mass function given the precise radio ephemeris. Our modelling aims to unite the photometric and spectroscopic measurements available by fitting each simultaneously to the same underlying physical model, ensuring self-consistency. This targets centre-of-light radial velocity corrections necessitated by the irradiation endemic to spider systems. Depending on the gravity darkening prescription used, we find a moderate neutron star mass of either $1.6 \pm 0.2$ or $1.4 \pm 0.2 M_\odot$. The companion mass of either $0.45 \pm 0.04$ or $0.43_{-0.03}^{+0.04} M_\odot$ also further confirms PSR J1910−5320 as an irradiated redback spider pulsar.

**Key words:** pulsars: general – pulsars: individual: J1910−5320 – techniques: photometric – techniques: radial velocities.

## 1 INTRODUCTION

The fastest subset of pulsars are known as millisecond pulsars (MSPs), quite simply due to their millisecond spin periods. In addition to their blistering rotations, MSP periods also decay slowly relative to other pulsars due to surface magnetic fields several

orders of magnitudes lower than the general pulsar population. Their extreme characteristics are thought to be attained in a suitably exotic manner; the recycling scenario ascribes the ‘spin-up’ of an old, slow neutron star to the accretion of mass from a binary companion. This transfers angular momentum on to the neutron star, accelerating its spin speed. Given a suitably long period of mass transfer, the neutron star may be spun up to millisecond periods (Alpar et al. 1982; Bhattacharya & van den Heuvel 1991).

<footnote n="1">¹Jodrell Bank Centre for Astrophysics, The University of Manchester, Alan Turing Building, Manchester M13 9PL, UK</footnote>

<footnote n="2">²Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-30167 Hannover, Germany</footnote>

<footnote n="3">³Leibniz Universität Hannover, D-30167 Hannover, Germany</footnote>

<footnote n="4">⁴INAF - Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius (CA), Italy</footnote>

<footnote n="5">⁵Center for Data Intensive and Time Domain Astronomy, Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA</footnote>

<footnote n="6">⁶Department of Physics, National Cheng Kung University, No. 1 University Road, Tainan City 70101, Taiwan</footnote>

<footnote n="7">⁷Max-Planck-Institut für Radioastronomie, D-53121 Bonn, Germany</footnote>

<footnote n="8">⁸South African Radio Astronomy Observatory, 2 Fir Street, Cape Town 7925, South Africa</footnote>

<footnote n="9">⁹Department of Physics and Astronomy, University of Sheffield, Sheffield S3 7RH, UK</footnote>

<footnote n="10">¹⁰Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="11">¹¹Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="12">¹²Center for Research and Exploration in Space Science & Technology II (CRESST II), NASA/GSFC, Greenbelt, MD 20771, USA</footnote>

<footnote n="13">¹³NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="14">¹⁴LPC2E - Université d'Orléans / CNRS, F-45071 Orléans cedex 2, France</footnote>

<footnote n="15">¹⁵Observatoire Radioastronomique de Nançay (ORN), Observatoire de Paris, Université PSL, Univ Orléans, CNRS, F-18330 Nançay, France</footnote>

<footnote n="16">¹⁶School of Physics, Kane Building, University College Cork, Cork T12 K8AF, Ireland</footnote>

<footnote n="17">¹⁷Institute for Defense Analyses, 730 E Glebe Rd, Alexandria, VA 22305, USA</footnote>

<footnote n="18">¹⁸National Astronomical Research Institute of Thailand, Don Kaeo, Mae Rim, Chiang Mai 50180, Thailand</footnote>

<footnote>*E-mail: oli.dodge@gmail.com</footnote>
# Mass estimates from optical modelling of the new TRAPUM redback PSR J1910−5320

O. G. Dodge Ⓓ,<fnref n="1" />★ R. P. Breton Ⓓ,<fnref n="1" /> C. J. Clark Ⓓ,<fnref n="2" />,<fnref n="3" /> M. Burgay Ⓓ,<fnref n="4" /> J. Strader,<fnref n="5" /> K.-Y. Au,<fnref n="6" /> E. D. Barr Ⓓ,<fnref n="7" />

S. Buchner,<fnref n="8" /> V. S. Dhillon Ⓓ,<fnref n="9" />,<fnref n="10" /> E. C. Ferrara,<fnref n="11" />,<fnref n="12" />,<fnref n="13" /> P. C. C. Freire Ⓓ,<fnref n="7" /> J.-M. Griessmeier,<fnref n="14" />,<fnref n="15" />

M. R. Kennedy Ⓓ,<fnref n="1" />,<fnref n="16" /> M. Kramer Ⓓ,<fnref n="1" />,<fnref n="7" /> K.-L. Li,<fnref n="6" /> P. V. Padmanabh Ⓓ,<fnref n="2" />,<fnref n="3" />,<fnref n="7" /> A. Phosrisom,<fnref n="1" /> B. W. Stappers Ⓓ,<fnref n="1" />

S. J. Swihart<fnref n="17" /> and T. Thongmeearkom Ⓓ,<fnref n="1" />,<fnref n="18" />

Accepted 2024 January 17. Received 2024 January 12; in original form 2023 October 17

## ABSTRACT

Spider pulsars continue to provide promising candidates for neutron star mass measurements. Here we present the discovery of PSR J1910−5320, a new millisecond pulsar discovered in a MeerKAT observation of an unidentified Fermi-LAT gamma-ray source. This pulsar is coincident with a recently identified candidate redback binary, independently discovered through its periodic optical flux and radial velocity. New multicolour optical light curves obtained with ULTRACAM/New Technology Telescope in combination with MeerKAT timing and updated SOAR/Goodman spectroscopic radial velocity measurements allow a mass constraint for PSR J1910−5320. ICARUS optical light curve modelling, with streamlined radial velocity fitting, constrains the orbital inclination and companion velocity, unlocking the binary mass function given the precise radio ephemeris. Our modelling aims to unite the photometric and spectroscopic measurements available by fitting each simultaneously to the same underlying physical model, ensuring self-consistency. This targets centre-of-light radial velocity corrections necessitated by the irradiation endemic to spider systems. Depending on the gravity darkening prescription used, we find a moderate neutron star mass of either $1.6 \pm 0.2$ or $1.4 \pm 0.2 M_{\odot}$. The companion mass of either $0.45 \pm 0.04$ or $0.43_{-0.03}^{+0.04} M_{\odot}$ also further confirms PSR J1910−5320 as an irradiated redback spider pulsar.

**Key words:** pulsars: general – pulsars: individual: J1910−5320 – techniques: photometric – techniques: radial velocities.

## 1 INTRODUCTION

The fastest subset of pulsars are known as millisecond pulsars (MSPs), quite simply due to their millisecond spin periods. In addition to their blistering rotations, MSP periods also decay slowly relative to other pulsars due to surface magnetic fields several

orders of magnitudes lower than the general pulsar population. Their extreme characteristics are thought to be attained in a suitably exotic manner; the recycling scenario ascribes the ‘spin-up’ of an old, slow neutron star to the accretion of mass from a binary companion. This transfers angular momentum on to the neutron star, accelerating its spin speed. Given a suitably long period of mass transfer, the neutron star may be spun up to millisecond periods (Alpar et al. 1982; Bhattacharya & van den Heuvel 1991).

<footnote n="1">¹Jodrell Bank Centre for Astrophysics, The University of Manchester, Alan Turing Building, Manchester M13 9PL, UK</footnote>

<footnote n="2">²Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-30167 Hannover, Germany</footnote>

<footnote n="3">³Leibniz Universität Hannover, D-30167 Hannover, Germany</footnote>

<footnote n="4">⁴INAF - Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius (CA), Italy</footnote>

<footnote n="5">⁵Center for Data Intensive and Time Domain Astronomy, Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA</footnote>

<footnote n="6">⁶Department of Physics, National Cheng Kung University, No. 1 University Road, Tainan City 70101, Taiwan</footnote>

<footnote n="7">⁷Max-Planck-Institut für Radioastronomie, D-53121 Bonn, Germany</footnote>

<footnote n="8">⁸South African Radio Astronomy Observatory, 2 Fir Street, Cape Town 7925, South Africa</footnote>

<footnote n="9">⁹Department of Physics and Astronomy, University of Sheffield, Sheffield S3 7RH, UK</footnote>

<footnote n="10">¹⁰Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="11">¹¹Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="12">¹²Center for Research and Exploration in Space Science & Technology II (CRESST II), NASA/GSFC, Greenbelt, MD 20771, USA</footnote>

<footnote n="13">¹³NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="14">¹⁴LPC2E - Université d'Orléans / CNRS, F-45071 Orléans cedex 2, France</footnote>

<footnote n="15">¹⁵Observatoire Radioastronomique de Nançay (ORN), Observatoire de Paris, Université PSL, Univ Orléans, CNRS, F-18330 Nançay, France</footnote>

<footnote n="16">¹⁶School of Physics, Kane Building, University College Cork, Cork T12 K8AF, Ireland</footnote>

<footnote n="17">¹⁷Institute for Defense Analyses, 730 E Glebe Rd, Alexandria, VA 22305, USA</footnote>

<footnote n="18">¹⁸National Astronomical Research Institute of Thailand, Don Kaeo, Mae Rim, Chiang Mai 50180, Thailand</footnote>

<footnote>*E-mail: oli.dodge@gmail.com</footnote>

<footnote>© 2024 The Author(s)</footnote>

<footnote>Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.</footnote>
Given the recycling scenario, spinning up an MSP requires a companion. However, since around 20 per cent of known MSPs are isolated (Jiang et al. 2020), one needs to explore how these seemingly lost their companion. The discovery of the first 'black widow' MSP by Fruchter, Stinebring & Taylor (1988) presented one possible formation mechanism, and established the 'spider' pulsar sub-class of MSPs. Typically a spider system pairs a low-mass, non-degenerate companion with an MSP in a compact ($P_B < 24$ h) orbit. The companion is tidally locked to the pulsar, thus the irradiating pulsar wind heats one face, while the opposite side remains cooler (Djorgovski & Evans 1988). This irradiation ablates material from the companion which often results in eclipsing of the pulsar's beam at radio frequencies (see e.g. Polzin et al. 2020), as well as leading to their nicknames – associating their cannibalistic tendencies with arachnid analogues. Though spiders initially appeared a promising route to isolated MSPs, it still remain highly uncertain as to whether full evaporation within a Hubble time is a realistic option (see e.g. Stappers et al. 1996; Polzin et al. 2020; Kandel, Romani & An 2021). In any case, they provide fascinating environments to study the pulsar wind and high-energy particle physics.

Spider pulsars are typically split into two categories based on their companion mass: black widows with extremely low mass ($M_C < 0.05 M_\odot$) and rebacks with higher companion masses ($M_C \gtrsim 0.1 M_\odot$) (Roberts 2013). Black widows normally have single-peaked light curves over an orbital period, as the impinging irradiation flux dominates the companion star's base temperature. Redbacks light curves can also often exhibit strong irradiation, though unlike black widows it is not ever-present as their base temperatures are higher. Thus, the relative contribution to their light curves of ellipsoidal modulation caused by the tidal distortion of the star is important and produces two peaks per orbital period (see Turchetta et al. 2023, for discussion on the interplay between irradiation and tidal effects in rebacks). Three rebacks, known as transitional millisecond pulsars (tMSPs), were witnessed to switch between MSP (radio-loud) and accreting low-mass X-ray binary states, with each state typically lasting a few years or more. tMSPs are hailed as providing clear evidence for the recycling scenario described above (Archibald et al. 2009; Papitto et al. 2013; Bassa et al. 2014; Stappers et al. 2014).

Constraining the neutron star equation of state (EoS), through neutron star mass measurements (Özel & Freire 2016), fuels a great deal of interest in spider pulsars. Linares (2019) has demonstrated that spiders often host particularly massive neutron stars, with several contending to be the most massive neutron star observed. The original black widow, PSR B1957+20, for a time seemed the heaviest known neutron star, clocking in at 2.4 $M_\odot$ (van Kerkwijk, Breton & Kulkarni 2011). Improved knowledge and data around γ-ray eclipsing in spiders has since revised this measurement down significantly (Clark et al. 2023a), but the promise of massive neutron stars in spider systems remains. There are many EoS model contenders, each predicting a maximum possible neutron star mass. Thus, by observing and measuring massive neutron stars, any EoS predicting a maximum mass below that of the most massive known neutron star can be discarded. The binary nature of spiders where both components can be studied separately therefore provides a convenient avenue to constraining neutron star masses. Radio timing provides the orbital period and pulsar radial velocity, while optical observations can determine inclination and companion radial velocity from photometric and spectroscopic modelling, respectively. Once put together, these can constrain the masses in the system. This then motivates the work in this paper: any new spider to be characterized provides valuable mass

measurements and a potential to constrain the EoS. Whilst for optical modelling one must make certain assumptions regarding the companion heating, producing systematics when compared with other neutron star mass measurements (see Özel & Freire 2016), Romani et al. (2021); Kennedy et al. (2022) and Clark et al. (2021) clearly demonstrate the potential spiders have for precise mass determinations.

Spectroscopic modelling of spiders is relatively novel field, certainly when compared with its photometric counterpart. Both sides of spider modelling are far from complete providing complete descriptions of the companion, with spectroscopic modelling in particular suffering from its extreme computational expense. Aside from technical concerns, the fundamental complications when measuring the radial velocity in spider binaries from observations are summarized as 'centre-of-light' effects. Determining the binary mass ratio requires to combine the well-measured pulsar's projected semimajor axis with a value of the companion's projected centre-of-mass radial velocity. However the radial velocities derived from observed spectroscopy track the centre of light of the particular line or set of lines observed. Indeed, the non-uniform temperature and non-spherical shape of the companion imply that the strength of a line may vary greatly across its surface, which translates into a line velocity that is offset from the centre of mass, therefore producing a different projected radial velocity amplitude but also an orbital profile which may depart slightly from the perfect function expected from a circular orbit.

Several approaches have been used to connect the observed radial velocities to the correct centre-of-mass radial velocity amplitude. van Kerkwijk, Breton & Kulkarni (2011) and Romani et al. (2021) both produced synthetic radial velocity curves which are then fitted to the observed curve to estimate the correction factor. Linares, Shahbaz & Casares (2018), on the other hand, takes a more empirical approach in which observed line species are assessed to originate from the hotter dayside or colder nightside of the companion based on the temperature at which they are produced. In this way, they can 'bracket' line velocities to lie between the true centre of mass and the maximal extent of the star in either direction. Finally, Kennedy et al. (2022) implemented the ultimate step in producing full synthetic spectra which are directly fitted to the raw observed spectroscopy. This modelling of the photometry and spectroscopy ensures the necessary centre-of-light corrections are intrinsically embedded in the line profile which is self-consistent with the heating model at any given parameters.

Follow-up observations are fruitful in various wavelengths; Ray et al. (2013) reported the discovery of 43 new MSPs, many of which were spiders, from the first generation of deep radio searches targeting unassociated Fermi-LAT sources. The population has kept growing since, with the latest Fermi-LAT survey reporting at least 110 MSPs discovered in this fashion (Smith et al. 2023). In addition to these, Clark et al. (2023b) detailed a new MeerKAT L-band survey of LAT sources in which nine new MSPs were found among 79 Fermi-LAT sources, including two new rebacks. Optical searching of similar fields, with or without prior radio search, can also produce new spider candidates by looking for the signature orbital modulation in the light curves described earlier, with spectroscopy possibly providing further evidence through the system's mass function (see e.g. Strader et al. 2015, 2016; Swihart et al. 2022).

One such recent discovery is that of a candidate redback binary system within the previously unidentified gamma-ray source 4FGL J1910.7–5320 (Au et al. 2023). The discovery is a fruit of cross-matching the 4FGL-DR3 catalogue against sub-24 h period optical variables in Catalina Real-Time Transient Surveys (Drake
Given the recycling scenario, spinning up an MSP requires a companion. However, since around 20 per cent of known MSPs are isolated (Jiang et al. 2020), one needs to explore how these seemingly lost their companion. The discovery of the first ‘black widow’ MSP by Fruchter, Stinebring & Taylor (1988) presented one possible formation mechanism, and established the ‘spider’ pulsar sub-class of MSPs. Typically a spider system pairs a low-mass, non-degenerate companion with an MSP in a compact ($P_B < 24$ h) orbit. The companion is tidally locked to the pulsar, thus the irradiating pulsar wind heats one face, while the opposite side remains cooler (Djorgovski & Evans 1988). This irradiation ablates material from the companion which often results in eclipsing of the pulsar’s beam at radio frequencies (see e.g. Polzin et al. 2020), as well as leading to their nicknames – associating their cannibalistic tendencies with arachnid analogues. Though spiders initially appeared a promising route to isolated MSPs, it still remain highly uncertain as to whether full evaporation within a Hubble time is a realistic option (see e.g. Stappers et al. 1996; Polzin et al. 2020; Kandel, Romani & An 2021). In any case, they provide fascinating environments to study the pulsar wind and high-energy particle physics.

Spider pulsars are typically split into two categories based on their companion mass: black widows with extremely low mass ($M_C < 0.05 M_\odot$) and redbacks with higher companion masses ($M_C \gtrsim 0.1 M_\odot$) (Roberts 2013). Black widows normally have single-peaked light curves over an orbital period, as the impinging irradiation flux dominates the companion star’s base temperature. Redbacks light curves can also often exhibit strong irradiation, though unlike black widows it is not ever-present as their base temperatures are higher. Thus, the relative contribution to their light curves of ellipsoidal modulation caused by the tidal distortion of the star is important and produces two peaks per orbital period (see Turchetta et al. 2023, for discussion on the interplay between irradiation and tidal effects in redbacks). Three redbacks, known as transitional millisecond pulsars (tMSPs), were witnessed to switch between MSP (radio-loud) and accreting low-mass X-ray binary states, with each state typically lasting a few years or more. tMSPs are hailed as providing clear evidence for the recycling scenario described above (Archibald et al. 2009; Papitto et al. 2013; Bassa et al. 2014; Stappers et al. 2014).

Constraining the neutron star equation of state (EoS), through neutron star mass measurements (Özel & Freire 2016), fuels a great deal of interest in spider pulsars. Linares (2019) has demonstrated that spiders often host particularly massive neutron stars, with several contending to be the most massive neutron star observed. The original black widow, PSR B1957+20, for a time seemed the heaviest known neutron star, clocking in at 2.4 $M_\odot$ (van Kerkwijk, Breton & Kulkarni 2011). Improved knowledge and data around γ-ray eclipsing in spiders has since revised this measurement down significantly (Clark et al. 2023a), but the promise of massive neutron stars in spider systems remains. There are many EoS model contenders, each predicting a maximum possible neutron star mass. Thus, by observing and measuring massive neutron stars, any EoS predicting a maximum mass below that of the most massive known neutron star can be discarded. The binary nature of spiders where both components can be studied separately therefore provides a convenient avenue to constraining neutron star masses. Radio timing provides the orbital period and pulsar radial velocity, while optical observations can determine inclination and companion radial velocity from photometric and spectroscopic modelling, respectively. Once put together, these can constrain the masses in the system. This then motivates the work in this paper: any new spider to be characterized provides valuable mass

measurements and a potential to constrain the EoS. Whilst for optical modelling one must make certain assumptions regarding the companion heating, producing systematics when compared with other neutron star mass measurements (see Özel & Freire 2016), Romani et al. (2021); Kennedy et al. (2022) and Clark et al. (2021) clearly demonstrate the potential spiders have for precise mass determinations.

Spectroscopic modelling of spiders is relatively novel field, certainly when compared with its photometric counterpart. Both sides of spider modelling are far from complete providing complete descriptions of the companion, with spectroscopic modelling in particular suffering from its extreme computational expense. Aside from technical concerns, the fundamental complications when measuring the radial velocity in spider binaries from observations are summarized as ‘centre-of-light’ effects. Determining the binary mass ratio requires to combine the well-measured pulsar’s projected semimajor axis with a value of the companion’s projected centre-of-mass radial velocity. However the radial velocities derived from observed spectroscopy track the centre of light of the particular line or set of lines observed. Indeed, the non-uniform temperature and non-spherical shape of the companion imply that the strength of a line may vary greatly across its surface, which translates into a line velocity that is offset from the centre of mass, therefore producing a different projected radial velocity amplitude but also an orbital profile which may depart slightly from the perfect function expected from a circular orbit.

Several approaches have been used to connect the observed radial velocities to the correct centre-of-mass radial velocity amplitude. van Kerkwijk, Breton & Kulkarni (2011) and Romani et al. (2021) both produced synthetic radial velocity curves which are then fitted to the observed curve to estimate the correction factor. Linares, Shahbaz & Casares (2018), on the other hand, takes a more empirical approach in which observed line species are assessed to originate from the hotter dayside or colder nightside of the companion based on the temperature at which they are produced. In this way, they can ‘bracket’ line velocities to lie between the true centre of mass and the maximal extent of the star in either direction. Finally, Kennedy et al. (2022) implemented the ultimate step in producing full synthetic spectra which are directly fitted to the raw observed spectroscopy. This modelling of the photometry and spectroscopy ensures the necessary centre-of-light corrections are intrinsically embedded in the line profile which is self-consistent with the heating model at any given parameters.

Follow-up observations are fruitful in various wavelengths; Ray et al. (2013) reported the discovery of 43 new MSPs, many of which were spiders, from the first generation of deep radio searches targeting unassociated Fermi-LAT sources. The population has kept growing since, with the latest Fermi-LAT survey reporting at least 110 MSPs discovered in this fashion (Smith et al. 2023). In addition to these, Clark et al. (2023b) detailed a new MeerKAT L-band survey of LAT sources in which nine new MSPs were found among 79 Fermi-LAT sources, including two new redbacks. Optical searching of similar fields, with or without prior radio search, can also produce new spider candidates by looking for the signature orbital modulation in the light curves described earlier, with spectroscopy possibly providing further evidence through the system’s mass function (see e.g. Strader et al. 2015, 2016; Swihart et al. 2022).

One such recent discovery is that of a candidate redback binary system within the previously unidentified gamma-ray source 4FGL J1910.7–5320 (Au et al. 2023). The discovery is a fruit of cross-matching the 4FGL-DR3 catalogue against sub-24 h period optical variables in Catalina Real-Time Transient Surveys (Drake
et al. 2017). 4FGL J1910.7−5320 was one of two spiders found in this way (the other being PSR J0955−3947; Li et al. 2018). SOAR/Goodman spectroscopy was also obtained, from which a sinusoidal radial velocity curve confirmed the binary nature of the system with an orbital period $P_B = 0.34847592$ d. The observed radial velocity amplitude, $K_{2,obs} = 218 \pm 8$ km s$^{-1}$, is in line with what is seen in many redback systems, thus favoured as a redback candidate. Independently of this optical discovery, we detected radio pulsations from this source as part of an ongoing survey for new pulsars in *Fermi*-LAT sources (Clark et al. 2023b) being performed as part of the TRAnsients and Pulsars Using MeerKAT (TRAPUM) large survey project (Stappers & Kramer 2016). This confirmed the redback prediction of Au et al. (2023).

In this paper, we present the TRAPUM discovery of radio pulsations from the neutron star associated with 4FGL J1910.7−5320 using the MeerKAT telescope. In Section 2 we describe the radio discovery and timing of the new pulsar, PSR J1910−5320, as well as multiband optical photometry obtained with ULTRACAM on the ESO New Technology Telescope (NTT). Section 3 details the optical modelling of the optical light curves. In particular, we introduce a novel method to utilize values provided by radial velocity measurements made from optical spectroscopy. This modelling provides constraints on component masses, through the inclination and companion velocity, further confirming J1910 as a redback. Section 4 discusses the physical interpretation of our modelling, including an analysis of the impact of different gravity darkening prescriptions on the final results and an assessment of centre-of-light location where the absorption features are produced. A summary and conclusion is provided in Section 5.

## 2 OBSERVATIONS

### 2.1 Radio discovery and timing

In Clark et al. (2023b), we presented the first results from an ongoing survey being performed as part of the TRAPUM large survey project (Stappers & Kramer 2016) using the MeerKAT radio telescope (Jonas 2009; Jonas & the MeerKAT Team 2016) to search for new pulsars in unassociated pulsar-like *Fermi*-LAT sources. The survey presented therein consisted of two 10-min observations of 79 sources from the 4FGL catalogue (Abdollahi et al. 2020), conducted using MeerKAT’s L-band receiver (at observing frequencies in the range 856–1712 MHz). This project has since been extended with a further two-pass survey (Thongmeearkom et al. in preparation) being performed with the UHF receiver (544–1088 MHz). Tied-array beams cover a larger solid angle at this lower frequency band, and so a small number of additional *Fermi*-LAT sources whose localization regions were too uncertain to cover in single observations at L band were added to this UHF survey. One of these new sources was 4FGL J1910.7−5320.

TRAPUM observed this source on 2022 May 31, and detected highly significant radio pulsations with signal-to-noise ratio, S/N ≈380. The signal had a spin period of 2.33 ms and significant acceleration of $4.12 \pm 0.02$ m s$^{-2}$ indicative of an MSP in a short-period binary system. We used *SEEKAT*<fnref n="1" /> (Bezuidenhout et al. 2023) to localize this signal to a position less than 0.5 arcsec from an optical star detected in the *Gaia* DR3 (Gaia Collaboration 2016, 2023) and Catalina Surveys Southern (CSS) periodic variable star catalogues (Drake et al. 2017). The CSS catalogue lists this source

as having a 16.8 h periodicity, with a double-peaked light curve of 1.1 mag amplitude. However, such a light curve is inconsistent with that of a pulsar binary companion, as the ellipsoidal modulation that gives rise to a double-peaked light curve has a maximum amplitude of around 0.3 mag. However, folding the CSS data with half this period leaves a single-peaked light curve that is consistent with an irradiated binary pulsar companion star. Unknown to us at the time, this 8.4 h orbital period was independently confirmed by the optical spectroscopy presented in Au et al. (2023) through the measurement of Doppler-shifted spectral.

We therefore proceeded under the assumption that this star was indeed an irradiated redback counterpart to our newly detected MSP, and used the CSS ephemeris to schedule follow-up timing observations with both MeerKAT and Murriyang, the Parkes 64 m telescope, during the half of the orbit centred on the companion star’s superior conjunction (i.e. orbital phases between 0.5 and 1.0) when the pulsar should not be eclipsed by wind from the companion.

Our timing campaign with MeerKAT consisted of 15 pseudo-logarithmically spaced observations between 2022 June 29 and 2022 September 29 with several observations on the first days (2022 June 29 and 2022 June 30) and increasing intervals between subsequent observations to facilitate phase connection. These observations each lasted 5 min, and were taken using the Pulsar Timing User Supplied Instrument (PTUSE, Bailes et al. 2020) with coherent de-dispersion. The first eight observations were taken with MeerKAT’s UHF receiver, the rest were performed at L band. A second pseudo-logarithmic timing campaign began with Parkes on 2022 September 6 until 2023 March 25. These observations each lasted 1.5 h using the Ultra-wide-band Low receiver (Hobbs et al. 2020), covering a frequency range from 0.7 to 4 GHz, with coherent de-dispersion. The resulting data were reduced using standard radio timing techniques, as described by Clark et al. (2023b); additional details will be provided elsewhere.

The resultant pulse times of arrival at the location of the radio telescope (ToAs) were analysed using the TEMPO (Nice et al. 2015) timing package. To model the motion of the radio telescope relative to the Solar system barycentre, we used the Jet Propulsion laboratory’s DE421 Solar system ephemeris (Folkner, Williams & Boggs 2009). To model the pulsar’s orbit, we used the BTX orbital model, which allows for the measurement of multiple orbital frequency derivatives. This is necessary because, like in most other redback systems, the ToAs revealed unpredictable deviations in the times of the pulsar’s ascending node of the order of a few seconds, thought to be due to orbital period variations caused by variability of the companion star’s gravitational quadrupole moment via the Applegate mechanism (Applegate 1992). The parameters of the timing solution are presented in Table 1, where the numbers in parentheses indicated the $1\sigma$ uncertainties on the last digits of the nominal values. These parameters are presented in the Dynamic Barycentric time (TDB).

The determination of the timing solution was greatly assisted by previous knowledge of the orbital period (from CSS photometry) and the *Gaia* astrometry, which was assumed for this solution.

### 2.2 Optical photometry

We obtained multiband light curves of J1910 on two nights, 2022 June 28 and 30, using the ULTRACAM high-speed multiband photometer (Dhillon et al. 2007), mounted on the 3.50 m NTT at the European Southern Observatory (ESO) La Silla, Chile. These observations are shown in Fig. 1, with the corresponding times and length provided in Table 2. ULTRACAM utilizes three CCDs simultaneously, each using a different Super Sloan Digital Sky

<footnote n="1">¹https://github.com/BezuidenhoutMC/SeeKAT</footnote>
et al. 2017). 4FGL J1910.7−5320 was one of two spiders found in this way (the other being PSR J0955−3947; Li et al. 2018). SOAR/Goodman spectroscopy was also obtained, from which a sinusoidal radial velocity curve confirmed the binary nature of the system with an orbital period $P_B = 0.34847592$ d. The observed radial velocity amplitude, $K_{2,obs} = 218 \pm 8$ km s⁻<fnref n="1" />, is in line with what is seen in many redback systems, thus favoured as a redback candidate. Independently of this optical discovery, we detected radio pulsations from this source as part of an ongoing survey for new pulsars in *Fermi*–LAT sources (Clark et al. 2023b) being performed as part of the TRAnsients and Pulsars Using MeerKAT (TRAPUM) large survey project (Stappers & Kramer 2016). This confirmed the redback prediction of Au et al. (2023).

In this paper, we present the TRAPUM discovery of radio pulsations from the neutron star associated with 4FGL J1910.7−5320 using the MeerKAT telescope. In Section 2 we describe the radio discovery and timing of the new pulsar, PSR J1910−5320, as well as multiband optical photometry obtained with ULTRACAM on the ESO New Technology Telescope (NTT). Section 3 details the optical modelling of the optical light curves. In particular, we introduce a novel method to utilize values provided by radial velocity measurements made from optical spectroscopy. This modelling provides constraints on component masses, through the inclination and companion velocity, further confirming J1910 as a redback. Section 4 discusses the physical interpretation of our modelling, including an analysis of the impact of different gravity darkening prescriptions on the final results and an assessment of centre-of-light location where the absorption features are produced. A summary and conclusion is provided in Section 5.

## 2 OBSERVATIONS

### 2.1 Radio discovery and timing

In Clark et al. (2023b), we presented the first results from an ongoing survey being performed as part of the TRAPUM large survey project (Stappers & Kramer 2016) using the MeerKAT radio telescope (Jonas 2009; Jonas & the MeerKAT Team 2016) to search for new pulsars in unassociated pulsar-like *Fermi*–LAT sources. The survey presented therein consisted of two 10-min observations of 79 sources from the 4FGL catalogue (Abdollahi et al. 2020), conducted using MeerKAT’s L-band receiver (at observing frequencies in the range 856–1712 MHz). This project has since been extended with a further two-pass survey (Thongmeearkom et al. in preparation) being performed with the UHF receiver (544–1088 MHz). Tied-array beams cover a larger solid angle at this lower frequency band, and so a small number of additional *Fermi*–LAT sources whose localization regions were too uncertain to cover in single observations at L band were added to this UHF survey. One of these new sources was 4FGL J1910.7−5320.

TRAPUM observed this source on 2022 May 31, and detected highly significant radio pulsations with signal-to-noise ratio, S/N ≈380. The signal had a spin period of 2.33 ms and significant acceleration of $4.12 \pm 0.02$ m s⁻² indicative of an MSP in a short-period binary system. We used *SEEKAT*<fnref n="1" /> (Bezuidenhout et al. 2023) to localize this signal to a position less than 0.5 arcsec from an optical star detected in the *Gaia* DR3 (Gaia Collaboration 2016, 2023) and Catalina Surveys Southern (CSS) periodic variable star catalogues (Drake et al. 2017). The CSS catalogue lists this source

as having a 16.8 h periodicity, with a double-peaked light curve of 1.1 mag amplitude. However, such a light curve is inconsistent with that of a pulsar binary companion, as the ellipsoidal modulation that gives rise to a double-peaked light curve has a maximum amplitude of around 0.3 mag. However, folding the CSS data with half this period leaves a single-peaked light curve that is consistent with an irradiated binary pulsar companion star. Unknown to us at the time, this 8.4 h orbital period was independently confirmed by the optical spectroscopy presented in Au et al. (2023) through the measurement of Doppler-shifted spectral.

We therefore proceeded under the assumption that this star was indeed an irradiated redback counterpart to our newly detected MSP, and used the CSS ephemeris to schedule follow-up timing observations with both MeerKAT and Murriyang, the Parkes 64 m telescope, during the half of the orbit centred on the companion star’s superior conjunction (i.e. orbital phases between 0.5 and 1.0) when the pulsar should not be eclipsed by wind from the companion.

Our timing campaign with MeerKAT consisted of 15 pseudo-logarithmically spaced observations between 2022 June 29 and 2022 September 29 with several observations on the first days (2022 June 29 and 2022 June 30) and increasing intervals between subsequent observations to facilitate phase connection. These observations each lasted 5 min, and were taken using the Pulsar Timing User Supplied Instrument (PTUSE, Bailes et al. 2020) with coherent de-dispersion. The first eight observations were taken with MeerKAT’s UHF receiver, the rest were performed at L band. A second pseudo-logarithmic timing campaign began with Parkes on 2022 September 6 until 2023 March 25. These observations each lasted 1.5 h using the Ultra-wide-band Low receiver (Hobbs et al. 2020), covering a frequency range from 0.7 to 4 GHz, with coherent de-dispersion. The resulting data were reduced using standard radio timing techniques, as described by Clark et al. (2023b); additional details will be provided elsewhere.

The resultant pulse times of arrival at the location of the radio telescope (ToAs) were analysed using the TEMPO (Nice et al. 2015) timing package. To model the motion of the radio telescope relative to the Solar system barycentre, we used the Jet Propulsion laboratory’s DE421 Solar system ephemeris (Folkner, Williams & Boggs 2009). To model the pulsar’s orbit, we used the BTX orbital model, which allows for the measurement of multiple orbital frequency derivatives. This is necessary because, like in most other redback systems, the ToAs revealed unpredictable deviations in the times of the pulsar’s ascending node of the order of a few seconds, thought to be due to orbital period variations caused by variability of the companion star’s gravitational quadrupole moment via the Applegate mechanism (Applegate 1992). The parameters of the timing solution are presented in Table 1, where the numbers in parentheses indicated the $1\sigma$ uncertainties on the last digits of the nominal values. These parameters are presented in the Dynamic Barycentric time (TDB).

The determination of the timing solution was greatly assisted by previous knowledge of the orbital period (from CSS photometry) and the *Gaia* astrometry, which was assumed for this solution.

### 2.2 Optical photometry

We obtained multiband light curves of J1910 on two nights, 2022 June 28 and 30, using the ULTRACAM high-speed multiband photometer (Dhillon et al. 2007), mounted on the 3.50 m NTT at the European Southern Observatory (ESO) La Silla, Chile. These observations are shown in Fig. 1, with the corresponding times and length provided in Table 2. ULTRACAM utilizes three CCDs simultaneously, each using a different Super Sloan Digital Sky

<footnote n="1">¹https://github.com/BezuidenhoutMC/SeeKAT</footnote>
<table><caption>**Table 1.** Timing solution for PSR J1910−5320, obtained using the BTX orbital model.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td colspan="2"><em>Gaia</em> DR3 astrometry</td></tr><tr><td>R.A., α (J2000)</td><td>19<sup>h</sup>10<sup>m</sup>49<sup>s</sup>12053(1)</td></tr><tr><td>Dec., δ (J2000)</td><td>−53°20′57″1205(2)</td></tr><tr><td>Proper motion in α, μαcosδ (mas yr<sup>−1</sup>)</td><td>1.7 ± 0.2</td></tr><tr><td>Proper motion in δ, μδ (mas yr<sup>−1</sup>)</td><td>−6.8 ± 0.2</td></tr><tr><td>Parallax, ω (mas)</td><td>−0.42 ± 0.26</td></tr><tr><td>Epoch of position measurement (MJD)</td><td>57388.0</td></tr><tr><td colspan="2"><em>Timing parameters</em></td></tr><tr><td>Solar-system ephemeris</td><td>DE421</td></tr><tr><td>Time-scale</td><td>TDB</td></tr><tr><td>Data span (MJD)</td><td>59759.8−59978.5</td></tr><tr><td>Epoch of spin period measurement (MJD)</td><td>59760</td></tr><tr><td>Number of ToAs</td><td>939</td></tr><tr><td>Residual rms (μs)</td><td>5.07</td></tr><tr><td>Reduced χ<sup>2</sup></td><td>2.1</td></tr><tr><td>Spin frequency, ν (Hz)</td><td>428.7490184657(3)</td></tr><tr><td>Spin-down rate, ṽ (Hz s<sup>−1</sup>)</td><td>−6.80(7) × 10<sup>−15</sup></td></tr><tr><td>Dispersion measure, DM (pc cm<sup>−3</sup>)</td><td>24.42</td></tr><tr><td>Binary model</td><td>BTX</td></tr><tr><td>Orbital frequency, ν<sub>orb</sub> (Hz)</td><td>3.32132606(1) × 10<sup>−5</sup></td></tr><tr><td>First orbital frequency derivative, ṽ<sub>orb</sub> (Hz s<sup>−1</sup>)</td><td>−3.45(4) × 10<sup>−18</sup></td></tr><tr><td>Second orbital frequency derivative, ṽ<sub>orb</sub> (Hz s<sup>−2</sup>)</td><td>1.85(4) × 10<sup>−25</sup></td></tr><tr><td>Projected semimajor axis, x (lt s)</td><td>0.969183(6)</td></tr><tr><td>Epoch of ascending node, T<sub>asc</sub> (MJD)</td><td>59759.9208124(3)</td></tr><tr><td colspan="2"><em>Derived parameters</em></td></tr><tr><td>Spin period, p (ms)</td><td>2.33236685551(2)</td></tr><tr><td>Spin period derivative, ṝ</td><td>3.70(4) × 10<sup>−20</sup></td></tr><tr><td>Orbital period, P<sub>orb</sub> (d)</td><td>0.348477501(1)</td></tr><tr><td>Spin-down power, Ẽ (erg/s)</td><td>1.15 × 10<sup>35</sup></td></tr><tr><td>Surface magnetic field strength, B<sub>S</sub> (G)</td><td>3.0 × 10<sup>8</sup></td></tr><tr><td>Light-cylinder magnetic field strength, B<sub>LC</sub> (G)</td><td>2.2 × 10<sup>5</sup></td></tr></tbody></table>

Survey (Super-SDSS) $u_s g_s r_s i_s z_s$ filter (Dhillon et al. 2021). For these observations CCDs 1, 2, and 3 used the $r_s$, $g_s$, and $u_s$ filters, respectively. The data were taken under photometric conditions, with seeing varying between 1 and 1.5 arcsec. The observations were reduced using the HiPERCAM (Dhillon et al. 2016, 2018) pipeline.<fnref n="2" /> Ensemble photometry (Honeycutt 1992) was used to calibrate the $r_s$ and $g_s$ bands. 12 nearby stars with known Gaia magnitudes were chosen as reference apertures. In order to use the Gaia magnitudes, they were transformed first into the SDSS prime r' and g' bands, then again into the corresponding HiPERCAM filters (appendix A Brown et al. 2022). Due to a lack of Gaia transform, and the unreliable transform between the HiPERCAM and SDSS filters, the u' band was calibrated using the instrumental zero-point determined by observing the known SDSS standard PG1323-086D. After processing the data we were left with 3746 data points: 1608 and 1291 from the $r_s$ and $g_s$ bands, respectively (20 s exposures), and 530 from the $u_s$ (60 s exposures). Co-addition of $u_s$-band exposures, maximizing S/N, leaves fewer $u_s$ data points relative to the other bands. The orbital phase of each point was calculated using the ephemeris given in Table 1. Here the light curve phases have been folded as assumed in our ephemeris, with $\phi = 0$ corresponding to the ascending node of the pulsar. Phases 0.25 and 0.75 therefore correspond with the companion's inferior and superior conjunctions, respectively.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Phased light curves of J1910 in the $r_s$, $g_s$, and $u_s$ bands. The phasing is calculated according to the solution presented in Table 1, with the companion's inferior and superior conjunctions occurring at phases 0.25 and 0.75, respectively. Two repeated cycles are shown for clarity. The gap in coverage in all bands around $\phi = 0.90$ is due to poor focus during the initial stages of the 2022 June 28 observing run. A larger portion of the $g_s$-band light curve is excluded due to irreducible artifacts in the data.</figcaption></figure>

## 2.3 SOAR/Goodman spectroscopy

The SOAR/Goodman spectroscopic data set for PSR J1910−5320 is identical to that described in Au et al. (2023). However, we found that the orbital ephemerides inferred from these data show relatively modest but nevertheless quite statistically significant discrepancies with the ephemerides derived from pulsar timing. An investigation of these discrepancies led to the conclusion that a greater than expected degree of flexure was present in the previous SOAR/Goodman observations. Despite having calibration arc lamp observations continually interspersed throughout the object observations, and using night sky lines for an additional zero-point correction, some residual effects of flexure remained. This could perhaps be associated with spatial flexure somewhere along the light path in the instrument, or instead with imperfect guiding that led to miscentring of the source in the slit.

Therefore, we have re-derived the PSR J1910−5320 radial velocities through a process that differs in some details from the method used in Au et al. (2023). To improve the wavelength zero-point corrections, we use the TELFIT code (Gullikson, Dodson-Robinson & Kraus 2014) to generate a telluric absorption spectrum based on the airmass, the local humidity, pressure, and temperature, and the 3-h Global Data Assimilation System atmospheric model closest in time to each object spectrum. This model spectrum, smoothed to the resolution of the SOAR data and binned to the same pixel scale, is then fit to the object spectrum in the region of the Fraunhofer A band (7580–7700 Å) to determine the wavelength zero-point correction. While other telluric features are also present in some spectra, this is the only telluric feature measurable in essentially all usable spectra, even those of low signal-to-noise, so we restrict the fit to this feature. Comparisons over a number of data sets show that the

<footnote n="2">²http://deneb.astro.warwick.ac.uk/phsaap/hipercam/docs/html/</footnote>
<table><caption>**Table 1.** Timing solution for PSR J1910−5320, obtained using the BTX orbital model.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td colspan="2"><em>Gaia</em> DR3 astrometry</td></tr><tr><td>R.A., α (J2000)</td><td>19<sup>h</sup>10<sup>m</sup>49<sup>s</sup>12053(1)</td></tr><tr><td>Dec., δ (J2000)</td><td>−53°20′57″1205(2)</td></tr><tr><td>Proper motion in α, μαcosδ (mas yr<sup>−1</sup>)</td><td>1.7 ± 0.2</td></tr><tr><td>Proper motion in δ, μδ (mas yr<sup>−1</sup>)</td><td>−6.8 ± 0.2</td></tr><tr><td>Parallax, ω (mas)</td><td>−0.42 ± 0.26</td></tr><tr><td>Epoch of position measurement (MJD)</td><td>57388.0</td></tr><tr><td colspan="2"><em>Timing parameters</em></td></tr><tr><td>Solar-system ephemeris</td><td>DE421</td></tr><tr><td>Time-scale</td><td>TDB</td></tr><tr><td>Data span (MJD)</td><td>59759.8−59978.5</td></tr><tr><td>Epoch of spin period measurement (MJD)</td><td>59760</td></tr><tr><td>Number of ToAs</td><td>939</td></tr><tr><td>Residual rms (μs)</td><td>5.07</td></tr><tr><td>Reduced χ<sup>2</sup></td><td>2.1</td></tr><tr><td>Spin frequency, ν (Hz)</td><td>428.7490184657(3)</td></tr><tr><td>Spin-down rate, ṷ (Hz s<sup>−1</sup>)</td><td>−6.80(7) × 10<sup>−15</sup></td></tr><tr><td>Dispersion measure, DM (pc cm<sup>−3</sup>)</td><td>24.42</td></tr><tr><td>Binary model</td><td>BTX</td></tr><tr><td>Orbital frequency, ν<sub>orb</sub> (Hz)</td><td>3.32132606(1) × 10<sup>−5</sup></td></tr><tr><td>First orbital frequency derivative, ṽ<sub>orb</sub> (Hz s<sup>−1</sup>)</td><td>−3.45(4) × 10<sup>−18</sup></td></tr><tr><td>Second orbital frequency derivative, ṽ<sub>orb</sub> (Hz s<sup>−2</sup>)</td><td>1.85(4) × 10<sup>−25</sup></td></tr><tr><td>Projected semimajor axis, x (lt s)</td><td>0.969183(6)</td></tr><tr><td>Epoch of ascending node, T<sub>asc</sub> (MJD)</td><td>59759.9208124(3)</td></tr><tr><td colspan="2"><em>Derived parameters</em></td></tr><tr><td>Spin period, p (ms)</td><td>2.33236685551(2)</td></tr><tr><td>Spin period derivative, ṝ</td><td>3.70(4) × 10<sup>−20</sup></td></tr><tr><td>Orbital period, P<sub>orb</sub> (d)</td><td>0.348477501(1)</td></tr><tr><td>Spin-down power, Ẽ (erg/s)</td><td>1.15 × 10<sup>35</sup></td></tr><tr><td>Surface magnetic field strength, B<sub>S</sub> (G)</td><td>3.0 × 10<sup>8</sup></td></tr><tr><td>Light-cylinder magnetic field strength, B<sub>LC</sub> (G)</td><td>2.2 × 10<sup>5</sup></td></tr></tbody></table>

Survey (Super-SDSS) $u_s g_s r_s i_s z_s$ filter (Dhillon et al. 2021). For these observations CCDs 1, 2, and 3 used the $r_s$, $g_s$, and $u_s$ filters, respectively. The data were taken under photometric conditions, with seeing varying between 1 and 1.5 arcsec. The observations were reduced using the HiPERCAM (Dhillon et al. 2016, 2018) pipeline.<fnref n="2" /> Ensemble photometry (Honeycutt 1992) was used to calibrate the $r_s$ and $g_s$ bands. 12 nearby stars with known Gaia magnitudes were chosen as reference apertures. In order to use the Gaia magnitudes, they were transformed first into the SDSS prime r' and g' bands, then again into the corresponding HiPERCAM filters (appendix A Brown et al. 2022). Due to a lack of Gaia transform, and the unreliable transform between the HiPERCAM and SDSS filters, the u' band was calibrated using the instrumental zero-point determined by observing the known SDSS standard PG1323-086D. After processing the data we were left with 3746 data points: 1608 and 1291 from the $r_s$ and $g_s$ bands, respectively (20 s exposures), and 530 from the $u_s$ (60 s exposures). Co-addition of $u_s$-band exposures, maximizing S/N, leaves fewer $u_s$ data points relative to the other bands. The orbital phase of each point was calculated using the ephemeris given in Table 1. Here the light curve phases have been folded as assumed in our ephemeris, with $\phi = 0$ corresponding to the ascending node of the pulsar. Phases 0.25 and 0.75 therefore correspond with the companion's inferior and superior conjunctions, respectively.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Phased light curves of J1910 in the $r_s$, $g_s$, and $u_s$ bands. The phasing is calculated according to the solution presented in Table 1, with the companion's inferior and superior conjunctions occurring at phases 0.25 and 0.75, respectively. Two repeated cycles are shown for clarity. The gap in coverage in all bands around $\phi = 0.90$ is due to poor focus during the initial stages of the 2022 June 28 observing run. A larger portion of the $g_s$-band light curve is excluded due to irreducible artifacts in the data.</figcaption></figure>

## 2.3 SOAR/Goodman spectroscopy

The SOAR/Goodman spectroscopic data set for PSR J1910−5320 is identical to that described in Au et al. (2023). However, we found that the orbital ephemerides inferred from these data show relatively modest but nevertheless quite statistically significant discrepancies with the ephemerides derived from pulsar timing. An investigation of these discrepancies led to the conclusion that a greater than expected degree of flexure was present in the previous SOAR/Goodman observations. Despite having calibration arc lamp observations continually interspersed throughout the object observations, and using night sky lines for an additional zero-point correction, some residual effects of flexure remained. This could perhaps be associated with spatial flexure somewhere along the light path in the instrument, or instead with imperfect guiding that led to miscentring of the source in the slit.

Therefore, we have re-derived the PSR J1910−5320 radial velocities through a process that differs in some details from the method used in Au et al. (2023). To improve the wavelength zero-point corrections, we use the TELFIT code (Gullikson, Dodson-Robinson & Kraus 2014) to generate a telluric absorption spectrum based on the airmass, the local humidity, pressure, and temperature, and the 3-h Global Data Assimilation System atmospheric model closest in time to each object spectrum. This model spectrum, smoothed to the resolution of the SOAR data and binned to the same pixel scale, is then fit to the object spectrum in the region of the Fraunhofer A band (7580–7700 Å) to determine the wavelength zero-point correction. While other telluric features are also present in some spectra, this is the only telluric feature measurable in essentially all usable spectra, even those of low signal-to-noise, so we restrict the fit to this feature. Comparisons over a number of data sets show that the

<footnote n="2">²http://deneb.astro.warwick.ac.uk/phsaap/hipercam/docs/html/</footnote>
<table><caption>**Table 2.** Time and phase coverage for ULTRACAM photometry obtained of J1910. The phase coverage, calculated with the timing ephemeris provided in Table 1, corresponds with the phase axis of Fig. 1. The *g*<sub>s</sub> filter is split from the other two due to the exclusion of irreducible data for the 2022 June 28 night.</caption><thead><tr><th>Start time (UTC)</th><th>Observation length (h)</th><th colspan="2">Phase coverage</th></tr><tr><th></th><th></th><th>r<sub>s</sub>, u<sub>s</sub></th><th>g<sub>s</sub></th></tr></thead><tbody><tr><td>2022 June 29 (04:11:44)</td><td>4.25</td><td>0.90–1.38</td><td>1.10–1.38</td></tr><tr><td>2022 July 1 (01:31:58)</td><td>5.0</td><td>0.30–0.89</td><td>0.30–0.89</td></tr></tbody></table>

<table><caption>**Table 3.** Updated radial velocities (RV) of PSR J1910-5320 from SOAR for both the full spectrum and targeting just the Mg β triplet.</caption><thead><tr><th rowspan="2">BJD (d)</th><th colspan="2">Full spectrum</th><th colspan="2">Mg β</th></tr><tr><th>RV</th><th>ΔRV (km s<sup>-1</sup>)</th><th>RV</th><th>ΔRV (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>59679.32228</td><td>-17.8</td><td>21.7</td><td>-28.9</td><td>25.3</td></tr><tr><td>59679.34012</td><td>-58.4</td><td>20.6</td><td>-72.4</td><td>26.5</td></tr><tr><td>59679.35806</td><td>-141.8</td><td>18</td><td>-164.2</td><td>26.3</td></tr><tr><td>59680.32745</td><td>80.7</td><td>17.5</td><td>82.1</td><td>21.1</td></tr><tr><td>59680.34496</td><td>67.5</td><td>19.1</td><td>79</td><td>21.4</td></tr><tr><td>59680.36458</td><td>16.1</td><td>16.3</td><td>-7.5</td><td>19.3</td></tr><tr><td>59700.30193</td><td>-218.9</td><td>15.7</td><td>-219.2</td><td>19.7</td></tr><tr><td>59700.31948</td><td>-243.2</td><td>16.3</td><td>-240.5</td><td>20.4</td></tr><tr><td>59700.33924</td><td>-206.4</td><td>18.2</td><td>-240.7</td><td>22.1</td></tr><tr><td>59722.16463</td><td>76.9</td><td>14.8</td><td>98.8</td><td>24</td></tr><tr><td>59722.18212</td><td>42.6</td><td>15.5</td><td>43.8</td><td>18.5</td></tr><tr><td>59724.34289</td><td>-187.2</td><td>19.1</td><td>-239.4</td><td>32.6</td></tr><tr><td>59724.36042</td><td>-188.9</td><td>20.7</td><td>-203.4</td><td>32.4</td></tr><tr><td>59724.38436</td><td>-226.2</td><td>33.3</td><td>-253.4</td><td>40.6</td></tr><tr><td>59740.19389</td><td>150.5</td><td>22.2</td><td>160.2</td><td>33.4</td></tr><tr><td>59740.21175</td><td>146.4</td><td>28</td><td>164.2</td><td>33</td></tr><tr><td>59740.30736</td><td>-20</td><td>16.1</td><td>3.4</td><td>22.1</td></tr><tr><td>59740.32483</td><td>-60.9</td><td>17.1</td><td>-86.1</td><td>23</td></tr><tr><td>59740.38604</td><td>-149.9</td><td>19.2</td><td>-161.6</td><td>21.7</td></tr><tr><td>59740.40352</td><td>-228.8</td><td>23.8</td><td>-261.7</td><td>27.9</td></tr></tbody></table>

corrections from this method are generally similar to, but sometimes more accurate than, those from the night sky lines.

We also refit the object radial velocities with RVSPECFIT (Koposov et al. 2011; Koposov 2019), using a library of PHOENIX synthetic templates (Allard 2016) of varying metallicity, temperature, surface gravity, [α/Fe] abundance, as well as allowing for rotation. As described in Section 1 companion surface heating complicates the measurement; the inferred velocity does not necessarily track the true centre-of-mass velocity, rather the centre of light associated with a specific line. This is clearly reflected by the differing K<sub>2</sub> amplitudes determined in Au et al. (2023), and updated here in Table 3, when considering the full spectrum versus only the Mg β triplet (a similar treatment is given in Linares, Shahbaz & Casares 2018).

Hence for each spectrum we performed two fits: the first over the entire range of the optical spectrum with measurable absorption lines (4000–6800 Å) and the second solely in the region of the Mg β line. Overall, the inferred velocities from this method are consistent with those obtained from cross-correlation with an appropriate template over a comparable wavelength range.

## 3 PHOTOMETRIC MODELLING

The optical light-curve modelling performed here utilized the binary stellar synthesis code ICARUS (Breton et al. 2012), with some novel

modifications. As such, the procedure followed is comparable, though not identical, to the modelling performed in similar analyses (Breton et al. 2013; Draghis et al. 2019; Stringer et al. 2021; Kennedy et al. 2022; Mata Sánchez et al. 2023). Here the specific procedure and priors used for this system will be described (see Breton et al. 2012, for a more in-depth description of ICARUS).

### 3.1 Surface heating models

Compared to previous uses of ICARUS, not limited to those cited above, here we have amended the gravity darkening prescription applied to the companion’s surface. Previously the temperature of companion surface element *i*, $T_i$, **before** irradiation was calculated as

$$T_i = T_{\text{base}} \left( \frac{g_i}{g_{\text{pole}}} \right)^{\beta}, \quad (1)$$

where $T_{\text{base}}$ is the ICARUS input parameter specifying the temperature at the pole of the star, $g_i$ is the surface gravity at surface element *i*, $g_{\text{pole}}$ is the surface gravity at the pole of the star, and $\beta$ is the gravity darkening coefficient. This equation still applies here, though its deployment differs in two significant ways:

(i) We assume the companion’s atmosphere heat transfer *close to the surface* is radiative, as opposed to convective. A radiative gravity darkening coefficient ($\beta$) of 0.25 was used, as opposed to the usual 0.08 used for a convective atmosphere (Breton et al. 2013).

(ii) We include the option to apply gravity darkening after irradiation and heat redistribution on the heated companion surface. This differs from the previously standard ICARUS behaviour to gravity darken the base (singular temperature) companion surface before heating effects are considered.

We found that these changes improve our model fits substantially and are physically motivated by a number of new insights we gained on the stellar physics. For the first assumption, following Zilles et al. (2020), we expect the inner photosphere of the companion to be convective where the Schwarzschild criterion is satisfied, and radiative towards the surface. Therefore the gravity darkening prescription for the photosphere surface should follow the radiative law. Espinosa Lara & Rieutord (2012) also demonstrated that tidally distorted low-mass, convective stars should in fact present gravity darkening coefficients in the interval [0.20,0.25], with spider-like companions being at the upper end of this range.

Though this latter work does not include the effects of irradiation, there is a strong possibility that the irradiation impinging on to J1910’s companion, and other spider companions, leads to deep heating of their photosphere. This is in contrast to our previous application of gravity darkening before irradiation, which implicitly assumed it was only superficial. The fact that spectral lines in these systems are generally absorption features (except for a few emission-line features which are likely connected to outflowing material)
<table><caption>**Table 2.** Time and phase coverage for ULTRACAM photometry obtained of J1910. The phase coverage, calculated with the timing ephemeris provided in Table 1, corresponds with the phase axis of Fig. 1. The *g*<sub>s</sub> filter is split from the other two due to the exclusion of irreducible data for the 2022 June 28 night.</caption><thead><tr><th>Start time (UTC)</th><th>Observation length (h)</th><th colspan="2">Phase coverage</th></tr><tr><th></th><th></th><th>r<sub>s</sub>, u<sub>s</sub></th><th>g<sub>s</sub></th></tr></thead><tbody><tr><td>2022 June 29 (04:11:44)</td><td>4.25</td><td>0.90–1.38</td><td>1.10–1.38</td></tr><tr><td>2022 July 1 (01:31:58)</td><td>5.0</td><td>0.30–0.89</td><td>0.30–0.89</td></tr></tbody></table>

<table><caption>**Table 3.** Updated radial velocities (RV) of PSR J1910-5320 from SOAR for both the full spectrum and targeting just the Mg β triplet.</caption><thead><tr><th rowspan="2">BJD (d)</th><th colspan="2">Full spectrum</th><th colspan="2">Mg β</th></tr><tr><th>RV</th><th>ΔRV (km s<sup>-1</sup>)</th><th>RV</th><th>ΔRV (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>59679.32228</td><td>-17.8</td><td>21.7</td><td>-28.9</td><td>25.3</td></tr><tr><td>59679.34012</td><td>-58.4</td><td>20.6</td><td>-72.4</td><td>26.5</td></tr><tr><td>59679.35806</td><td>-141.8</td><td>18</td><td>-164.2</td><td>26.3</td></tr><tr><td>59680.32745</td><td>80.7</td><td>17.5</td><td>82.1</td><td>21.1</td></tr><tr><td>59680.34496</td><td>67.5</td><td>19.1</td><td>79</td><td>21.4</td></tr><tr><td>59680.36458</td><td>16.1</td><td>16.3</td><td>-7.5</td><td>19.3</td></tr><tr><td>59700.30193</td><td>-218.9</td><td>15.7</td><td>-219.2</td><td>19.7</td></tr><tr><td>59700.31948</td><td>-243.2</td><td>16.3</td><td>-240.5</td><td>20.4</td></tr><tr><td>59700.33924</td><td>-206.4</td><td>18.2</td><td>-240.7</td><td>22.1</td></tr><tr><td>59722.16463</td><td>76.9</td><td>14.8</td><td>98.8</td><td>24</td></tr><tr><td>59722.18212</td><td>42.6</td><td>15.5</td><td>43.8</td><td>18.5</td></tr><tr><td>59724.34289</td><td>-187.2</td><td>19.1</td><td>-239.4</td><td>32.6</td></tr><tr><td>59724.36042</td><td>-188.9</td><td>20.7</td><td>-203.4</td><td>32.4</td></tr><tr><td>59724.38436</td><td>-226.2</td><td>33.3</td><td>-253.4</td><td>40.6</td></tr><tr><td>59740.19389</td><td>150.5</td><td>22.2</td><td>160.2</td><td>33.4</td></tr><tr><td>59740.21175</td><td>146.4</td><td>28</td><td>164.2</td><td>33</td></tr><tr><td>59740.30736</td><td>-20</td><td>16.1</td><td>3.4</td><td>22.1</td></tr><tr><td>59740.32483</td><td>-60.9</td><td>17.1</td><td>-86.1</td><td>23</td></tr><tr><td>59740.38604</td><td>-149.9</td><td>19.2</td><td>-161.6</td><td>21.7</td></tr><tr><td>59740.40352</td><td>-228.8</td><td>23.8</td><td>-261.7</td><td>27.9</td></tr></tbody></table>

corrections from this method are generally similar to, but sometimes more accurate than, those from the night sky lines.

We also refit the object radial velocities with RVSPECFIT (Koposov et al. 2011; Koposov 2019), using a library of PHOENIX synthetic templates (Allard 2016) of varying metallicity, temperature, surface gravity, [α/Fe] abundance, as well as allowing for rotation. As described in Section 1 companion surface heating complicates the measurement; the inferred velocity does not necessarily track the true centre-of-mass velocity, rather the centre of light associated with a specific line. This is clearly reflected by the differing K<sub>2</sub> amplitudes determined in Au et al. (2023), and updated here in Table 3, when considering the full spectrum versus only the Mg β triplet (a similar treatment is given in Linares, Shahbaz & Casares 2018).

Hence for each spectrum we performed two fits: the first over the entire range of the optical spectrum with measurable absorption lines (4000–6800 Å) and the second solely in the region of the Mg β line. Overall, the inferred velocities from this method are consistent with those obtained from cross-correlation with an appropriate template over a comparable wavelength range.

## 3 PHOTOMETRIC MODELLING

The optical light-curve modelling performed here utilized the binary stellar synthesis code ICARUS (Breton et al. 2012), with some novel

modifications. As such, the procedure followed is comparable, though not identical, to the modelling performed in similar analyses (Breton et al. 2013; Draghis et al. 2019; Stringer et al. 2021; Kennedy et al. 2022; Mata Sánchez et al. 2023). Here the specific procedure and priors used for this system will be described (see Breton et al. 2012, for a more in-depth description of ICARUS).

### 3.1 Surface heating models

Compared to previous uses of ICARUS, not limited to those cited above, here we have amended the gravity darkening prescription applied to the companion’s surface. Previously the temperature of companion surface element *i*, $T_i$, **before** irradiation was calculated as

$$T_i = T_{\text{base}} \left( \frac{g_i}{g_{\text{pole}}} \right)^{\beta}, \qquad (1)$$

where $T_{\text{base}}$ is the ICARUS input parameter specifying the temperature at the pole of the star, $g_i$ is the surface gravity at surface element *i*, $g_{\text{pole}}$ is the surface gravity at the pole of the star, and $\beta$ is the gravity darkening coefficient. This equation still applies here, though its deployment differs in two significant ways:

(i) We assume the companion’s atmosphere heat transfer *close to the surface* is radiative, as opposed to convective. A radiative gravity darkening coefficient ($\beta$) of 0.25 was used, as opposed to the usual 0.08 used for a convective atmosphere (Breton et al. 2013).

(ii) We include the option to apply gravity darkening after irradiation and heat redistribution on the heated companion surface. This differs from the previously standard ICARUS behaviour to gravity darken the base (singular temperature) companion surface before heating effects are considered.

We found that these changes improve our model fits substantially and are physically motivated by a number of new insights we gained on the stellar physics. For the first assumption, following Zilles et al. (2020), we expect the inner photosphere of the companion to be convective where the Schwarzschild criterion is satisfied, and radiative towards the surface. Therefore the gravity darkening prescription for the photosphere surface should follow the radiative law. Espinosa Lara & Rieutord (2012) also demonstrated that tidally distorted low-mass, convective stars should in fact present gravity darkening coefficients in the interval [0.20,0.25], with spider-like companions being at the upper end of this range.

Though this latter work does not include the effects of irradiation, there is a strong possibility that the irradiation impinging on to J1910’s companion, and other spider companions, leads to deep heating of their photosphere. This is in contrast to our previous application of gravity darkening before irradiation, which implicitly assumed it was only superficial. The fact that spectral lines in these systems are generally absorption features (except for a few emission-line features which are likely connected to outflowing material)
indicates that irradiation is deposited deep enough for no substantial thermal inversion to occur as is seen in the case of cataclysmic variables where the shallow heating is caused by UV photons from a hot white dwarf. It then follows that the irradiating flux should be considered a fundamental aspect of the surface temperature profile, and as such gravity darkened along with the rest. As the exact depth of the heating in J1910 is unclear and a full theoretical treatment of its effect on gravity darkening not available at the moment, we opted to test both pre- and post-irradiation gravity darkening models for completeness.

The parameters fit for using *ICARUS* depended on the surface heating model applied. The most basic model, direct heating (DH), applies symmetrical irradiation on to the companion's inner face, locked towards the pulsar. The parameters fit for this model constitute our fundamental set: the systemic velocity $\gamma$, the interstellar reddening $E(B-V)$, the system inclination $i$, the Roche-lobe filling factor $f_{RL}^*$,<fnref n="3" /> the base and irradiating temperatures $T_{base}$ and $T_{irr}$, the distance $d$, and the projected radial velocity amplitude of the companion $K_2$.

Heat redistribution across the stellar surface was also considered, as set out in Voisin et al. (2020). For an irradiated companion face with temperature differences between the dayside and nightside, diffusion of heat from the irradiated face can be expected. In our models this is accounted for by adding two parameters to our 'fundamental' parameter set: $\kappa$, which parametrizes the amplitude of the diffusion effect, and $\Gamma$, which governs the temperature dependence of the diffusion (Stringer et al. 2021). In this case, we have elected not to include $\Gamma$. Trial fits including it regularly found very little constraint on it, and those without obtained a better Bayesian evidence without significant effect on other parameters.

Heat redistribution models can also account for asymmetrical light curves, found for a number of spiders (Stappers et al. 2001; Romani & Sanchez 2016; Linares, Shahbaz & Casares 2018; Kandel & Romani 2020; Romani et al. 2021; Stringer et al. 2021), whereby light curves at not symmetric between the half orbits centred on the companion's ascending and descending nodes. Three main approaches have usually been implemented to account for this:

(i) A convective wind following a certain latitudinal profile, with strength parametrized by $C_{amp}$.

(ii) A surface hot/cold spot with fitted temperature, size, and position (e.g. Clark et al. 2021).

(iii) Re-distribution of irradiating flux by an extended, swept-back intra-binary shock (Romani & Sanchez 2016) and/or magnetic ducting (Sanchez & Romani 2017).

These models account for asymmetry by shifting or adding flux on to one side of the companion's inner face, such that more/less flux is seen at ingress/egress to the companion's superior conjunction. In this work we have focused on using diffusion and convection (D+C) models to redistribute heat across the companion's surface. While hot spots are well supported in literature and physically (Sanchez & Romani 2017), in the present case spot models invariably placed the spot, given the modelled inclination, largely out of sight on the companion's surface at all orbital phases. We took this as an indication that a spot model was not suitable for J1910.

The parameters set for each model were sampled and constrained by channelling *ICARUS* through *DYNESTY* (Speagle 2020), a Python implementation of a dynamic nested sampling Bayesian parameter and evidence estimation algorithm (Skilling 2004; Feroz & Hobson

2008; Feroz, Hobson & Bridges 2009; Feroz et al. 2019). Nested sampling algorithms provide the Bayesian evidence of a model, Z, a useful advantage over a classic Monte Carlo Markov Chain algorithms. Allowing for the calculation of the Bayes factor,

$$B_{1,2} = \frac{Z_1}{Z_2},$$

between two models enables one to determine which is favoured; $B_{1,2} > 1$ suggests model 1 is preferable, whereas $B_{1,2} < 1$ would prefer model 2 (Jeffreys 1939). The basic procedure on a given iteration of the nested sampler, using only the optical photometry, first selects a set of samples from the provided priors, passing them into *ICARUS*. The likelihood is calculated from the $\chi^2$ fit of the observed photometry and the simulated light curves generated given sampled parameters.

## 3.2 Priors

Careful consideration must be given to the choice of priors for our models and, where possible, they should be strongly motivated by physical or geometric constraints or, in the case of $K_2$, the use of complementary independent data (Au et al. 2023, Section 3.3). The main priors used here were as follows:

(i) A Gaussian prior applied to $E(B - V)$, centred on the reddening provided by the dust maps of Schlafy & Finkbeiner (2011): $0.0596 \pm 0.0033$.

(ii) A simple $\sin(i)$ prior applied to $i$, corresponding to an isotropic distribution of orbital angular momentum vectors.

(iii) A distance prior constructed using the same procedure as in Clark et al. (2021) and Kennedy et al. (2022). This combines the expected density of Galactic MSPs along the line of sight to J1910 (Levin et al. 2013), the transverse velocity distribution for binary MSPs in the ATNF Pulsar catalogue (Manchester et al. 2005), and the *Gaia* DR3 parallax (Gaia Collaboration 2023). Additional constraint can be provided by the DM inferred from radio timing using the Galactic electron density model (Yao, Manchester & Wang 2017, YMW16). In the present case, we have opted not to employ it. The DM distance is not equally reliable for all lines of sight, and the distance inferred from the DM ($0.92 \pm 0.49$ kpc) is much smaller, and less reliable, than that from the *Gaia* parallax ($6.8 \pm 3.9$ kpc). Yao, Manchester & Wang (2017) themselves compiled a list of pulsars with independent distance measurements both underestimated and overestimated by their model, therefore an underestimation from it for J1910 is not entirely unexpected.

## 3.3 Spectroscopic $K_2$ constraint

Given the very high-precision timing measurement of the pulsar's projected velocity amplitude, any measurement of the companion's $K_2$ determines the mass ratio $q$, and then provides a constraint on the masses via the mass function of the system. $K_2$ is typically measured from the Doppler motion of absorption lines over the orbit, to which a centre-of-light correction must be applied.

Previous iterations of *ICARUS* have allowed for the incorporation of spectroscopic data in various ways. Clark et al. (2021) calculated an average of companion surface element velocities (simulated as part of *ICARUS*) over the orbit, weighted by their flux to compensate for centre-of-light effects in an approximate manner. The resulting model radial velocities were subtracted from the observed radial velocities, and the overall model penalized according to the resulting likelihood. Kennedy et al. (2022) used a self-consistent procedure,

<footnote n="3">³Calculated as $r_N/r_{L1}$, where $r_N$ is the distance from the companion's barycentre to its nose, and $r_{L1}$ is the distance from the barycentre to the L1 point.</footnote>
indicates that irradiation is deposited deep enough for no substantial thermal inversion to occur as is seen in the case of cataclysmic variables where the shallow heating is caused by UV photons from a hot white dwarf. It then follows that the irradiating flux should be considered a fundamental aspect of the surface temperature profile, and as such gravity darkened along with the rest. As the exact depth of the heating in J1910 is unclear and a full theoretical treatment of its effect on gravity darkening not available at the moment, we opted to test both pre- and post-irradiation gravity darkening models for completeness.

The parameters fit for using *ICARUS* depended on the surface heating model applied. The most basic model, direct heating (DH), applies symmetrical irradiation on to the companion's inner face, locked towards the pulsar. The parameters fit for this model constitute our fundamental set: the systemic velocity $\gamma$, the interstellar reddening $E(B-V)$, the system inclination $i$, the Roche-lobe filling factor $f_{RL}^*$,<fnref n="3" /> the base and irradiating temperatures $T_{base}$ and $T_{irr}$, the distance $d$, and the projected radial velocity amplitude of the companion $K_2$.

Heat redistribution across the stellar surface was also considered, as set out in Voisin et al. (2020). For an irradiated companion face with temperature differences between the dayside and nightside, diffusion of heat from the irradiated face can be expected. In our models this is accounted for by adding two parameters to our 'fundamental' parameter set: $\kappa$, which parametrizes the amplitude of the diffusion effect, and $\Gamma$, which governs the temperature dependence of the diffusion (Stringer et al. 2021). In this case, we have elected not to include $\Gamma$. Trial fits including it regularly found very little constraint on it, and those without obtained a better Bayesian evidence without significant effect on other parameters.

Heat redistribution models can also account for asymmetrical light curves, found for a number of spiders (Stappers et al. 2001; Romani & Sanchez 2016; Linares, Shahbaz & Casares 2018; Kandel & Romani 2020; Romani et al. 2021; Stringer et al. 2021), whereby light curves at not symmetric between the half orbits centred on the companion's ascending and descending nodes. Three main approaches have usually been implemented to account for this:

(i) A convective wind following a certain latitudinal profile, with strength parametrized by $C_{amp}$.

(ii) A surface hot/cold spot with fitted temperature, size, and position (e.g. Clark et al. 2021).

(iii) Re-distribution of irradiating flux by an extended, swept-back intra-binary shock (Romani & Sanchez 2016) and/or magnetic ducting (Sanchez & Romani 2017).

These models account for asymmetry by shifting or adding flux on to one side of the companion's inner face, such that more/less flux is seen at ingress/egress to the companion's superior conjunction. In this work we have focused on using diffusion and convection (D+C) models to redistribute heat across the companion's surface. While hot spots are well supported in literature and physically (Sanchez & Romani 2017), in the present case spot models invariably placed the spot, given the modelled inclination, largely out of sight on the companion's surface at all orbital phases. We took this as an indication that a spot model was not suitable for J1910.

The parameters set for each model were sampled and constrained by channelling *ICARUS* through *DYNESTY* (Speagle 2020), a Python implementation of a dynamic nested sampling Bayesian parameter and evidence estimation algorithm (Skilling 2004; Feroz & Hobson

2008; Feroz, Hobson & Bridges 2009; Feroz et al. 2019). Nested sampling algorithms provide the Bayesian evidence of a model, Z, a useful advantage over a classic Monte Carlo Markov Chain algorithms. Allowing for the calculation of the Bayes factor,

$$B_{1,2} = \frac{Z_1}{Z_2},$$

between two models enables one to determine which is favoured; $B_{1,2} > 1$ suggests model 1 is preferable, whereas $B_{1,2} < 1$ would prefer model 2 (Jeffreys 1939). The basic procedure on a given iteration of the nested sampler, using only the optical photometry, first selects a set of samples from the provided priors, passing them into *ICARUS*. The likelihood is calculated from the $\chi^2$ fit of the observed photometry and the simulated light curves generated given sampled parameters.

## 3.2 Priors

Careful consideration must be given to the choice of priors for our models and, where possible, they should be strongly motivated by physical or geometric constraints or, in the case of $K_2$, the use of complementary independent data (Au et al. 2023, Section 3.3). The main priors used here were as follows:

(i) A Gaussian prior applied to $E(B-V)$, centred on the reddening provided by the dust maps of Schlafy & Finkbeiner (2011): $0.0596 \pm 0.0033$.

(ii) A simple $\sin(i)$ prior applied to $i$, corresponding to an isotropic distribution of orbital angular momentum vectors.

(iii) A distance prior constructed using the same procedure as in Clark et al. (2021) and Kennedy et al. (2022). This combines the expected density of Galactic MSPs along the line of sight to J1910 (Levin et al. 2013), the transverse velocity distribution for binary MSPs in the ATNF Pulsar catalogue (Manchester et al. 2005), and the *Gaia* DR3 parallax (Gaia Collaboration 2023). Additional constraint can be provided by the DM inferred from radio timing using the Galactic electron density model (Yao, Manchester & Wang 2017, YMW16). In the present case, we have opted not to employ it. The DM distance is not equally reliable for all lines of sight, and the distance inferred from the DM ($0.92 \pm 0.49$ kpc) is much smaller, and less reliable, than that from the *Gaia* parallax ($6.8 \pm 3.9$ kpc). Yao, Manchester & Wang (2017) themselves compiled a list of pulsars with independent distance measurements both underestimated and overestimated by their model, therefore an underestimation from it for J1910 is not entirely unexpected.

## 3.3 Spectroscopic $K_2$ constraint

Given the very high-precision timing measurement of the pulsar's projected velocity amplitude, any measurement of the companion's $K_2$ determines the mass ratio $q$, and then provides a constraint on the masses via the mass function of the system. $K_2$ is typically measured from the Doppler motion of absorption lines over the orbit, to which a centre-of-light correction must be applied.

Previous iterations of *ICARUS* have allowed for the incorporation of spectroscopic data in various ways. Clark et al. (2021) calculated an average of companion surface element velocities (simulated as part of *ICARUS*) over the orbit, weighted by their flux to compensate for centre-of-light effects in an approximate manner. The resulting model radial velocities were subtracted from the observed radial velocities, and the overall model penalized according to the resulting likelihood. Kennedy et al. (2022) used a self-consistent procedure,

<footnote n="3">³Calculated as $r_N/r_{L1}$, where $r_N$ is the distance from the companion's barycentre to its nose, and $r_{L1}$ is the distance from the barycentre to the L1 point.</footnote>
<table><caption>**Table 4.** Posterior parameter results from photometry and radial velocity curve fitting. Results are split into the two key models used: DH, which employs no heat redistribution, and diffusion + convection (D+C). These are subsequently split by the gravity darkening prescription, pre- or post-irradiation (including heat redistribution effects). Note the *ICARUS* parameters $T_B$ and $T_I$ do not reflect the 'true' physical conditions on the companion's surface. Rather, $T_S$ and $T_I$ average the visible surface element temperatures at the companion's superior and inferior conjunctions, respectively. $T_{day}$ and $T_{night}$ then provide the intrinsic temperatures of the day and night sides, again averaging surface element temperatures assuming an edge on inclination. Quoted uncertainties correspond to 68 per cent confidence intervals. The (log) Bayesian evidence (ln Z) produced by *DYNESTY* is used to calculate the Bayes factors between a given model and a reference one, chosen to be DH with pre-irradiation gravity darkening.</caption><thead><tr><th rowspan="2">Param</th><th colspan="2">DH</th><th colspan="2">D+C</th></tr><tr><th>Pre</th><th>Post</th><th>Pre</th><th>Post</th></tr></thead><tbody><tr><td colspan="5" style="text-align:center; font-weight:bold;">Icarus parameters</td></tr><tr><td>Interstellar reddening, E(B-V)</td><td>0.060<sup>+0.003</sup><sub>-0.003</sub></td><td>0.060<sup>+0.003</sup><sub>-0.003</sub></td><td>0.060<sup>+0.003</sup><sub>-0.003</sub></td><td>0.060<sup>+0.003</sup><sub>-0.003</sub></td></tr><tr><td>Inclination angle, i (deg)</td><td>54<sup>+3</sup><sub>-3</sub></td><td>52<sup>+2</sup><sub>-2</sub></td><td>46<sup>+1</sup><sub>-1</sub></td><td>45<sup>+1</sup><sub>-1</sub></td></tr><tr><td>Roche-lobe filling factor, f<sub>RL</sub></td><td>0.818<sup>+0.007</sup><sub>-0.008</sub></td><td>0.762<sup>+0.005</sup><sub>-0.005</sub></td><td>0.838<sup>+0.008</sup><sub>-0.010</sub></td><td>0.782<sup>+0.005</sup><sub>-0.005</sub></td></tr><tr><td>Base temperature T<sub>base</sub> (K)</td><td>5460<sup>+40</sup><sub>-40</sub></td><td>5360<sup>+42</sup><sub>-50</sub></td><td>5310<sup>+40</sup><sub>-40</sub></td><td>5200<sup>+40</sup><sub>-40</sub></td></tr><tr><td>Irradiating temperature T<sub>irr</sub> (K)</td><td>6060<sup>+60</sup><sub>-60</sub></td><td>6700<sup>+70</sup><sub>-70</sub></td><td>6240<sup>+80</sup><sub>-70</sub></td><td>6760<sup>+70</sup><sub>-70</sub></td></tr><tr><td>Distance, d (kpc)</td><td>4.2<sup>+0.2</sup><sub>-0.2</sub></td><td>4.0<sup>+0.1</sup><sub>-0.1</sub></td><td>4.8<sup>+0.2</sup><sub>-0.2</sub></td><td>4.43<sup>+0.2</sup><sub>-0.1</sub></td></tr><tr><td>Companion radial velocity amplitude, K<sub>2</sub> (km s<sup>-1</sup>)</td><td>206<sup>+7</sup><sub>-7</sub></td><td>197<sup>+8</sup><sub>-7</sub></td><td>216<sup>+8</sup><sub>-7</sub></td><td>200<sup>+8</sup><sub>-8</sub></td></tr><tr><td colspan="5" style="text-align:center; font-weight:bold;">Heating parameters</td></tr><tr><td>Diffusion coefficient, κ</td><td>—</td><td>—</td><td>900<sup>+1200</sup><sub>-600</sub></td><td>80<sup>+100</sup><sub>-50</sub></td></tr><tr><td>Convection amplitude C<sub>amp</sub></td><td>—</td><td>—</td><td>-1700<sup>+100</sup><sub>-100</sub></td><td>-1760<sup>+60</sup><sub>-60</sub></td></tr><tr><td colspan="5" style="text-align:center; font-weight:bold;">Derived parameters</td></tr><tr><td>Mass ratio, q</td><td>3.4<sup>+0.1</sup><sub>-0.1</sub></td><td>3.2<sup>+0.1</sup><sub>-0.1</sub></td><td>3.6<sup>+0.1</sup><sub>-0.1</sub></td><td>3.3<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>Pulsar mass, M<sub>P</sub> (M<sub>☉</sub>)</td><td>1.0<sup>+0.1</sup><sub>-0.1</sub></td><td>1.0<sup>+0.1</sup><sub>-0.1</sub></td><td>1.6<sup>+0.2</sup><sub>-0.2</sub></td><td>1.4<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>Companion mass, M<sub>C</sub> (M<sub>☉</sub>)</td><td>0.29<sup>+0.03</sup><sub>-0.04</sub></td><td>0.29<sup>+0.03</sup><sub>-0.03</sub></td><td>0.45<sup>+0.04</sup><sub>-0.04</sub></td><td>0.43<sup>+0.04</sup><sub>-0.03</sub></td></tr><tr><td>Observed temp. at companion superior conjunction, T<sub>S</sub> (K)</td><td>6170<sup>+50</sup><sub>-50</sub></td><td>6110<sup>+50</sup><sub>-60</sub></td><td>6060<sup>+40</sup><sub>-40</sub></td><td>5940<sup>+50</sup><sub>-60</sub></td></tr><tr><td>Observed temp. at companion inferior conjunction, T<sub>I</sub> (K)</td><td>5140<sup>+30</sup><sub>-30</sub></td><td>5100<sup>+30</sup><sub>-40</sub></td><td>5070<sup>+30</sup><sub>-30</sub></td><td>4990<sup>+40</sup><sub>-40</sub></td></tr><tr><td>Day-side temperature, T<sub>day</sub> (K)</td><td>6350<sup>+40</sup><sub>-40</sub></td><td>6300<sup>+50</sup><sub>-50</sub></td><td>6340<sup>+40</sup><sub>-40</sub></td><td>6220<sup>+50</sup><sub>-60</sub></td></tr><tr><td>Night-side temperature, T<sub>night</sub> (K)</td><td>5050<sup>+40</sup><sub>-40</sub></td><td>5000<sup>+40</sup><sub>-40</sub></td><td>4930<sup>+30</sup><sub>-30</sub></td><td>4840<sup>+40</sup><sub>-40</sub></td></tr><tr><td>Volume-averaged filling factor, f<sub>V</sub></td><td>0.947<sup>+0.004</sup><sub>-0.004</sub></td><td>0.911<sup>+0.004</sup><sub>-0.004</sub></td><td>0.957<sup>+0.004</sup><sub>-0.006</sub></td><td>0.924<sup>+0.004</sup><sub>-0.003</sub></td></tr><tr><td>Irradiation efficiency, ε</td><td>0.22<sup>+0.03</sup><sub>-0.03</sub></td><td>0.32<sup>+0.03</sup><sub>-0.03</sub></td><td>0.34<sup>+0.03</sup><sub>-0.03</sub></td><td>0.42<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>Transverse velocity, v<sub>T</sub> (km s<sup>-1</sup>)</td><td>140<sup>+6</sup><sub>-7</sub></td><td>132<sup>+5</sup><sub>-5</sub></td><td>158<sup>+6</sup><sub>-5</sub></td><td>144<sup>+5</sup><sub>-5</sub></td></tr><tr><td colspan="5" style="text-align:center; font-weight:bold;">Fit statistics</td></tr><tr><td>Photometry χ<sup>2</sup> (3446 datapoints)</td><td>4889.585</td><td>5106.263</td><td>3611.385</td><td>3637.874</td></tr><tr><td>Radial velocity χ<sup>2</sup> (20 datapoints)</td><td>29.144</td><td>29.830</td><td>30.456</td><td>29.132</td></tr><tr><td>log-Evidence, Z</td><td>-2366.3 ± 0.2</td><td>-2442.1 ± 0.2</td><td>-1851.6 ± 0.2</td><td>-1865.8 ± 0.2</td></tr><tr><td>log Bayes Factor versus Convective</td><td>119</td><td>43</td><td>634</td><td>619</td></tr><tr><td>log Bayes Factor versus Radiative</td><td>0</td><td>-76</td><td>514</td><td>501</td></tr></tbody></table>

where observed spectra were directly fitted to simulated spectra generated by *ICARUS* from ATLAS9 (Castelli & Kurucz 2003) atmosphere grids to produce a likelihood. This method intrinsically overcomes the centre-of-light issue, as irradiation is implicit in the generated model spectra. There is, however, a significant computational cost associated with simulating full model spectra and a potential risk for the fitting to try and reproduce features of the spectrum which are not well accounted for by the atmosphere model.

In this work a middle ground between the two methods described above was used, balancing adequate simulation of the spectra with computational expense. As with the self-consistent spectroscopy modelling of Kennedy et al. (2022), here *ICARUS* is used to simulate spectra for each sample. However, these spectra were not directly compared with their observed counterparts, rather the radial velocities of the models were determined and compared to their experimental analogues. Specifically narrow, and thus inexpensive, spectra centred around the 5183 Å Mg β triplet were generated for each orbital phase covered by the SOAR/Goodman data set. The radial velocity for each phase was determined by cross-correlating the spectrum at a reference orbital phase (chosen to be that showing

the strongest line feature), thus providing a relative projected radial velocity curve. The likelihood between the observed and modelled radial velocities was then incorporated into the fitting procedure.

## 4 MODELLING RESULTS

Table 4 contains the results for the models considered and discussed above. These are split by heating model (DH or D+C) and subsequently by the prescription used to apply gravity darkening (pre- versus post-irradiation and heat redistribution effects). Select model fits for the photometry and radial velocity curve are shown in Fig. 2 and Fig. 3 respectively. Corner plots for the D+C models are available in Fig. B1 and Fig. B2. In both heating models a consistent trend emerges: post-irradiation gravity darkening finds a smaller projected companion velocity $K_2$. Before dissecting the differences between the pre- and post-irradiation gravity darkening, we can first get an overall picture of the parameters determined for this newly modelled system.

The DH models are presented for completeness; they do not constitute favourable models. The left-hand panels of Fig. 2 show the
<table><caption>**Table 4.** Posterior parameter results from photometry and radial velocity curve fitting. Results are split into the two key models used: DH, which employs no heat redistribution, and diffusion + convection (D+C). These are subsequently split by the gravity darkening prescription, pre- or post-irradiation (including heat redistribution effects). Note the *ICARUS* parameters $T_B$ and $T_I$ do not reflect the 'true' physical conditions on the companion's surface. Rather, $T_S$ and $T_I$ average the visible surface element temperatures at the companion's superior and inferior conjunctions, respectively. $T_{day}$ and $T_{night}$ then provide the intrinsic temperatures of the day and night sides, again averaging surface element temperatures assuming an edge on inclination. Quoted uncertainties correspond to 68 per cent confidence intervals. The (log) Bayesian evidence (ln Z) produced by *DYNESTY* is used to calculate the Bayes factors between a given model and a reference one, chosen to be DH with pre-irradiation gravity darkening.</caption><thead><tr><th rowspan="2">Param</th><th colspan="2">DH</th><th colspan="2">D+C</th></tr><tr><th>Pre</th><th>Post</th><th>Pre</th><th>Post</th></tr></thead><tbody><tr><td colspan="5" style="text-align:center; font-weight:bold;">Icarus parameters</td></tr><tr><td>Interstellar reddening, E(B-V)</td><td>0.060<sup>+0.003</sup><sub>-0.003</sub></td><td>0.060<sup>+0.003</sup><sub>-0.003</sub></td><td>0.060<sup>+0.003</sup><sub>-0.003</sub></td><td>0.060<sup>+0.003</sup><sub>-0.003</sub></td></tr><tr><td>Inclination angle, i (deg)</td><td>54<sup>+3</sup><sub>-3</sub></td><td>52<sup>+2</sup><sub>-2</sub></td><td>46<sup>+1</sup><sub>-1</sub></td><td>45<sup>+1</sup><sub>-1</sub></td></tr><tr><td>Roche-lobe filling factor, f<sub>RL</sub></td><td>0.818<sup>+0.007</sup><sub>-0.008</sub></td><td>0.762<sup>+0.005</sup><sub>-0.005</sub></td><td>0.838<sup>+0.008</sup><sub>-0.010</sub></td><td>0.782<sup>+0.005</sup><sub>-0.005</sub></td></tr><tr><td>Base temperature T<sub>base</sub> (K)</td><td>5460<sup>+40</sup><sub>-40</sub></td><td>5360<sup>+42</sup><sub>-50</sub></td><td>5310<sup>+40</sup><sub>-40</sub></td><td>5200<sup>+40</sup><sub>-40</sub></td></tr><tr><td>Irradiating temperature T<sub>irr</sub> (K)</td><td>6060<sup>+60</sup><sub>-60</sub></td><td>6700<sup>+70</sup><sub>-70</sub></td><td>6240<sup>+80</sup><sub>-70</sub></td><td>6760<sup>+70</sup><sub>-70</sub></td></tr><tr><td>Distance, d (kpc)</td><td>4.2<sup>+0.2</sup><sub>-0.2</sub></td><td>4.0<sup>+0.1</sup><sub>-0.1</sub></td><td>4.8<sup>+0.2</sup><sub>-0.2</sub></td><td>4.43<sup>+0.2</sup><sub>-0.1</sub></td></tr><tr><td>Companion radial velocity amplitude, K<sub>2</sub> (km s<sup>-1</sup>)</td><td>206<sup>+7</sup><sub>-7</sub></td><td>197<sup>+8</sup><sub>-7</sub></td><td>216<sup>+8</sup><sub>-7</sub></td><td>200<sup>+8</sup><sub>-8</sub></td></tr><tr><td colspan="5" style="text-align:center; font-weight:bold;">Heating parameters</td></tr><tr><td>Diffusion coefficient, κ</td><td>—</td><td>—</td><td>900<sup>+1200</sup><sub>-600</sub></td><td>80<sup>+100</sup><sub>-50</sub></td></tr><tr><td>Convection amplitude C<sub>amp</sub></td><td>—</td><td>—</td><td>-1700<sup>+100</sup><sub>-100</sub></td><td>-1760<sup>+60</sup><sub>-60</sub></td></tr><tr><td colspan="5" style="text-align:center; font-weight:bold;">Derived parameters</td></tr><tr><td>Mass ratio, q</td><td>3.4<sup>+0.1</sup><sub>-0.1</sub></td><td>3.2<sup>+0.1</sup><sub>-0.1</sub></td><td>3.6<sup>+0.1</sup><sub>-0.1</sub></td><td>3.3<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>Pulsar mass, M<sub>P</sub> (M<sub>☉</sub>)</td><td>1.0<sup>+0.1</sup><sub>-0.1</sub></td><td>1.0<sup>+0.1</sup><sub>-0.1</sub></td><td>1.6<sup>+0.2</sup><sub>-0.2</sub></td><td>1.4<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>Companion mass, M<sub>C</sub> (M<sub>☉</sub>)</td><td>0.29<sup>+0.03</sup><sub>-0.04</sub></td><td>0.29<sup>+0.03</sup><sub>-0.03</sub></td><td>0.45<sup>+0.04</sup><sub>-0.04</sub></td><td>0.43<sup>+0.04</sup><sub>-0.03</sub></td></tr><tr><td>Observed temp. at companion superior conjunction, T<sub>S</sub> (K)</td><td>6170<sup>+50</sup><sub>-50</sub></td><td>6110<sup>+50</sup><sub>-60</sub></td><td>6060<sup>+40</sup><sub>-40</sub></td><td>5940<sup>+50</sup><sub>-60</sub></td></tr><tr><td>Observed temp. at companion inferior conjunction, T<sub>I</sub> (K)</td><td>5140<sup>+30</sup><sub>-30</sub></td><td>5100<sup>+30</sup><sub>-40</sub></td><td>5070<sup>+30</sup><sub>-30</sub></td><td>4990<sup>+40</sup><sub>-40</sub></td></tr><tr><td>Day-side temperature, T<sub>day</sub> (K)</td><td>6350<sup>+40</sup><sub>-40</sub></td><td>6300<sup>+50</sup><sub>-50</sub></td><td>6340<sup>+40</sup><sub>-40</sub></td><td>6220<sup>+50</sup><sub>-60</sub></td></tr><tr><td>Night-side temperature, T<sub>night</sub> (K)</td><td>5050<sup>+40</sup><sub>-40</sub></td><td>5000<sup>+40</sup><sub>-40</sub></td><td>4930<sup>+30</sup><sub>-30</sub></td><td>4840<sup>+40</sup><sub>-40</sub></td></tr><tr><td>Volume-averaged filling factor, f<sub>V</sub></td><td>0.947<sup>+0.004</sup><sub>-0.004</sub></td><td>0.911<sup>+0.004</sup><sub>-0.004</sub></td><td>0.957<sup>+0.004</sup><sub>-0.006</sub></td><td>0.924<sup>+0.004</sup><sub>-0.003</sub></td></tr><tr><td>Irradiation efficiency, ε</td><td>0.22<sup>+0.03</sup><sub>-0.03</sub></td><td>0.32<sup>+0.03</sup><sub>-0.03</sub></td><td>0.34<sup>+0.03</sup><sub>-0.03</sub></td><td>0.42<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>Transverse velocity, v<sub>T</sub> (km s<sup>-1</sup>)</td><td>140<sup>+6</sup><sub>-7</sub></td><td>132<sup>+5</sup><sub>-5</sub></td><td>158<sup>+6</sup><sub>-5</sub></td><td>144<sup>+5</sup><sub>-5</sub></td></tr><tr><td colspan="5" style="text-align:center; font-weight:bold;">Fit statistics</td></tr><tr><td>Photometry χ<sup>2</sup> (3446 datapoints)</td><td>4889.585</td><td>5106.263</td><td>3611.385</td><td>3637.874</td></tr><tr><td>Radial velocity χ<sup>2</sup> (20 datapoints)</td><td>29.144</td><td>29.830</td><td>30.456</td><td>29.132</td></tr><tr><td>log-Evidence, Z</td><td>-2366.3 ± 0.2</td><td>-2442.1 ± 0.2</td><td>-1851.6 ± 0.2</td><td>-1865.8 ± 0.2</td></tr><tr><td>log Bayes Factor versus Convective</td><td>119</td><td>43</td><td>634</td><td>619</td></tr><tr><td>log Bayes Factor versus Radiative</td><td>0</td><td>-76</td><td>514</td><td>501</td></tr></tbody></table>

where observed spectra were directly fitted to simulated spectra generated by *ICARUS* from ATLAS9 (Castelli & Kurucz 2003) atmosphere grids to produce a likelihood. This method intrinsically overcomes the centre-of-light issue, as irradiation is implicit in the generated model spectra. There is, however, a significant computational cost associated with simulating full model spectra and a potential risk for the fitting to try and reproduce features of the spectrum which are not well accounted for by the atmosphere model.

In this work a middle ground between the two methods described above was used, balancing adequate simulation of the spectra with computational expense. As with the self-consistent spectroscopy modelling of Kennedy et al. (2022), here *ICARUS* is used to simulate spectra for each sample. However, these spectra were not directly compared with their observed counterparts, rather the radial velocities of the models were determined and compared to their experimental analogues. Specifically narrow, and thus inexpensive, spectra centred around the 5183 Å Mg β triplet were generated for each orbital phase covered by the SOAR/Goodman data set. The radial velocity for each phase was determined by cross-correlating the spectrum at a reference orbital phase (chosen to be that showing

the strongest line feature), thus providing a relative projected radial velocity curve. The likelihood between the observed and modelled radial velocities was then incorporated into the fitting procedure.

## 4 MODELLING RESULTS

Table 4 contains the results for the models considered and discussed above. These are split by heating model (DH or D+C) and subsequently by the prescription used to apply gravity darkening (pre- versus post-irradiation and heat redistribution effects). Select model fits for the photometry and radial velocity curve are shown in Fig. 2 and Fig. 3 respectively. Corner plots for the D+C models are available in Fig. B1 and Fig. B2. In both heating models a consistent trend emerges: post-irradiation gravity darkening finds a smaller projected companion velocity $K_2$. Before dissecting the differences between the pre- and post-irradiation gravity darkening, we can first get an overall picture of the parameters determined for this newly modelled system.

The DH models are presented for completeness; they do not constitute favourable models. The left-hand panels of Fig. 2 show the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Photometry fits produced by post-irradiation gravity darkening models presented in Table 4. The maximum a posteriori likelihood models have been selected. The left-hand panel shows the DH model, while the right-hand panel is diffusion + convection (D+C). The light curve data (Fig. 1) are shown in the corresponding colours, with model fits overlaid in black. Residuals for each band are shown below. Clearly visible between the two panels is the improvement in the residuals with the introduction of diffusion + convection to address the asymmetry in the light curve.</figcaption></figure>

post-irradiation gravity darkened DH model fit to the data. Paying attention to the residuals, the asymmetry in the light curve becomes clear. The model both overestimates the flux at the ingress to the optical maximum and underestimates the flux at the egress. The 12 reference stars used in ensemble photometry show no consistent excess corresponding to these orbital phases, thus it is safe to assume this is intrinsic light curve asymmetry. As such, the extremely low pulsar masses determined for both DH models can be safely discarded.

Our D+C models are much better than DH models at capturing the behaviour of the data and can account well for the asymmetry. The inferred $C_{\text{amp}}$ implies a convective surface wind blowing in the direction of the companion's rotation, and thus depositing heat towards the companion's leading edge. The improvement in the fit is reflected in the statistics provided in Table 4. The underlying reasons for changes in parameter values are far from trivial to pin down, but notable is a shift in *i* between the DH and D+C models, which implies a different inferred pulsar mass. Given a DH model will struggle to
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Photometry fits produced by post-irradiation gravity darkening models presented in Table 4. The maximum a posteriori likelihood models have been selected. The left-hand panel shows the DH model, while the right-hand panel is diffusion + convection (D+C). The light curve data (Fig. 1) are shown in the corresponding colours, with model fits overlaid in black. Residuals for each band are shown below. Clearly visible between the two panels is the improvement in the residuals with the introduction of diffusion + convection to address the asymmetry in the light curve.</figcaption></figure>

post-irradiation gravity darkened DH model fit to the data. Paying attention to the residuals, the asymmetry in the light curve becomes clear. The model both overestimates the flux at the ingress to the optical maximum and underestimates the flux at the egress. The 12 reference stars used in ensemble photometry show no consistent excess corresponding to these orbital phases, thus it is safe to assume this is intrinsic light curve asymmetry. As such, the extremely low pulsar masses determined for both DH models can be safely discarded.

Our D+C models are much better than DH models at capturing the behaviour of the data and can account well for the asymmetry. The inferred $C_{\text{amp}}$ implies a convective surface wind blowing in the direction of the companion's rotation, and thus depositing heat towards the companion's leading edge. The improvement in the fit is reflected in the statistics provided in Table 4. The underlying reasons for changes in parameter values are far from trivial to pin down, but notable is a shift in *i* between the DH and D+C models, which implies a different inferred pulsar mass. Given a DH model will struggle to
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Mg *β* radial velocity curve fit for post-irradiation gravity darkening diffusion + convection model. The top panel shows our model radial velocity points, blue, against the observed curve, red. The corresponding dashed lines are sinusoidal fits through each set, giving the parameters in the top corners. The grey solid line is the centre-of-mass radial velocity curve, using the underlying *K*₂ for the best-fitting model. Point wise residuals between the model and observed points are shown in the bottom panel.</figcaption></figure>

fit the amplitude of a asymmetric light curve it is unsurprising that *i*, which directly modulates the amplitude of an optical light curve, will be affected once heat redistribution is incorporated.

When compared with similar *ICARUS* modelling results involving asymmetric heat redistribution, J1910 is the only redback in which the heat is transferred to the leading edge (i.e. excess flux near descending node of the companion). PSRs J2215+5135 (Voisin et al. 2020), J1227–4853, and J1023+0038 (Stringer et al. 2021) all show excess flux towards the trailing edge of the light curve (i.e. excess flux near ascending node of the companion). Though we draw no major assertions from it, J1910 marks a notable departure from previously modelled redbacks.

## 4.1 Overall constraints

Considering now only the D+C models, a number of parameters agree across both gravity darkening options. The inclination remains consistent around 45°, with both models agreeing within their respective 68 per cent confidence interval. The irradiating temperatures in both models are consistently above 6000 K. More importantly, both models find average temperatures – where the temperatures across the visible surface are averaged in their fourth power, i.e. according to their bolometric luminosity, and weighted by the projected surface area – at the observed superior and inferior conjunctions that agree within their 68 per cent confidence intervals. This means that both models essentially reproduce the same colours in these parts of the light curves. From the lowest and highest points of the 68 per cent confidence regions, we find $4950 < T_{\text{l}} < 5100$. This is slightly lower than our expectation from the broad-band spectral energy distribution (SED) but within the allowed uncertainty (Au et al. 2023). $C_{\text{amp}}$ also agrees well for both which is expected given this parameters controls the asymmetry in the light curve.

Several parameters are not consistent between models, though we can still produce ‘ballpark’ educated guesses at their values. The

filling factors do vary between the models, but not over a large range, with both implying a significantly underfilling companion. Moving from the *ICARUS* parameter $f_{\text{RL}}$ to the volume-averaged filling factor we find an even smaller interval. Though significantly higher than the *ICARUS* parameter $f_{\text{RL}}$, these should still be interpreted as underfilling, particularly the post-irradiation gravity darkening case.

A key aim of light-curve modelling in spider systems is to constrain the pulsar mass. Fig. 4 shows a collection of spider mass measurements, with the masses determined for our D+C models shown in purple; the square and triangle denote the pre- and post-irradiation gravity darkening models, respectively. In this case we get a two moderate masses depending on the model chosen – none threaten the upper end observed pulsar masses and thus are useful to constrain the dense matter EoS on their own.

Linares (2019) collated a number of ‘supermassive’ neutron star mass measurements. The quality of our measurement is at a similar level to other spiders in this sample – especially those without independent constraints on either the inclination or companion mass. For example, PSR B1957+20’s recently updated mass constraint uses *γ*-ray eclipsing to provide hard constraints on the inclination (Clark et al. 2023a). We do not reach the same mass precision as Kennedy et al. (2022) or Romani et al. (2021), where the full, high-S/N spectroscopy has been used in constraining the model. The high-precision masses determined for relativistic NS–NS binaries, utilizing post-Keplerian parameters measured through pulsar timing, outperform the measurement here as do measurements for NS–WD binaries (see Lattimer 2012). The systematics inherent to spider light-curve modelling, namely the reliance on inferring a heating model for the surface, somewhat limit the precision we can expect to achieve. As these systematics are chiefly driven by irradiation, they are typically assumed to be lessened in redbacks when compared with black widows (Strader et al. 2019). However as J1910 is an irradiation-driven redback, significant surface heating must be accounted for. The precision of J1910’s mass measurement, as
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Mg *β* radial velocity curve fit for post-irradiation gravity darkening diffusion + convection model. The top panel shows our model radial velocity points, blue, against the observed curve, red. The corresponding dashed lines are sinusoidal fits through each set, giving the parameters in the top corners. The grey solid line is the centre-of-mass radial velocity curve, using the underlying *K*₂ for the best-fitting model. Point wise residuals between the model and observed points are shown in the bottom panel.</figcaption></figure>

fit the amplitude of a asymmetric light curve it is unsurprising that *i*, which directly modulates the amplitude of an optical light curve, will be affected once heat redistribution is incorporated.

When compared with similar *ICARUS* modelling results involving asymmetric heat redistribution, J1910 is the only redback in which the heat is transferred to the leading edge (i.e. excess flux near descending node of the companion). PSRs J2215+5135 (Voisin et al. 2020), J1227–4853, and J1023+0038 (Stringer et al. 2021) all show excess flux towards the trailing edge of the light curve (i.e. excess flux near ascending node of the companion). Though we draw no major assertions from it, J1910 marks a notable departure from previously modelled redbacks.

## 4.1 Overall constraints

Considering now only the D+C models, a number of parameters agree across both gravity darkening options. The inclination remains consistent around 45°, with both models agreeing within their respective 68 per cent confidence interval. The irradiating temperatures in both models are consistently above 6000 K. More importantly, both models find average temperatures – where the temperatures across the visible surface are averaged in their fourth power, i.e. according to their bolometric luminosity, and weighted by the projected surface area – at the observed superior and inferior conjunctions that agree within their 68 per cent confidence intervals. This means that both models essentially reproduce the same colours in these parts of the light curves. From the lowest and highest points of the 68 per cent confidence regions, we find $4950 < T_{\text{eff}} < 5100$. This is slightly lower than our expectation from the broad-band spectral energy distribution (SED) but within the allowed uncertainty (Au et al. 2023). $C_{\text{amp}}$ also agrees well for both which is expected given this parameters controls the asymmetry in the light curve.

Several parameters are not consistent between models, though we can still produce ‘ballpark’ educated guesses at their values. The

filling factors do vary between the models, but not over a large range, with both implying a significantly underfilling companion. Moving from the *ICARUS* parameter $f_{\text{RL}}$ to the volume-averaged filling factor we find an even smaller interval. Though significantly higher than the *ICARUS* parameter $f_{\text{RL}}$, these should still be interpreted as underfilling, particularly the post-irradiation gravity darkening case.

A key aim of light-curve modelling in spider systems is to constrain the pulsar mass. Fig. 4 shows a collection of spider mass measurements, with the masses determined for our D+C models shown in purple; the square and triangle denote the pre- and post-irradiation gravity darkening models, respectively. In this case we get a two moderate masses depending on the model chosen – none threaten the upper end observed pulsar masses and thus are useful to constrain the dense matter EoS on their own.

Linares (2019) collated a number of ‘supermassive’ neutron star mass measurements. The quality of our measurement is at a similar level to other spiders in this sample – especially those without independent constraints on either the inclination or companion mass. For example, PSR B1957+20’s recently updated mass constraint uses *γ*-ray eclipsing to provide hard constraints on the inclination (Clark et al. 2023a). We do not reach the same mass precision as Kennedy et al. (2022) or Romani et al. (2021), where the full, high-S/N spectroscopy has been used in constraining the model. The high-precision masses determined for relativistic NS–NS binaries, utilizing post-Keplerian parameters measured through pulsar timing, outperform the measurement here as do measurements for NS–WD binaries (see Lattimer 2012). The systematics inherent to spider light-curve modelling, namely the reliance on inferring a heating model for the surface, somewhat limit the precision we can expect to achieve. As these systematics are chiefly driven by irradiation, they are typically assumed to be lessened in redbacks when compared with black widows (Strader et al. 2019). However as J1910 is an irradiation-driven redback, significant surface heating must be accounted for. The precision of J1910’s mass measurement, as
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Companion ($M_C$) and pulsar ($M_P$) masses for a selection of redbacks (red) and black widows (black). Pre and post gravity darkening diffusion + convection models for J1910 are shown by the purple square and triangle, respectively. Spider mass demographics sourced from Strader et al. (2019), Nieder et al. (2020), Romani et al. (2021), Romani et al. (2022), Clark et al. (2021), Kennedy et al. (2022), and references therein.</figcaption></figure>

well as other irradiation-dominated spiders, is closely tied to our understanding of the irradiation in these systems (see Romani & Sanchez 2016; Sanchez & Romani 2017; Voisin et al. 2020; Zilles et al. 2020). In addition to full spectroscopy modelling, using high signal-to-noise spectra, and independent constraints would allow for a more precise mass measurement. Unfortunately here the inferred inclination is too low for a γ-ray eclipse, removing one independent constraint we might appeal to (Clark et al. 2023a).

## 4.2 Gravity darkening

Changing the gravity darkening prescription, as detailed in Section 3.1, has a notable effect on the inferred pulsar mass in J1910; a higher $M_P$ for pre-irradiation gravity darkening, and a lower one for post-irradiation. Masses in the system are not directly fitted for; they are derived from other parameters, and most specifically from $i$ and $K_2$. Given the high-precision binary mass function determined from the radio timing, the pulsar mass should roughly scale with the cube of the companion's centre-of-mass velocity and inversely with the cube of $\sin i$. As $i$ does not change significantly between the two prescriptions, $K_2$ must primarily drive the variation in pulsar mass. From the ratio $K_2$ between the two models, we would expect a ~25 per cent change in mass, while the actual difference is ~15 per cent. This implies that the changes cannot be entirely treated in isolation and that correlations between these two key parameters, and other ones from the model, contribute to dictating the masses.

Separately, we also observe that going from the pre-irradiation to the post-irradiation prescription causes the inferred values of $f_{RL}$, $T_{irr}$, and $q$ to decrease, and $T_{irr}$ to increase. Allowing for the irradiated face of the companion to be gravity darkened changes the balance between the irradiating flux and the star's size (mediated by $f_{RL}$). The

exact interplay between these parameters is difficult to disentangle and, while we cannot summarize it with a single effect, we can suggest a few correlations.

Changing the gravity darkening prescription naturally changes the heating pattern on the companion's surface. Temperature maps produced post-irradiation gravity darkening appears to shift heat, and thus flux, away from the centre of the irradiated face and towards the sides of the companion. This will shift the centre of light for any spectral lines, in our case the Mg β triplet, towards the centre of mass. Therefore, to match the observed line velocities, the sampled centre-of-mass $K_2$ must decrease to compensate. This effect is explored further in Section 4.3. This shifting of flux to the sides is likely linked to the smaller diffusion coefficient $\kappa$ found for the post-irradiation gravity darkening model.

$K_2$ directly constrains the mass ratio, which in turn changes the size of the companion's Roche lobe. Decreasing the companion's size lowers the overall flux we expect to receive. As $K_2$ has also decreased, the orbital separation must have also decreased to keep the period constant. A smaller separation and smaller companion mass would suggest the companion's Roche lobe become smaller. The filling factor must then reflect the size of the companion; to find both a lower filling factor and $K_2$ compared to the pre-irradiation models the companion must decrease in size. The nightside temperature remains similar for both approaches, so the lower flux expected from a smaller star on the nightside is compensated for by finding a lower distance.

The filling factor and $K_2$ (through the derived mass ratio) both affect the ellipsoidal component of the companion's optical variability. For example, a larger filling factor produces a more ellipsoidal star, adding flux at the orbital quadrature points ($\phi = 0.25, 0.75$). If the post-irradiation gravity darkening is moving flux from the centre to the sides of the companion, this in effect removes flux from the superior conjunction while adding it to the quadrature points, mimicking ellipsoidal modulation. This relieves the need for a large filling factor to reproduce the observed ellipsoidal component.

The irradiation efficiency, $\epsilon$, is also higher in the post-irradiation model, which is not surprising as heat is more effectively redistributed to the sides but the front of the star still needs to achieve the same temperature in order to reproduce the colours and amplitude at superior conjunction of the companion. For an irradiation-driven redback the irradiation component in the light curve must overcome the comparatively large ellipsoidal component, thus obtaining a high efficiency is not too surprising. Higher efficiencies have only previously been determined for PSR J1810+1744, an extremely irradiated black widow (Breton et al. 2012). Our pre-irradiation gravity darkening $\epsilon$ is comparable to that found for PSR J1555-2908 (Kennedy et al. 2022). However, much past ICARUS modelling has assumed a convective gravity darkening coefficient (0.08) which fundamentally affects the temperature on the companion's irradiated face. The stronger gravity darkening produced by the radiative coefficient deployed here requires more irradiation to achieve the same dayside temperature. In short, irradiation efficiencies of models with varying gravity darkening coefficients should not be directly compared. Post-irradiation gravity darkening then exacerbates this further, as the irradiation itself is gravity darkened. Yet more irradiating flux is then required to reproduce the temperature pattern. This quite naturally accounts for the increased $T_{irr}$ and $\epsilon$ for the post gravity darkening models.

Our modelling does not decisively indicate whether pre- or post-irradiation gravity darkening is preferred. Comparing our D+C mod-
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Companion ($M_C$) and pulsar ($M_P$) masses for a selection of redbacks (red) and black widows (black). Pre and post gravity darkening diffusion + convection models for J1910 are shown by the purple square and triangle, respectively. Spider mass demographics sourced from Strader et al. (2019), Nieder et al. (2020), Romani et al. (2021), Romani et al. (2022), Clark et al. (2021), Kennedy et al. (2022), and references therein.</figcaption></figure>

well as other irradiation-dominated spiders, is closely tied to our understanding of the irradiation in these systems (see Romani & Sanchez 2016; Sanchez & Romani 2017; Voisin et al. 2020; Zilles et al. 2020). In addition to full spectroscopy modelling, using high signal-to-noise spectra, and independent constraints would allow for a more precise mass measurement. Unfortunately here the inferred inclination is too low for a γ-ray eclipse, removing one independent constraint we might appeal to (Clark et al. 2023a).

## 4.2 Gravity darkening

Changing the gravity darkening prescription, as detailed in Section 3.1, has a notable effect on the inferred pulsar mass in J1910; a higher $M_P$ for pre-irradiation gravity darkening, and a lower one for post-irradiation. Masses in the system are not directly fitted for; they are derived from other parameters, and most specifically from $i$ and $K_2$. Given the high-precision binary mass function determined from the radio timing, the pulsar mass should roughly scale with the cube of the companion's centre-of-mass velocity and inversely with the cube of $\sin i$. As $i$ does not change significantly between the two prescriptions, $K_2$ must primarily drive the variation in pulsar mass. From the ratio $K_2$ between the two models, we would expect a ~25 per cent change in mass, while the actual difference is ~15 per cent. This implies that the changes cannot be entirely treated in isolation and that correlations between these two key parameters, and other ones from the model, contribute to dictating the masses.

Separately, we also observe that going from the pre-irradiation to the post-irradiation prescription causes the inferred values of $f_{RL}$, $T_{irr}$, and $q$ to decrease, and $T_{irr}$ to increase. Allowing for the irradiated face of the companion to be gravity darkened changes the balance between the irradiating flux and the star's size (mediated by $f_{RL}$). The

exact interplay between these parameters is difficult to disentangle and, while we cannot summarize it with a single effect, we can suggest a few correlations.

Changing the gravity darkening prescription naturally changes the heating pattern on the companion's surface. Temperature maps produced post-irradiation gravity darkening appears to shift heat, and thus flux, away from the centre of the irradiated face and towards the sides of the companion. This will shift the centre of light for any spectral lines, in our case the Mg β triplet, towards the centre of mass. Therefore, to match the observed line velocities, the sampled centre-of-mass $K_2$ must decrease to compensate. This effect is explored further in Section 4.3. This shifting of flux to the sides is likely linked to the smaller diffusion coefficient $\kappa$ found for the post-irradiation gravity darkening model.

$K_2$ directly constrains the mass ratio, which in turn changes the size of the companion's Roche lobe. Decreasing the companion's size lowers the overall flux we expect to receive. As $K_2$ has also decreased, the orbital separation must have also decreased to keep the period constant. A smaller separation and smaller companion mass would suggest the companion's Roche lobe become smaller. The filling factor must then reflect the size of the companion; to find both a lower filling factor and $K_2$ compared to the pre-irradiation models the companion must decrease in size. The nightside temperature remains similar for both approaches, so the lower flux expected from a smaller star on the nightside is compensated for by finding a lower distance.

The filling factor and $K_2$ (through the derived mass ratio) both affect the ellipsoidal component of the companion's optical variability. For example, a larger filling factor produces a more ellipsoidal star, adding flux at the orbital quadrature points ($\phi = 0.25, 0.75$). If the post-irradiation gravity darkening is moving flux from the centre to the sides of the companion, this in effect removes flux from the superior conjunction while adding it to the quadrature points, mimicking ellipsoidal modulation. This relieves the need for a large filling factor to reproduce the observed ellipsoidal component.

The irradiation efficiency, $\epsilon$, is also higher in the post-irradiation model, which is not surprising as heat is more effectively redistributed to the sides but the front of the star still needs to achieve the same temperature in order to reproduce the colours and amplitude at superior conjunction of the companion. For an irradiation-driven redback the irradiation component in the light curve must overcome the comparatively large ellipsoidal component, thus obtaining a high efficiency is not too surprising. Higher efficiencies have only previously been determined for PSR J1810+1744, an extremely irradiated black widow (Breton et al. 2012). Our pre-irradiation gravity darkening $\epsilon$ is comparable to that found for PSR J1555-2908 (Kennedy et al. 2022). However, much past ICARUS modelling has assumed a convective gravity darkening coefficient (0.08) which fundamentally affects the temperature on the companion's irradiated face. The stronger gravity darkening produced by the radiative coefficient deployed here requires more irradiation to achieve the same dayside temperature. In short, irradiation efficiencies of models with varying gravity darkening coefficients should not be directly compared. Post-irradiation gravity darkening then exacerbates this further, as the irradiation itself is gravity darkened. Yet more irradiating flux is then required to reproduce the temperature pattern. This quite naturally accounts for the increased $T_{irr}$ and $\epsilon$ for the post gravity darkening models.

Our modelling does not decisively indicate whether pre- or post-irradiation gravity darkening is preferred. Comparing our D+C mod-
els the Bayesian evidence as provided by the *dynesty* sampler is higher for the pre-irradiation gravity darkening case. The photometric fit is also better. However, post-irradiation gravity darkening models find a much tighter fit to the radial velocity curve. We tentatively support the post-irradiation gravity darkening case over the pre-irradiation gravity darkening due to the improved radial velocity fit in addition to our work as well as similar conclusions obtained by other authors (see Romani et al. 2021). This is also driven from the fact that it probably replicates the physical conditions on the companion's surface, though full scale simulations of an irradiated atmosphere would be required to settle this. In conclusion, we suggest that our post-irradiation gravity darkening D+C model is our 'best-fitting model' to characterize the companion in this system.

## 4.3 Centre-of-light corrections

As described in Section 1, surface heating of the companion is expected to affect where a given spectral line is emitted. Thus, a centre-of-light correction is needed to get the radial velocity determined for that line to reflect the true centre-of-mass radial velocity,

$$k = \frac{K_{\text{CoL}}}{K_{\text{CoM}}} \qquad (2)$$

Depending on where exactly the line is emitted, we should expect either a larger or smaller centre-of-light radial velocity than that the centre-of-mass radial velocity; larger if the line is preferentially emitted towards the nightside of the star (effectively orbiting at a larger radius than the CoM), or smaller if the line is stronger on the irradiated dayside. Linares, Shahbaz & Casares (2018) (hereafter L18) models PSR J2215+5135, as in this work, using Balmer-dominated and Mg β radial velocity curves. They calculate the expected equivalent width (EW) of each line across the companion's surface. They conclude the lower temperature Mg β line tracks the nightside and the high temperature Balmer series the dayside, 'bracketing' the K<sub>CoM</sub> between them.

Appendix A of Kandel & Romani (2020) adds some nuance to the 'bracketing' scenario. They assert that, while the EW of the Mg β triplet is indeed highest across the nightside, the raw EW is not the correct metric to use to measure the brightness of a given line. Rather, the EW must be weighted by the continuum flux at that point. A stronger line is not necessarily brighter, the local brightness dominates over the varying line strength over the surface. When weighting the EW by the local flux, the Mg β triplet is expected to be brightest towards the dayside, rejecting the 'bracketing' scenario.

Fig. 3 lends credence to the conclusion of Kandel & Romani (2020). The amplitude of our modelled radial velocity curve supports the Mg β feature being stronger towards the dayside, or at least does not support observing it towards the nightside, given it has a lower amplitude than the centre-of-mass velocity sampled to generate it. Table 5 displays the correction needed for the observed (red) curve to posterior K<sub>2</sub> for the pre- and post- irradiation gravity darkening D+C models. The exact value determined is clearly affected by the prescription chosen. Here we can appeal to our physical model. As in L18, we have calculated the EW of the H β and Mg β triplet across the companion's surface. To standardize our calculation we follow the procedure of Trager et al. (1998). Here, the flux-weighted EW (wEW) is calculated

<table><caption>**Table 5.** Centre-of-light corrections implied by pre- and post- irradiation gravity darkening D+C models.</caption><thead><tr><th>Gravity darkening</th><th>Centre-of-light correction</th></tr></thead><tbody><tr><td>Pre-irradiation</td><td>0.91 ± 0.05</td></tr><tr><td>Post-irradiation</td><td>0.98 ± 0.06</td></tr></tbody></table>

as

$$wEW = F_C \int_{\lambda_I}^{\lambda_H} \left(1 - \frac{F_\lambda}{F_C}\right) d\lambda, \qquad (3)$$

relative to a continuum level calculated either side of the spectral feature within pre-determined wavelength ranges, and weighted by the continuum level. The wEW for a given line can then be determined for every ICARUS surface element, producing an EW map of the surface.

Fig. 5 shows several absorption line surface maps produced for our D+C models, most notably the temperature and wEW. The temperature maps immediately reinforce differing heating patterns between the two options: applying gravity darkening after irradiation effectively removes flux from the centre of the dayside, while adding it to the sides of the companion as compared to the pre-irradiation gravity darkening case. The effect this has on the centre-of-light correction is then somewhat predictable. The broader irradiation of the post-irradiation model naturally reduces the correction needed, meaning the Mg β triplet more closely tracks the centre of mass. Conversely, the sharply heated dayside for the pre-irradiation gravity darkening case concentrates the line flux towards the companion's nose, exacerbating the correction needed.

Naturally the two line species can also be compared. For H β the wEW is clearly higher towards the dayside. The Mg β triplet is slightly stronger on the dayside, but relative to H β sees a fairly uniform distribution across the surface at all phases. This nicely reflects the expected interplay between the EW and continuum flux; For Mg β between the two distributions the whole surface is covered. By weighting the surface element velocities by their wEW we can make an estimate of the correction needed between the centre-of-mass and centre-of-light velocities. A physical interpretation of this is shown on the wEW map for each line: the red dashed line shows the effective centre-of-light position of the line relative to the centre of mass. For H β, matching the concentration of wEW on the dayside, the centre of light moves much closer to the nose of the star. For Mg β, we find the centre of light is actually nearly coincident with the centre of mass. Full surface plots including the EW and continuum flux maps are available in Fig. B3.

## 5 CONCLUSIONS

In this work we have presented the discovery, radio timing, and multiwavelength optical photometry of the redback PSR J1910-5320, as well as updated the radial velocity curve reported in Au et al. (2023). These data sets have been modelled using ICARUS, providing a new neutron star mass measurement. We have also tested our assumptions about the heating in spider systems, in particular examining whether the surface should be gravity darkened before or after the irradiation is applied to the companion.

Our modelling has constrained a number of system parameters. All our models find an inclination consistent with ~46°, and similar base temperatures consistent with our expectation from the SED. The remaining parameters vary bimodally, depending on whether gravity darkening is applied before or after irradiation. In particular the filling factor, irradiating temperature (and thus efficiency), companion
els the Bayesian evidence as provided by the *dynesty* sampler is higher for the pre-irradiation gravity darkening case. The photometric fit is also better. However, post-irradiation gravity darkening models find a much tighter fit to the radial velocity curve. We tentatively support the post-irradiation gravity darkening case over the pre-irradiation gravity darkening due to the improved radial velocity fit in addition to our work as well as similar conclusions obtained by other authors (see Romani et al. 2021). This is also driven from the fact that it probably replicates the physical conditions on the companion's surface, though full scale simulations of an irradiated atmosphere would be required to settle this. In conclusion, we suggest that our post-irradiation gravity darkening D+C model is our 'best-fitting model' to characterize the companion in this system.

## 4.3 Centre-of-light corrections

As described in Section 1, surface heating of the companion is expected to affect where a given spectral line is emitted. Thus, a centre-of-light correction is needed to get the radial velocity determined for that line to reflect the true centre-of-mass radial velocity,

$$k = \frac{K_{\text{CoL}}}{K_{\text{CoM}}} \qquad (2)$$

Depending on where exactly the line is emitted, we should expect either a larger or smaller centre-of-light radial velocity than that the centre-of-mass radial velocity; larger if the line is preferentially emitted towards the nightside of the star (effectively orbiting at a larger radius than the CoM), or smaller if the line is stronger on the irradiated dayside. Linares, Shahbaz & Casares (2018) (hereafter L18) models PSR J2215+5135, as in this work, using Balmer-dominated and Mg β radial velocity curves. They calculate the expected equivalent width (EW) of each line across the companion's surface. They conclude the lower temperature Mg β line tracks the nightside and the high temperature Balmer series the dayside, 'bracketing' the K<sub>CoM</sub> between them.

Appendix A of Kandel & Romani (2020) adds some nuance to the 'bracketing' scenario. They assert that, while the EW of the Mg β triplet is indeed highest across the nightside, the raw EW is not the correct metric to use to measure the brightness of a given line. Rather, the EW must be weighted by the continuum flux at that point. A stronger line is not necessarily brighter, the local brightness dominates over the varying line strength over the surface. When weighting the EW by the local flux, the Mg β triplet is expected to be brightest towards the dayside, rejecting the 'bracketing' scenario.

Fig. 3 lends credence to the conclusion of Kandel & Romani (2020). The amplitude of our modelled radial velocity curve supports the Mg β feature being stronger towards the dayside, or at least does not support observing it towards the nightside, given it has a lower amplitude than the centre-of-mass velocity sampled to generate it. Table 5 displays the correction needed for the observed (red) curve to posterior K<sub>2</sub> for the pre- and post- irradiation gravity darkening D+C models. The exact value determined is clearly affected by the prescription chosen. Here we can appeal to our physical model. As in L18, we have calculated the EW of the H β and Mg β triplet across the companion's surface. To standardize our calculation we follow the procedure of Trager et al. (1998). Here, the flux-weighted EW (wEW) is calculated

<table><caption>**Table 5.** Centre-of-light corrections implied by pre- and post- irradiation gravity darkening D+C models.</caption><thead><tr><th>Gravity darkening</th><th>Centre-of-light correction</th></tr></thead><tbody><tr><td>Pre-irradiation</td><td>0.91 ± 0.05</td></tr><tr><td>Post-irradiation</td><td>0.98 ± 0.06</td></tr></tbody></table>

as

$$wEW = F_C \int_{\lambda_I}^{\lambda_H} \left(1 - \frac{F_\lambda}{F_C}\right) d\lambda, \qquad (3)$$

relative to a continuum level calculated either side of the spectral feature within pre-determined wavelength ranges, and weighted by the continuum level. The wEW for a given line can then be determined for every ICARUS surface element, producing an EW map of the surface.

Fig. 5 shows several absorption line surface maps produced for our D+C models, most notably the temperature and wEW. The temperature maps immediately reinforce differing heating patterns between the two options: applying gravity darkening after irradiation effectively removes flux from the centre of the dayside, while adding it to the sides of the companion as compared to the pre-irradiation gravity darkening case. The effect this has on the centre-of-light correction is then somewhat predictable. The broader irradiation of the post-irradiation model naturally reduces the correction needed, meaning the Mg β triplet more closely tracks the centre of mass. Conversely, the sharply heated dayside for the pre-irradiation gravity darkening case concentrates the line flux towards the companion's nose, exacerbating the correction needed.

Naturally the two line species can also be compared. For H β the wEW is clearly higher towards the dayside. The Mg β triplet is slightly stronger on the dayside, but relative to H β sees a fairly uniform distribution across the surface at all phases. This nicely reflects the expected interplay between the EW and continuum flux; For Mg β between the two distributions the whole surface is covered. By weighting the surface element velocities by their wEW we can make an estimate of the correction needed between the centre-of-mass and centre-of-light velocities. A physical interpretation of this is shown on the wEW map for each line: the red dashed line shows the effective centre-of-light position of the line relative to the centre of mass. For H β, matching the concentration of wEW on the dayside, the centre of light moves much closer to the nose of the star. For Mg β, we find the centre of light is actually nearly coincident with the centre of mass. Full surface plots including the EW and continuum flux maps are available in Fig. B3.

## 5 CONCLUSIONS

In this work we have presented the discovery, radio timing, and multiwavelength optical photometry of the redback PSR J1910-5320, as well as updated the radial velocity curve reported in Au et al. (2023). These data sets have been modelled using ICARUS, providing a new neutron star mass measurement. We have also tested our assumptions about the heating in spider systems, in particular examining whether the surface should be gravity darkened before or after the irradiation is applied to the companion.

Our modelling has constrained a number of system parameters. All our models find an inclination consistent with ~46°, and similar base temperatures consistent with our expectation from the SED. The remaining parameters vary bimodally, depending on whether gravity darkening is applied before or after irradiation. In particular the filling factor, irradiating temperature (and thus efficiency), companion
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Surface maps for pre (top) and post (bottom) irradiation gravity darkening diffusion + convection models. The leftmost plot shows the surface temperature over the companion surface. The two plots on the right show the normalized flux weighted equivalent width (wEW) from each surface element (see equation 3). These are split into the Mg β triplet, which corresponds with our radial velocity curve, and the H β feature. In the picture of Linares, Shahbaz & Casares (2018) these track the companion nightside and dayside, respectively. The dashed lines on the wEW maps indicate the centre-of-mass (black) and centre-of-light (red) positions for the given line. Recall that a centre of light towards the companion's nose should correspond with a lower radial velocity determined for that line than the true centre-of-mass radial velocity (sampled by *ICARUS* and used to calculate *M*<sub>P</sub>).</figcaption></figure>

velocity, distance, and component masses change depending on
our gravity darkening prescription. For both models a moderate
pulsar mass is found, constrained to better than 15 per cent fractional
uncertainty at the 68 per cent level.

The novel radial velocity modelling deployed here has also
provided evidence that, as advanced in Kandel & Romani (2020), the
centre-of-light position of absorption species is not solely determined
by its activation temperature. We find the metallic, low temperature
Mg β triplet closely tracks the centre-of-mass velocity, balancing
the temperature dependence of the EW and continuum flux. This
is currently only verified for J1910, an irradiation-driven redback,
though our findings should also apply to other systems presenting
milder irradiation effects.

The modelling performed here aims to be widely applicable to all
spiders where photometry can be supplemented with radial velocity
curves. Further spider discovery and follow-up, particularly spectro-
scopic, is then desirable to provide more reliable measurement, taping
on better self-consistency in the way that the centre of mass is inferred
from spectral lines. While J1910 did not yield a 'supermassive'
neutron star, which can directly constrain the neutron star EoS, this
work adds to the tally of spider masses and can help understand
better the evolution landscape not only between black widows and
redbacks, but also across the other types of neutron star binaries.

ACKNOWLEDGEMENTS

RPB acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 715051; Spiders). OGD acknowledges the support of a Science and Technology Facilities

Council (STFC) stipend (Grant number: 2487578) to permit work as
a post graduate researcher. ECF is supported by NASA under award
number 80GSFC21M0002. VSD and ULTRACAM are funded by the
Science and Technology Facilities Council (grant ST/V000853/1).
JS acknowledges support by NASA grants 80NSSC22K1583 and
80NSSC23K1350, NSF grant AST-2205550, and the Packard Foun-
dation. CJC, EDB, PCCF, MK, and PVP acknowledge continuing
valuable support from the Max-Planck Society. KLL and KYA are
supported by the National Science and Technology Council of the
Republic of China (Taiwan) through grants 111-2636-M-006-024
and 112-2636-M-006-009. KLL is also a Yushan Young Fellow
supported by the Ministry of Education of the Republic of China
(Taiwan).

The MeerKAT telescope is operated by the South African Ra-
dio Astronomy Observatory, which is a facility of the National
Research Foundation, an agency of the Department of Science
and Innovation. We thank the staff at SARAO for their help with
observations and commissioning. TRAPUM observations used the
FBFUSE and APSUSE computing clusters for data acquisition,
storage, and analysis. These clusters were funded, designed, and
installed by the Max-Planck-Institut-für-Radioastronomie (MPIfR)
and the Max-Planck-Gesellschaft. FBFUSE performs beamforming
operations in real-time using the MOSAIC software stack (Chen
et al. 2021). Observations used the Pulsar Timing User Supplied
Equipment (PTUSE) servers at MeerKAT which were funded by
the MeerTime Collaboration members ASTRON, AUT, CSIRO,
ICRAR-Curtin, MPIfR, INAF, NRAO, Swinburne University of
Technology, the University of Oxford, UBC, and the University of
Manchester. The system design and integration was led by Swinburne
University of Technology and Auckland University of Technology
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Surface maps for pre (top) and post (bottom) irradiation gravity darkening diffusion + convection models. The leftmost plot shows the surface temperature over the companion surface. The two plots on the right show the normalized flux weighted equivalent width (wEW) from each surface element (see equation 3). These are split into the Mg β triplet, which corresponds with our radial velocity curve, and the H β feature. In the picture of Linares, Shahbaz & Casares (2018) these track the companion nightside and dayside, respectively. The dashed lines on the wEW maps indicate the centre-of-mass (black) and centre-of-light (red) positions for the given line. Recall that a centre of light towards the companion's nose should correspond with a lower radial velocity determined for that line than the true centre-of-mass radial velocity (sampled by *ICARUS* and used to calculate *M*<sub>P</sub>).</figcaption></figure>

velocity, distance, and component masses change depending on
our gravity darkening prescription. For both models a moderate
pulsar mass is found, constrained to better than 15 per cent fractional
uncertainty at the 68 per cent level.

The novel radial velocity modelling deployed here has also
provided evidence that, as advanced in Kandel & Romani (2020), the
centre-of-light position of absorption species is not solely determined
by its activation temperature. We find the metallic, low temperature
Mg β triplet closely tracks the centre-of-mass velocity, balancing
the temperature dependence of the EW and continuum flux. This
is currently only verified for J1910, an irradiation-driven redback,
though our findings should also apply to other systems presenting
milder irradiation effects.

The modelling performed here aims to be widely applicable to all
spiders where photometry can be supplemented with radial velocity
curves. Further spider discovery and follow-up, particularly spectro-
scopic, is then desirable to provide more reliable measurement, taping
on better self-consistency in the way that the centre of mass is inferred
from spectral lines. While J1910 did not yield a 'supermassive'
neutron star, which can directly constrain the neutron star EoS, this
work adds to the tally of spider masses and can help understand
better the evolution landscape not only between black widows and
redbacks, but also across the other types of neutron star binaries.

ACKNOWLEDGEMENTS

RPB acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 715051; Spiders). OGD acknowledges the support of a Science and Technology Facilities

Council (STFC) stipend (Grant number: 2487578) to permit work as
a post graduate researcher. ECF is supported by NASA under award
number 80GSFC21M0002. VSD and ULTRACAM are funded by the
Science and Technology Facilities Council (grant ST/V000853/1).
JS acknowledges support by NASA grants 80NSSC22K1583 and
80NSSC23K1350, NSF grant AST-2205550, and the Packard Foun-
dation. CJC, EDB, PCCF, MK, and PVP acknowledge continuing
valuable support from the Max-Planck Society. KLL and KYA are
supported by the National Science and Technology Council of the
Republic of China (Taiwan) through grants 111-2636-M-006-024
and 112-2636-M-006-009. KLL is also a Yushan Young Fellow
supported by the Ministry of Education of the Republic of China
(Taiwan).

The MeerKAT telescope is operated by the South African Ra-
dio Astronomy Observatory, which is a facility of the National
Research Foundation, an agency of the Department of Science
and Innovation. We thank the staff at SARAO for their help with
observations and commissioning. TRAPUM observations used the
FBFUSE and APSUSE computing clusters for data acquisition,
storage, and analysis. These clusters were funded, designed, and
installed by the Max-Planck-Institut-für-Radioastronomie (MPIfR)
and the Max-Planck-Gesellschaft. FBFUSE performs beamforming
operations in real-time using the MOSAIC software stack (Chen
et al. 2021). Observations used the Pulsar Timing User Supplied
Equipment (PTUSE) servers at MeerKAT which were funded by
the MeerTime Collaboration members ASTRON, AUT, CSIRO,
ICRAR-Curtin, MPIfR, INAF, NRAO, Swinburne University of
Technology, the University of Oxford, UBC, and the University of
Manchester. The system design and integration was led by Swinburne
University of Technology and Auckland University of Technology
in collaboration with SARAO and supported by the ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav) under grant CE170100004.

The Parkes 'Murriyang' radio telescope is part of the Australia Telescope National Facility (https://ror.org/05qajvd42) 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.

This work has used data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Based on observations collected at the European Southern Observatory under ULTRACAM GTO program 0109.D-0678. We are grateful to the astronomy and technical staff support at La Silla Observatory who enabled the remote observation program and observer John A. Paice.

Based on observations obtained at the Southern Astrophysical Research (SOAR) telescope, which is a joint project of the Ministério da Ciência, Tecnologia e Inovações (MCTI/LNA) do Brasil, the US National Science Foundation's NOIRLab, the University of North Carolina at Chapel Hill (UNC), and Michigan State University (MSU).

This research has used 'Aladin sky atlas' developed at CDS, Strasbourg Observatory, France (Bonnarel et al. 2000).

## DATA AVAILABILITY

TRAPUM and ULTRACAM observations are available upon reasonable request to the contact author. SOAR/Goodman observations are available from the NOIRLab Astro Data Archive https://astroarc.hive.noirlab.edu/.

## REFERENCES

* Abdollahi S. et al., 2020, ApJS, 247, 33

* Allard F., 2016, in Réylé C., Richard J., Cambrésy L., Deleuil M., Pérontal E., Tresse L., Vauglin I., eds, SF2A-2016: Proc. of the Annual meeting of the French Society of Astronomy and Astrophysics, Société Française d'Astronomie et d'Astrophysique, Lyon, p. 223

* Alpar M. A., Cheng A. F., Ruderman M. A., Shaham J., 1982, Nature, 300, 728

* Applegate J. H., 1992, ApJ, 385, 621

* Archibald A. M. et al., 2009, Science, 324, 1411

* Au K.-Y. et al., 2023, ApJ, 943, 103

* Bailes M. et al., 2020, Publ. Astron. Soc. Aust., 37, e028

* Bassa C. G. et al., 2014, MNRAS, 441, 1825

* Bezuidenhout M. C. et al., 2023, RAS Tech. Instrum., 2, 114

* Bhattacharya D., van den Heuvel E. P. J., 1991, Phys. Rep., 203, 1

* Bonnarel F. et al., 2000, A&AS, 143, 33

* Breton R. P., Rappaport S. A., van Kerkwijk M. H., Carter J. A., 2012, ApJ, 748, 115

* Breton R. P. et al., 2013, ApJ, 769, 108

* Brown A. J. et al., 2022, MNRAS, 513, 3050

* Castelli F., Kurucz R. L., 2003, in Piskunov N., Weiss W. W., Gray D. F., eds, Proc. IAU Symp. 210, Modelling of Stellar Atmospheres, Kluwer, Dordrecht, p. A20

* Chen W., Barr E., Karuppusamy R., Kramer M., Stappers B., 2021, J. Astron. Instrum., 10, 2150013

* Clark C. J. et al., 2021, MNRAS, 502, 915

* Clark C. J. et al., 2023a, Nat. Astron., 7, 451

* Clark C. J. et al., 2023b, MNRAS, 519, 5590

* Dhillon V. S. et al., 2007, MNRAS, 378, 825

* Dhillon V. S. et al., 2016, in Evans C. J., Simard L., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI. SPIE, Bellingham, p. 99080Y

* Dhillon V. et al., 2018, in Evans C. J., Simard L., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII. SPIE, Bellingham, p. 107020L

* Dhillon V. S. et al., 2021, MNRAS, 507, 350

* Djorgovski S., Evans C. R., 1988, ApJ, 335, L61

* Draghis P., Romani R. W., Filippenko A. V., Brink T. G., Zheng W., Halpern J. P., Camilo F., 2019, ApJ, 883, 108

* Drake A. J. et al., 2017, MNRAS, 469, 3688

* Espinosa Lara F., Rieutord M., 2012, A&A, 547, A32

* Feroz F., Hobson M. P., 2008, MNRAS, 384, 449

* Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601

* Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2019, Open J. Astrophys., 2, 10

* Folkner W. M., Williams J. G., Boggs D. H., 2009, Interplanetary Network Progress Report, 42-178, 1

* Fruchter A. S., Stinebring D. R., Taylor J. H., 1988, Nature, 333, 237

* Gaia Collaboration, 2016, A&A, 595, A1

* Gaia Collaboration, 2023, A&A, 674, A1

* Gullikson K., Dodson-Robinson S., Kraus A., 2014, AJ, 148, 53

* Hobbs G. et al., 2020, Publ. Astron. Soc. Aust., 37, e012

* Honeycutt R. K., 1992, PASP, 104, 435

* Jeffreys H., 1939, Theory of Probability. Oxford Univ. Press, Oxford

* Jiang L., Wang N., Chen W.-C., Li X.-D., Liu W.-M., Gao Z.-F., 2020, A&A, 633, A45

* Jonas J. L., 2009, IEEE Proc., 97, 1522

* Jonas J., the MeerKAT Team, 2016, in Taylor R., Camilo F., Leeuw L., Moodley K., eds, Proc. MeerKAT Science: On the Pathway to the SKA – PoS(MeerKAT2016). SISSA, Trieste, p. 001

* Kandel D., Romani R. W., 2020, ApJ, 892, 101

* Kandel D., Romani R. W., An H., 2021, ApJ, 917, L13

* Kennedy M. R. et al., 2022, MNRAS, 512, 3001

* Koposov S. E., 2019, Astrophysics Source Code Library, record ascl:1907.013

* Koposov S. E. et al., 2011, ApJ, 736, 146

* Lattimer J. M., 2012, Annu. Rev. Nucl. Part. Sci., 62, 485

* Levin L. et al., 2013, MNRAS, 434, 1387

* Li K.-L. et al., 2018, ApJ, 863, 194

* Linares M., 2019, Proc. Sci., Multifrequency Behaviour of High Energy Cosmic Sources - XIII. 3-8 June 2019. SISSA, Trieste, PoS#23

* Linares M., Shahbaz T., Casares J., 2018, ApJ, 859, 54 (L18)

* Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993

* Mata Sánchez D. et al., 2023, MNRAS, 520, 2217

* Nice D. et al., 2015, Astrophysics Source Code Library, record ascl:1509.002

* Nieder L. et al., 2020, ApJ, 902, L46

* Özel F., Freire P., 2016, ARA&A, 54, 401

* Papitto A. et al., 2013, Nature, 501, 517

* Polzin E. J., Breton R. P., Bhattacharyya B., Scholte D., Sobey C., Stappers B. W., 2020, MNRAS, 494, 2948

* Ray P. S. et al., 2013, ApJ, 763, L13

* Roberts M. S. E., 2013, in van Leeuwen J., ed., Proc. IAU Symp. 291, Neutron Stars and Pulsars: Challenges and Opportunities after 80 years. Kluwer, Dordrecht, p. 127

* Romani R. W., Sanchez N., 2016, ApJ, 828, 7

* Romani R. W., Kandel D., Filippenko A. V., Brink T. G., Zheng W., 2021, ApJ, 908, L46

* Romani R. W., Kandel D., Filippenko A. V., Brink T. G., Zheng W., 2022, ApJ, 934, L17

* Sanchez N., Romani R. W., 2017, ApJ, 845, 42

* Schlaflay E. F., Finkbeiner D. P., 2011, ApJ, 737, 103

* Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, AIP Conf. Proc. Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering. Am. Inst. Phys., New York, p. 395

* Smith D. A. et al., 2023, ApJ, 958, 191
in collaboration with SARAO and supported by the ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav) under grant CE170100004.

The Parkes 'Murriyang' radio telescope is part of the Australia Telescope National Facility (https://ror.org/05qajvd42) 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.

This work has used data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Based on observations collected at the European Southern Observatory under ULTRACAM GTO program 0109.D-0678. We are grateful to the astronomy and technical staff support at La Silla Observatory who enabled the remote observation program and observer John A. Paice.

Based on observations obtained at the Southern Astrophysical Research (SOAR) telescope, which is a joint project of the Ministério da Ciência, Tecnologia e Inovações (MCTI/LNA) do Brasil, the US National Science Foundation's NOIRLab, the University of North Carolina at Chapel Hill (UNC), and Michigan State University (MSU).

This research has used 'Aladin sky atlas' developed at CDS, Strasbourg Observatory, France (Bonnarel et al. 2000).

## DATA AVAILABILITY

TRAPUM and ULTRACAM observations are available upon reasonable request to the contact author. SOAR/Goodman observations are available from the NOIRLab Astro Data Archive https://astroarc.hive.noirlab.edu/.

## REFERENCES

Abdollahi S. et al., 2020, ApJS, 247, 33

Allard F., 2016, in Réylé C., Richard J., Cambrésy L., Deleuil M., Pécontal E., Tresse L., Vauglin I., eds, SF2A-2016: Proc. of the Annual meeting of the French Society of Astronomy and Astrophysics, Société Française d'Astronomie et d'Astrophysique, Lyon, p. 223

Alpar M. A., Cheng A. F., Ruderman M. A., Shaham J., 1982, Nature, 300, 728

Applegate J. H., 1992, ApJ, 385, 621

Archibald A. M. et al., 2009, Science, 324, 1411

Au K.-Y. et al., 2023, ApJ, 943, 103

Bailes M. et al., 2020, Publ. Astron. Soc. Aust., 37, e028

Bassa C. G. et al., 2014, MNRAS, 441, 1825

Beuzuidenhout M. C. et al., 2023, RAS Tech. Instrum., 2, 114

Bhattacharya D., van den Heuvel E. P. J., 1991, Phys. Rep., 203, 1

Bonnarel F. et al., 2000, A&AS, 143, 33

Breton R. P., Rappaport S. A., van Kerkwijk M. H., Carter J. A., 2012, ApJ, 748, 115

Breton R. P. et al., 2013, ApJ, 769, 108

Brown A. J. et al., 2022, MNRAS, 513, 3050

Castelli F., Kurucz R. L., 2003, in Piskunov N., Weiss W. W., Gray D. F., eds, Proc. IAU Symp. 210, Modelling of Stellar Atmospheres, Kluwer, Dordrecht, p. A20

Chen W., Barr E., Karuppusamy R., Kramer M., Stappers B., 2021, J. Astron. Instrum., 10, 2150013

Clark C. J. et al., 2021, MNRAS, 502, 915

Clark C. J. et al., 2023a, Nat. Astron., 7, 451

Clark C. J. et al., 2023b, MNRAS, 519, 5590

Dhillon V. S. et al., 2007, MNRAS, 378, 825

Dhillon V. S. et al., 2016, in Evans C. J., Simard L., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI. SPIE, Bellingham, p. 99080Y

Dhillon V. et al., 2018, in Evans C. J., Simard L., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII. SPIE, Bellingham, p. 107020L

Dhillon V. S. et al., 2021, MNRAS, 507, 350

Djorgovski S., Evans C. R., 1988, ApJ, 335, L61

Draghis P., Romani R. W., Filippenko A. V., Brink T. G., Zheng W., Halpern J. P., Camilo F., 2019, ApJ, 883, 108

Drake A. J. et al., 2017, MNRAS, 469, 3688

Espinosa Lara F., Rieutord M., 2012, A&A, 547, A32

Feroz F., Hobson M. P., 2008, MNRAS, 384, 449

Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601

Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2019, Open J. Astrophys., 2, 10

Folkner W. M., Williams J. G., Boggs D. H., 2009, Interplanetary Network Progress Report, 42-178, 1

Fruchter A. S., Stinebring D. R., Taylor J. H., 1988, Nature, 333, 237

Gaia Collaboration, 2016, A&A, 595, A1

Gaia Collaboration, 2023, A&A, 674, A1

Gullikson K., Dodson-Robinson S., Kraus A., 2014, AJ, 148, 53

Hobbs G. et al., 2020, Publ. Astron. Soc. Aust., 37, e012

Honeycutt R. K., 1992, PASP, 104, 435

Jeffreys H., 1939, Theory of Probability. Oxford Univ. Press, Oxford

Jiang L., Wang N., Chen W.-C., Li X.-D., Liu W.-M., Gao Z.-F., 2020, A&A, 633, A45

Jonas J. L., 2009, IEEE Proc., 97, 1522

Jonas J., the MeerKAT Team, 2016, in Taylor R., Camilo F., Leeuw L., Moodley K., eds, Proc. MeerKAT Science: On the Pathway to the SKA – PoS(MeerKAT2016). SISSA, Trieste, p. 001

Kandel D., Romani R. W., 2020, ApJ, 892, 101

Kandel D., Romani R. W., An H., 2021, ApJ, 917, L13

Kennedy M. R. et al., 2022, MNRAS, 512, 3001

Koposov S. E., 2019, Astrophysics Source Code Library, record ascl:1907.013

Koposov S. E. et al., 2011, ApJ, 736, 146

Lattimer J. M., 2012, Annu. Rev. Nucl. Part. Sci., 62, 485

Levin L. et al., 2013, MNRAS, 434, 1387

Li K.-L. et al., 2018, ApJ, 863, 194

Linares M., 2019, Proc. Sci., Multifrequency Behaviour of High Energy Cosmic Sources - XIII. 3-8 June 2019. SISSA, Trieste, PoS#23

Linares M., Shahbaz T., Casares J., 2018, ApJ, 859, 54 (L18)

Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993

Mata Sánchez D. et al., 2023, MNRAS, 520, 2217

Nice D. et al., 2015, Astrophysics Source Code Library, record ascl:1509.002

Nieder L. et al., 2020, ApJ, 902, L46

Özel F., Freire P., 2016, ARA&A, 54, 401

Papitto A. et al., 2013, Nature, 501, 517

Polzin E. J., Breton R. P., Bhattacharyya B., Scholte D., Sobey C., Stappers B. W., 2020, MNRAS, 494, 2948

Ray P. S. et al., 2013, ApJ, 763, L13

Roberts M. S. E., 2013, in van Leeuwen J., ed., Proc. IAU Symp. 291, Neutron Stars and Pulsars: Challenges and Opportunities after 80 years. Kluwer, Dordrecht, p. 127

Romani R. W., Sanchez N., 2016, ApJ, 828, 7

Romani R. W., Kandel D., Filippenko A. V., Brink T. G., Zheng W., 2021, ApJ, 908, L46

Romani R. W., Kandel D., Filippenko A. V., Brink T. G., Zheng W., 2022, ApJ, 934, L17

Sanchez N., Romani R. W., 2017, ApJ, 845, 42

Schlaflay E. F., Finkbeiner D. P., 2011, ApJ, 737, 103

Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, AIP Conf. Proc. Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering. Am. Inst. Phys., New York, p. 395

Smith D. A. et al., 2023, ApJ, 958, 191
Speagle J. S., 2020, MNRAS, 493, 3132

Stappers B., Kramer M., 2016, in Taylor R., Camilo F., Leeuw L., Moodley K., eds, Proc. MeerKAT Science: On the Pathway to the SKA – PoS(MeerKAT2016), Vol. 277. SISSA, Trieste, p. 009

Stappers B. W. et al., 1996, ApJ, 465, L119

Stappers B. W., van Kerkwijk M. H., Bell J. F., Kulkarni S. R., 2001, ApJ, 548, L183

Stappers B. W. et al., 2014, ApJ, 790, 39

Strader J. et al., 2015, ApJ, 804, L12

Strader J., Li K.-L., Chomiuk L., Heinke C. O., Udalski A., Peacock M., Shishkovsky L., Tremou E., 2016, ApJ, 831, 89

Strader J. et al., 2019, ApJ, 872, 42

Stringer J. G. et al., 2021, MNRAS, 507, 2174

Swihart S. J., Strader J., Chomiuk L., Aydi E., Sokolovsky K. V., Ray P. S., Kerr M., 2022, ApJ, 941, 199

Trager S. C., Worthey G., Faber S. M., Burstein D., González J. J., 1998, ApJS, 116, 1

Turchetta M., Linares M., Koljonen K., Sen B., 2023, MNRAS, 525, 2565

van Kerkwijk M. H., Breton R. P., Kulkarni S. R., 2011, ApJ, 728, 95

Voisin G., Kennedy M. R., Breton R. P., Clark C. J., Mata-Sánchez D., 2020, MNRAS, 499, 1758

Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29

Zilles A., Kotera K., Rohrmann R., Althaus L., 2020, MNRAS, 492, 1579

the depth of lines would change quite dramatically with metallicity. Thus, we can be relatively confident that radial velocities derived from a model can be reliable, even if some of the assumptions

<figure><img src="image_7.png" /><figcaption>**Figure A1.** Template (green) Mg *β* spectrum shifted (blue) to determine the effective radial velocity at various orbital phases (red).</figcaption></figure>

## APPENDIX A: RADIAL VELOCITY FITTING

The radial velocity fitting technique employed here fundamentally aims to take only the most essential information from *ICARUS* spectroscopy modelling. Comparing the model spectra with full observed spectra seems on the surface appealing as the fit can be informed both by the position and depth/profile of a set of lines. Not only is the radial velocity constrained but, in principle, also the temperature. However, systematic effects such as the exact elemental abundances can greatly complicate the situation and drive parameter estimation to compensate by modifying other parameters away from their ‘true’ values. Photometry modelling is not really affected by such considerations as line contribution to the total flux is negligible. Another important challenge to overcome is the considerable computational expense connected to the full modelling of a spectral data set.

The most essential, model constraining information to extract from a spectrum is the radial velocity, encoded in the Doppler shift of individual lines. This is highlighted particularly in the case of J1910, where we add a likelihood term according to the radial velocity curve rather than the observed spectroscopy directly. Determining radial velocities is, in theory, quite simple: the Doppler shift in a line’s wavelength relative to its value at rest reflects the velocity it was emitted at. The wavelength shift should be relatively insensitive to the systematics mentioned above if the overall line shape is not too dissimilar to the template which is being used. For example, we would assume that underpinning our model spectra with atmospheres of differing metallicities should not result in differing radial velocity measurements if we consider one line species at a time. Conversely,

regarding abundances are off so long as the temperature profile and stellar and binary parameters are captured adequately (via the photometry), Moreover, as we are only interested in individual lines the computational cost is greatly reduced.

Fig. A1 demonstrates our simplified spectroscopy modelling and radial velocity fit. Given a radial velocity curve, we generate a synthetic *ICARUS* spectrum for the orbital phases at which radial velocity measurements are available. A reference orbital phase is picked as a template – either that with the strongest line feature or closest to a user-defined phase. This template is then cross-correlated with the others for the wavelength, and thus velocity, shift. This produces a relative radial velocity curve within our model, with the expected sinusoidal shape. We then fit this to the observed curve, analytically minimizing a velocity offset, to find the additional likelihood term to the model (via a $\chi^2$ penalty). Even though the radial velocity measurements extracted from the observed spectra in Section 2.3 adopted a standard template profile, our model fitting to the velocity should closely resembles them for the reasons that were explained above.

## APPENDIX B: SUPPLEMENTARY PLOTS
Speagle J. S., 2020, MNRAS, 493, 3132

Stappers B., Kramer M., 2016, in Taylor R., Camilo F., Leeuw L., Moodley K., eds, Proc. MeerKAT Science: On the Pathway to the SKA – PoS(MeerKAT2016), Vol. 277. SISSA, Trieste, p. 009

Stappers B. W. et al., 1996, ApJ, 465, L119

Stappers B. W., van Kerkwijk M. H., Bell J. F., Kulkarni S. R., 2001, ApJ, 548, L183

Stappers B. W. et al., 2014, ApJ, 790, 39

Strader J. et al., 2015, ApJ, 804, L12

Strader J., Li K.-L., Chomiuk L., Heinke C. O., Udalski A., Peacock M., Shishkovsky L., Tremou E., 2016, ApJ, 831, 89

Strader J. et al., 2019, ApJ, 872, 42

Stringer J. G. et al., 2021, MNRAS, 507, 2174

Swihart S. J., Strader J., Chomiuk L., Aydi E., Sokolovsky K. V., Ray P. S., Kerr M., 2022, ApJ, 941, 199

Trager S. C., Worthey G., Faber S. M., Burstein D., González J. J., 1998, ApJS, 116, 1

Turchetta M., Linares M., Koljonen K., Sen B., 2023, MNRAS, 525, 2565

van Kerkwijk M. H., Breton R. P., Kulkarni S. R., 2011, ApJ, 728, 95

Voisin G., Kennedy M. R., Breton R. P., Clark C. J., Mata-Sánchez D., 2020, MNRAS, 499, 1758

Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29

Zilles A., Kotera K., Rohrmann R., Althaus L., 2020, MNRAS, 492, 1579

the depth of lines would change quite dramatically with metallicity. Thus, we can be relatively confident that radial velocities derived from a model can be reliable, even if some of the assumptions

<figure><img src="image_7.png" /><figcaption>**Figure A1.** Template (green) Mg *β* spectrum shifted (blue) to determine the effective radial velocity at various orbital phases (red).</figcaption></figure>

## APPENDIX A: RADIAL VELOCITY FITTING

The radial velocity fitting technique employed here fundamentally aims to take only the most essential information from *ICARUS* spectroscopy modelling. Comparing the model spectra with full observed spectra seems on the surface appealing as the fit can be informed both by the position and depth/profile of a set of lines. Not only is the radial velocity constrained but, in principle, also the temperature. However, systematic effects such as the exact elemental abundances can greatly complicate the situation and drive parameter estimation to compensate by modifying other parameters away from their ‘true’ values. Photometry modelling is not really affected by such considerations as line contribution to the total flux is negligible. Another important challenge to overcome is the considerable computational expense connected to the full modelling of a spectral data set.

The most essential, model constraining information to extract from a spectrum is the radial velocity, encoded in the Doppler shift of individual lines. This is highlighted particularly in the case of J1910, where we add a likelihood term according to the radial velocity curve rather than the observed spectroscopy directly. Determining radial velocities is, in theory, quite simple: the Doppler shift in a line’s wavelength relative to its value at rest reflects the velocity it was emitted at. The wavelength shift should be relatively insensitive to the systematics mentioned above if the overall line shape is not too dissimilar to the template which is being used. For example, we would assume that underpinning our model spectra with atmospheres of differing metallicities should not result in differing radial velocity measurements if we consider one line species at a time. Conversely,

regarding abundances are off so long as the temperature profile and stellar and binary parameters are captured adequately (via the photometry), Moreover, as we are only interested in individual lines the computational cost is greatly reduced.

Fig. A1 demonstrates our simplified spectroscopy modelling and radial velocity fit. Given a radial velocity curve, we generate a synthetic *ICARUS* spectrum for the orbital phases at which radial velocity measurements are available. A reference orbital phase is picked as a template – either that with the strongest line feature or closest to a user-defined phase. This template is then cross-correlated with the others for the wavelength, and thus velocity, shift. This produces a relative radial velocity curve within our model, with the expected sinusoidal shape. We then fit this to the observed curve, analytically minimizing a velocity offset, to find the additional likelihood term to the model (via a $\chi^2$ penalty). Even though the radial velocity measurements extracted from the observed spectra in Section 2.3 adopted a standard template profile, our model fitting to the velocity should closely resembles them for the reasons that were explained above.

## APPENDIX B: SUPPLEMENTARY PLOTS
<figure><img src="image_8.png" /><figcaption>**Figure B1.** Corner plot showing ICARUS fit parameters for pre- (red) and post- (blue) irradiation gravity darkening diffusion + convection models. Contours outline the 68, 95, and 99.7 per cent confidence intervals. The dashed lines on the 1D posterior plots show the 0.025, 0.5, and 0.975 quantiles.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure B1.** Corner plot showing ICARUS fit parameters for pre- (red) and post- (blue) irradiation gravity darkening diffusion + convection models. Contours outline the 68, 95, and 99.7 per cent confidence intervals. The dashed lines on the 1D posterior plots show the 0.025, 0.5, and 0.975 quantiles.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure B2.** Corner plot showing derived parameters for pre- (red) and post- (blue) irradiation gravity darkening diffusion + convection models. Contours outline the 68, 95, and 99.7 per cent confidence intervals. The dashed lines on the 1D posterior plots show the 0.025, 0.5, and 0.975 quantiles.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure B2.** Corner plot showing derived parameters for pre- (red) and post- (blue) irradiation gravity darkening diffusion + convection models. Contours outline the 68, 95, and 99.7 per cent confidence intervals. The dashed lines on the 1D posterior plots show the 0.025, 0.5, and 0.975 quantiles.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

**Figure B3.** Surface maps for pre- (top) and post- (bottom) irradiation gravity darkening diffusion + convection models. The leftmost plot shows the surface temperature over the companion surface. The grids on the right show the normalized flux weighted equivalent width (wEW), equivalent width (EW), and continuum flux from each surface element (see equation 3). These are split into the Mg β triplet, which corresponds with our radial velocity curve, and the H β feature. In the picture of Linares, Shahbaz & Casares (2018) these lines should track the companion nightside and dayside, respectively. The dashed lines on the wEW maps indicate the centre of mass (black) and centre of light (red) positions for the given line.

This paper has been typeset from a TeX/LATEX file prepared by the author.
<figure><img src="image_10.png" /></figure>

**Figure B3.** Surface maps for pre- (top) and post- (bottom) irradiation gravity darkening diffusion + convection models. The leftmost plot shows the surface temperature over the companion surface. The grids on the right show the normalized flux weighted equivalent width (wEW), equivalent width (EW), and continuum flux from each surface element (see equation 3). These are split into the Mg β triplet, which corresponds with our radial velocity curve, and the H β feature. In the picture of Linares, Shahbaz & Casares (2018) these lines should track the companion nightside and dayside, respectively. The dashed lines on the wEW maps indicate the centre of mass (black) and centre of light (red) positions for the given line.

This paper has been typeset from a TeX/LATEX file prepared by the author.