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

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 \textit{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 multi-color optical light curves obtained with ULTRACAM/NTT in combination with MeerKAT timing and updated SOAR/Goodman spectroscopic radial velocity measurements allow a mass constraint for PSR~J1910$-$5320. \texttt{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\pm0.2$ or $1.4\pm0.2$ $M_\odot$. The companion mass of either $0.45\pm0.04$ or $0.43^{+0.04}_{-0.03}$ $M_\odot$ also further confirms PSR~J1910$-$5320 as an irradiated redback spider pulsar.radiated redback spider pulsar.


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 ★ E-mail: oliver.dodge@manchester.ac.uk 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 onto 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).
Given the recycling scenario, spinning up an MSP requires a companion.However, since around 20% 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 et al. (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 (  < 24 hours) orbit.The companion is tidally locked to the pulsar, thus the irradiating pulsar wind heats one face whilst 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 et al. 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 a categories based on their companion mass: black widows with extremely low mass (  < 0.05 ⊙ ) and redbacks with higher companion masses (  ≳ 0.1 ⊙ ) (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 irradiaton, 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-mas X-ray binary (LMXB) 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 ⊙ (van Kerkwĳk et al. 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 characterised provides valuable mass measurements and a potential to constrain the EoS.Whilst there are a number of systematics and assumptions inherent to optical modelling 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 summarised as "centre-of-light" effects.Determining the binary mass ratio, requires to combine the well-measured pulsar's projected semi-major 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 center 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 Kerkwĳk et al. (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 et al. (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 9 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. 2015Strader et al. , 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 hour 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 si-nusoidal radial velocity curve confirmed the binary nature of the system with an orbital period   = 0.34847592 days.The observed radial velocity amplitude,  2,obs = 218 ± 8km 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 §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.§3 details the optical modelling of the optical light curves.In particular, we introduce a novel method to utilise 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.§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 §5.

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-minute observations of 79 sources from the 4FGL catalogue (Abdollahi et al. 2020), conducted using MeerKAT's -band receiver (at observing frequencies between 856-1712 MHz).This project has since been extended with a further two-pass survey (Thongmeearkom, T., et al., in prep.)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 localisation regions were too uncertain to cover in single observations at -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 ± 0.02 m s −2 indicative of a millisecond pulsar in a short-period binary system.We used SeeKAT 1 (Bezuidenhout et al. 2023) to localise this signal to a position less than 0.5 ′′ from an optical star detected in the Gaia DR3 (Gaia Collaboration et al. 2016Collaboration et al. , 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 hr 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 hr 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 64m telescope, during the half of the orbit centered 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 pseudologarithmically 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 8 observations were taken with MeerKAT's UHF receiver, the rest were performed at -band.A second pseudo-logarithmic timing campaign began with Parkes on 2022 September 06 until 2023 March 25.These observations each lasted 1.5 hr using the Ultrawide-band Low (UWL) 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 analyzed 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 et al. 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 on 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- 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.

Optical Photometry
We obtained multi-band light curves of J1910 on two nights, 2022 June 28th and 30th, using the ULTRACAM high-speed multi-band photometer (Dhillon et al. 2007), mounted on the 3.50m New Technology Telescope (NTT) at the European Southern Observatory (ESO) La Silla, Chile.The times and length of these observations are provided in Table 2. ULTRACAM utilises 3 CCDs simultaneously, each using a different Super Sloan Digital Sky Survey (Super-SDSS) u  g  r  i  z  filter (Dhillon et al. 2021).For these observations CCDs 1, 2 and 3 used the r  , g  and u  filters respectively.The data were taken under photometric conditions, with seeing varying between 1 -1.5 ′′ .The observations were reduced using the HiPERCAM (Dhillon  (Honeycutt 1992) was used to calibrate the r  and g  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 (Brown et al. 2022, Appendix A).Due to a lack of Gaia transform, and the unreliable transform between the HiPERCAM and SDSS filters, the u ′ band was calibrated by 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  and g  bands respectively (20s exposures), and 530 from the u  (60s exposures).Co-addition of u  band exposures, maximising S/N, leaves fewer u  datapoints 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  = 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.

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 zeropoint 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 miscentering 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 zeropoint corrections, we use the TelFit code (Gullikson et al. 2014) to generate a telluric absorption spectrum based on the airmass, the local humidity, pressure, and temperature, and the 3-hour 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 zeropoint 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 datasets show that the corrections from this method are generally similar to, but sometimes more accurate than, those from the night sky lines.
We also re-fit 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 §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 2 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 et al. 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.

PHOTOMETRIC MODELLING
The optical light curve modelling performed here utilised 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).

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 ,   , before irradiation was calculated as where  base is the Icarus input parameter specifying the temperature at the pole of the star,   is the surface gravity at surface element ,  pole is the surface gravity at the pole of the star and  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 () 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 toward 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 lowmass, 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 onto 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 preand 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 onto the companion's inner face, locked toward the pulsar.The parameters fit for this model constitute our fundamental set: the systemic velocity , the interstellar reddening E(B-V), the system inclination , the Roche-lobe filling factor  * RL 3 , the base and irradiating temperatures  base and  irr , the distance  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: , which parameterises the amplitude of the diffusion effect, and Γ, which governs the temperature dependence of the diffusion (Stringer et al. 2021).In this case, we have elected not to include Γ. 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 et al. 2018;Kandel & Romani 2020;Romani et al. 2021;Stringer et al. 2021), whereby light curves at not symmetric between the half orbits centered 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 parameterized by  amp .
(ii) A surface hot/cold spot with fitted temperature, size and position (e.g.Clark et al. 2021).
These models account for asymmetry by shifting or adding flux onto 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.Whilst 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 channeling 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 et al. 2009Feroz et al. , 2019)).Nested sampling algorithms provide the Bayesian evidence of a model, , a useful advantage over a classic Monte Carlo Markov Chain (MCMC) algorithms.Allowing for the 3 * Calculated as     1 , where   is the distance from the companion's barycentre to its nose, and   1 is the distance from the barycentre to the L1 point.
calculation of the Bayes factor, between two models enables one to determine which is favoured;  1,2 > 1 suggests model 1 is preferable, whereas  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  2 fit of the observed photometry and the simulated light curves generated given sampled parameters.

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), §3.3).The main priors used here were as follows: • A Gaussian prior applied to E(B-V), centered on the reddening provided by the dust maps of Schlafly & Finkbeiner (2011): 0.0596 ± 0.0033.
• A simple sin() prior applied to , corresponding to an isotropic distribution of orbital angular momentum vectors.
• 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 et al. 2023).Additional constraint can be provided by the DM inferred from radio timing using the Galactic electron density model Yao et al. (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 ± 0.49 kpc) is much smaller, and less reliable, than that from the Gaia parallax (6.8 ± 3.9 kpc).Yao et al. (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.

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 , 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 penalised according to the resulting likelihood.Kennedy et al. (2022) used a self-consistent procedure, 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 dataset.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.

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 (prevs post-irradiation and heat redistribution effects).In both heating models a consistent trend emerges: post-irradiation gravity darkening finds a smaller projected companion velocity  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 2 show the postirradiation 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  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  between the DH and D+C models, which implies a different inferred pulsar mass.Given a DH model will struggle to fit the amplitude of a asymmetric light curve it is unsurprising that , 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 toward 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 an notable departure from previously modelled redbacks.

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% 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 4th 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% 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% confidence regions, we find 4950 <   < 5100.This is slightly lower than our expectation from the broadband spectral energy distribution (SED) but within the allowed uncertainty (Au et al. 2023). 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 under-filling companion.Moving from the Icarus parameter  RL to the volume averaged filling factor we find an even smaller interval.Though significantly higher than the Icarus parameter  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.Figure 4 shows a collection of spider mass measurements, with the masses determined for our D+C models shown in purple; the square and triangle denoting 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 equation of state on their own.Linares (2019) collated a number of 'super-massive' 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, utilising post-Keplerian parameters measured through pulsar timing, out perform 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 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).

Gravity darkening
Changing the gravity darkening prescription, as detailed in §3.1, has a notable effect on the inferred pulsar mass in J1910; a higher   for pre-irradiation gravity darkening, and a lower one for postirradiation.Masses in the system are not directly fitted for; they are derived from other parameters, and most specifically from  and  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 .As  does not change significantly between the two prescriptions,  2 must primarily drive the variation in pulsar mass.From the ratio  2 between the two models, we would expect a ∼ 25% change in mass, while the actual difference is ∼ 15%.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  RL ,  irr and  to decrease, and  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  RL ).The exact interplay between these parameters is difficult to disentangle and, while we cannot summarise 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 toward the sides of the companion.This will shift the centre of light for any spectral lines, in our case the Mg triplet, toward the centre-of-mass.Therefore, to match the observed line velocities, the sampled centreof-mass  2 must decrease to compensate.This effect is explored further in §4.3.This shifting of flux to the sides is likely linked to the smaller diffusion coefficient  found for the post-irradiation gravity darkening model.
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  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 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  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 ( = 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 whilst 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, , 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  is comparable to that found for PSR J1555−2908 (Kennedy et al. 2022).However, much past Icarus modelling has assumed a convective gravity darkening coefficent (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.Postirradiation gravity darkening then exacerbates this further, as the irradiation itself is gravity darkened.Yet more irradiating flux is then PSR J1910−5320 11  Table 5. Centre-of-light corrections implied by pre-and post-irradiation gravity darkening D+C models.
Gravity Darkening Centre-of-light correction Pre-irradiation 1.05 ± 0.06 Post-irradiation 1.00 ± 0.07 required to reproduce the temperature pattern.This quite naturally accounts for the increased T irr and  for the post gravity darkening models.
Our modelling does not decisively indicate whether pre-or postirradiation gravity darkening is preferred.Comparing our D+C models 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-fit model' to characterise the companion in this system.

Centre-of-light corrections
As described in §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, .
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 et al. (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  CoM between them.
Appendix A of Kandel & Romani (2020) adds some nuance to the 'bracketing' scenario.They assert that, whilst 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.
Figure 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 toward the nightside, given it has a lower amplitude than the centre-of-mass velocity sampled to generate it.Table 4.3 displays the correction needed for the observed (red) curve.For both gravity darkening prescriptions, the correction is within 1 of coincidence with the centre of mass.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 as relative to a continuum level calculated either side of the spectral feature within predetermined wavelength ranges, and weighted by the continuum level.The wEW for a given line can then be determined for every Icarus surface element, producing a EW map of the surface.
Figure 5 shows several absorption line surface maps produced for our 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 center of the dayside, whilst 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 lowers 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-ofmass 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 appendix B3.

CONCLUSIONS
In this work we have presented the discovery, radio timing and multiwavelength optical photometry of the redback PSR J1910−5320, as well as updating the radial velocity curve reported in Au et al. (2023).These datasets 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 spectral energy distribution.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 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% fractional uncertainty at the 68% 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 spectroscopic, 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.Whilst J1910 did not yield a 'supermassive' neutron star, which can directly constrain the neutron star EoS, the current work adds to the tally of spider masses and can help understand better the evolution landscape between black widows and redbacks, but also across to other types of neutron star binaries.

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 dataset.
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, 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 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.
Figure 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 minimising a velocity offset, to find the additional likelihood term to the model (via a  2 penalty).Even though the radial velocity measurements extracted from the observed spectra in §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
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 2 .
Figure 2. Photometry fits produced by post-irradiation gravity darkening models presented in Table 4.The maximum a posteriori likelihood models have been selected.Left shows the direct heating model (DH), while right is diffusion + convection (D+C).The light curve data (Fig 1) are shown in the corresponding colours, 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.

Figure 3 .
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  2 for the best-fit model.Point wise residuals between the model and observed points are shown in the bottom panel.

Figure 5 .
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 normalised flux weighted equivalent width (wEW) from each surface element (see (3)).These are split into the Mg triplet, which corresponds with our radial velocity curve, and the H feature.In the picture of Linares et al. (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   ).

Figure A1 .
Figure A1.Template (green) Mg spectrum shifted (blue) to determine the effective radial velocity at various orbital phases (red).

Figure B1 .
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% confidence intervals.The dashed lines on the 1D posterior plots show the 0.025, 0.5 and 0.975 quantiles.

Figure B2 .
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% confidence intervals.The dashed lines on the 1D posterior plots show the 0.025, 0.5 and 0.975 quantiles.

Figure B3 .
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 normalised flux weighted equivalent width (wEW), equivalent width (EW) and continuum flux from each surface element (see (3)).These are split into the Mg triplet, which corresponds with our radial velocity curve, and the H feature.In the picture of Linares et al. (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.

Table 2 .
Time and phase coverage for ULTRACAM photometry obtained of J1910.The phase coverage, calculated with the timing ephemeris provided in Table1, corresponds with the phase axis of Figure1.The   filter is split from the other two due to the exclusion of irreducible data for the 28/06/2022 night.

Table 3 .
Updated radial velocities (RV) of PSR J1910−5320 from SOAR for both the full spectrum and targeting just the Mg triplet.