Properties of the redback millisecond pulsar binary 3FGL J0212.1+5320
T. Shahbaz M. Linares R. P. Breton
 
Hyperbolic Relaxation Method for Elliptic Equations
    Hannes R. Rüter^1, David Hilditch^1,2, Marcus
Bugner^1, and Bernd Brügmann^1
    December 30, 2023
================================================================================

 

Linares et al. (2016) obtained quasi-simultaneous 
g', r' and i'-band light curves and an absorption line radial velocity curve of the 
secondary star in the redback system . The light curves showed two maxima
and minima primarily due to the secondary star's ellipsoidal 
modulation, but with unequal maxima and minima.
We fit these light curves and radial velocities 
with our X-ray binary model including either a 
dark solar-type star spot or a hot spot due to 
off-centre heating from an intrabinary shock, to account for the unequal maxima. 
Both models give a radial velocity semi-ampltiude
and rotational broadening that agree with the observations. 
The observed secondary star's effective temperature is best matched 
with the value obtained using the hot spot model,
which gives a  neutron star 
and secondary star mass of M_ 1=1.85 ^+0.32_-0.26  and 
M_ 2=0.50 ^+0.22_-0.19 , respectively.  


binaries: close – 
stars: fundamental parameters – 
stars: individual: 3FGL J0212.1+5320 –  
stars: neutron – 
X-rays: binaries





§ INTRODUCTION


Progenitors of binary millisecond pulsars (MSPs) are believed to be 
recycled dead pulsars in low-mass X-ray binaries (LMXBs). According to 
this scenario the neutron star in the LMXB accretes material and angular 
momentum from its late-type companion star and is thereby spun up for 
billions of years to ultimately evolve into a MSP <cit.>. Compact binary MSPs with orbital periods less than 1 d are commonly classified as 
either ”black widows” or ”redbacks” depending on the mass of the 
companion star, M_ 2. Black widows have relatively low-mass degenerate
companion stars 
0.02 M_⊙≤ M_ 2≤ 0.05 M_⊙, whereas
redbacks have relatively more massive non-degenerate
companions 0.2 M_⊙≤ M_ 2≤ 0.4 M_⊙
<cit.>. 
The evolutionary link between redback and black widow MSPs is still 
uncertain. According to <cit.>, black widows and redbacks are two 
distinct populations of MSPs with different evaporation efficiencies. On 
the other hand, <cit.> argue that redbacks with compact 
orbits evolve to black widows, while the ones with longer orbital periods 
evolve to MSP-He white dwarf systems, and that  black widows descend 
from redbacks but not all redbacks become black widows.

The Large Area Telescope (LAT) on the Fermi Gamma-Ray Space Telescope has 
been successful in uncovering binary MSPs because they are γ-ray 
emitters similar to young pulsars <cit.>. Subsequently, thanks to 
targeted radio and X-ray surveys where Fermi-LAT 
has localized sources, more than 30 black widow and 14 redback systems 
have been discovered <cit.>. 
To date  there are three (PSR J1023+0038, IGR J18245-2452 and XSS J12270-4859) 
redback systems that transition between 
accretion-powered LMXB states and rotation-powered radio pulsar states 
clearly indicates the close relationship between LMXBs and radio MSPs and 
has provided further support for the recycling scenario 
<cit.>.

                            
The optical light curves of black widow and redback binaries show 
large-amplitude variability due to irradiated ellipsoidal modulation of 
the near-Roche filling secondary star. Compared to the black widows, the 
redbacks, with their larger secondary stars and closer distance are 
relatively bright in the optical, allowing for detailed photometric and 
spectroscopic observations <cit.>. 
In the last few years, dynamical photometric and spectroscopic studies of 
binary MSPs have largely proven the neutron star masses in these 
systems are generally heavier than the 1.4 M_⊙ canonical value, 
as expected theoretically for the binary evolution of MSPs 
<cit.>.

In this paper we present the results of modelling the optical light and 
radial velocity curves of the redback candidate  discovered by 
<cit.> and <cit.>. First we briefly describe the X-ray binary model used, 
the fitting procedure and the results where we determine the binary system 
masses. Finally we discuss the impact of these results.




§ 



Recently, both <cit.> and <cit.> presented the discovery 
of a variable optical counterpart to the unidentified γ-ray source 
, and argued that it is a binary "redback" MSP candidate. 
<cit.> obtained quasi-simultaneous g', r' and i'-band 
light curves as well a low- and high resolution spectroscopy. The optical 
light curves obtained by <cit.> and <cit.> both show two 
maxima and minima, primarily due to the secondary star's ellipsoidal 
modulation. However, the light curve obtained by <cit.> shows 
unequal maxima, which is not seen in the light curve obtained by 
<cit.> due to very poor orbital phase coverage. From the combined 
photometry and radial velocities 
<cit.> determined an orbital period of 0.86955(15) d, the same 
orbital period was obtained by <cit.> from their R and g'-band 
light curves.

From high resolution spectra taken between orbital phases 0.17 and 0.22, 
<cit.> found the secondary star to have a spectral type of  F6±2.
They also found no notable changes in spectral type or colour  across the 
binary orbit, not surprising given the long orbital period and hence weak  X-ray
heating effects. They determined a radial velocity curve using the  Hα
absorption line and from  a sinusoidal fit to the curve obtained a radial
velocity semi-amplitude of  K_ 2=214.1±5.0. Finally using high
resolution spectroscopy  they estimate the projected rotational velocity of the
secondary star to  be =73.2±1.6.













§ THE X-RAY BINARY MODEL




We use the X-ray binary light curve model 
xrbcurve described in <cit.>, which  has 
successfully been used to model the light curves and radial velocity 
curves of neutron star and black hole X-ray binaries 
<cit.>.
to fit and interpret the photometric light curves and the 
Hα absorption-line 
radial velocity curve of  presented in 
<cit.>. 
Briefly, the model consists of a binary system in which the primary  
point-like compact object with mass M_ 1 and the secondary 
is a 
Roche lobe filling star with mass M_ 2. 
which is assumed to 
be in a circular orbit and in synchronous rotation. The binary mass ratio 
q is defined as the ratio M_ 2/M_ 1.


The binary geometry  is determined by the 
binary masses, the orbital inclination i, and the Roche lobe 
filling factor f of the secondary star, defined as the ratio of the 
radius of the star from its centre of mass to the inner Lagrangian point 
to the distance towards the inner Lagrangian point. 




 §.§ The equivalent volume radius



When the secondary stars in interacting binaries such as cataclysmic 
variables or X-ray binaries are tidally locked and in synchronous rotation, 
for a given orbital period the width of the rotationally-broadened absorption
line profile 
arising from the star scales with the size of its Roche lobe. One can 
show that the star's rotational broadening  and equivalent volume 
radius <cit.> r_ eq/a (a is the binary separation), 
defined as the radius 
of a sphere whose volume is the same as the volume of the secondary star, 
are related through the expression


    v_ rot sin i / K_ 2 = (1+q) r_ eq(q)  / a.


For a star that fully fills its Roche lobe, <cit.> numerically 
calculated the Roche lobe volume for different mass ratios and determined 
an analytical expression for r_ eq/a as a function of q (normally 
referred to as Eggleton's formula). However, for stars that do not fill 
their Roche lobe, one cannot use Eggleton's formula to determine r_ 
eq/a, instead one has to calculate r_ eq/a numerically for a given 
binary q and f configuration. The star's rotational broadening is then 
given by


    v_ rot sin i / K_ 2 = (1+q)r_ eq(f,q)/a.


Using our X-ray binary model we numerically determine r_ eq/a as a 
function of f and q, assuming synchronous rotation. For a given q and 
f configuration we determine the binary Roche potential and thus the 
star's Roche lobe. We sample the star's surface with 18340 quadrilaterals 
of approximately equal area and then perform a numerical integration to 
calculate its volume and hence equivalent radius. Fig. <ref> 
shows a contour plot of the r_ eq/a values for different q and f 
values. In Fig. <ref> we show the results for f in the range 0.5 
to 1.0 and q in the range 0.1 to 1.0. The model f=1.0 is equivalent to 
Eggleton's relation and we find our model and Eggleton's relation are 
consistent to within 0.1 per cent.










 §.§ Light curve




The secondary star's effective temperature the gravity-darkening exponent 
and reddening given by  T_ 2,  β and E(B-V), respectively 
determines the observed light arising from the secondary star. 
The model includes the effects of heating 
of the secondary star by a point source from the 
compact object. 
The irradiating flux F_ X  is thermalised in the 
secondary star's photosphere and is re-radiated locally at a higher 
effective temperature. 
For each point on the star's surface, the effective temperature is calculated 
by combining the intrinsic and incident fluxes. 
We assume the "deep 
heating" approximation where the irradiation 
does not affect the temperature structure of the atmosphere. This means that
 each element radiates as predicted by a 
model atmosphere for a single star. 
The scale of the system is set by 
the distance to the source in parsecs, the orbital period and the radial velocity 
amplitude of the secondary star (D_ pc, P_ orb and K_ 2, respectively).
We use NextGen model-atmosphere fluxes <cit.> to 
determine the intensity distribution on the secondary star and 
a quadratic limb-darkening law with coefficients taken from <cit.>, to correct 
the intensity for limb-darkening.












 §.§ Radial velocity curve



The optical light curves of binary MSPs show the secondary star's 
ellipsoidal modulation combined with the effects of heating. 
The spectrum shows intrinsic absorption lines arising from the secondary 
star as well as irradiation-induced absorption Balmer lines. Both radial 
velocity curves are distorted because of shift on the centre of mass of 
the lines that are used to determine the radial velocities. 


The model light curve is determined by integrating the observed flux from each 
element of area on the star in the observers line-of-sight. For the radial velocity 
curves we specify the strength of the absorption or emission lines over 
the secondary star's surface, and integrate to obtain the corresponding 
line-of-sight radial velocity.  We set the 
absorption line strength given by its equivalent width (EW) according to 
the effective temperature for each element on the star.  However, as 
mentioned in <cit.> we must also consider the 
consequences of external heating. The vertical temperature gradient in an 
  internally and externally heated atmosphere  produces weaker absorption 
lines than expected from the effective temperature. Since there is no 
decent models to treat the effects of external heating in atmospheres, we 
introduce the factor F_ AV, which represent the fraction of the 
external radiation flux that exceeds the unperturbed flux. A value of F_ AV=1.10 means that if the external radiation flux is greater 
than 10 per cent of the unperturbed flux, then we set the EW for that 
element to zero, otherwise, the absorption line strength takes the EW 
corresponding to the effective temperature of the element.


The secondary stars in binary MSPs and X-ray binaries are typically 
late-type stars, later than F, which contain metals absorption 
lines such as magnesium, calcium and iron as well as neutral hydrogen. 
The 
strongest absorption metal lines in the blue part of the optical spectrum 
are the Ca I 4226Å, Fe I 4383Å and 
Mg I b triplet 5167,5172,5183Å. The EW versus temperature 
relation we use in our model is determined from observed stars in the 04 
to K2 spectral type range, obtained from the VLT-UVES Paranal Observatory 
Project database <cit.>.  In Fig.<ref> we show the EW 
versus T_ eff relation for the Hα absorption line as also for 
comparison the Mg I triplet. 
Given that we do not really understand irradiation, the EW relation 
combined with F_ AV allows the model 
to fit the observed line strength distribution 
<cit.>.




 §.§ Additional sources of light



Some of the optical light curves of binary MSPs are asymmetric with 
unequal maxima <cit.>. Various models have been 
proposed to explain this, such as off-centre heating from an intrabinary 
shock <cit.>, a hot spot <cit.> or solar-type star spots 
<cit.>.

The collision between the pulsar wind and the mass outflow from the 
secondary star can produce an intrabinary shock. Therefore the high energy 
radiation can be mediated by the intrabinary shock producing off-centre 
heating <cit.>. Indeed the changes in colour and  spectral type 
across the face of the secondary star in  does not match heating patterns,  
suggesting that the heating is not direct. 
X-ray observations of PSR J2215+5135 show an X-ray 
minimum near orbital phase zero (inferior conjunction of the secondary star) 
which has been interpreted as due to variable 
obscuration from an intrabinary shock around the secondary star 
<cit.>. In PSR J2215+5135, an intrabinary shock is suggested to 
explain the significant phase shift of the optical maximum with 
respect to the radio-pulse ephemeris, as well as the asymmetric optical 
light curves <cit.>. 


In some MSPs, the transition from an accretion-powered to a 
rotation-powered pulsar has been observed, with the disappearance of 
features associated with an accretion disk, for example PSR J1023+0038 
<cit.>. However, as accretion onto the pulsar ceases, it is 
possible that a quiescent disk could remain between the companion and the 
light cylinder of the pulsar <cit.>, which would contribute to the 
overall spectral energy distribution. In PSR J2215+5135, the difference 
between the observed colour temperature and the model temperature has been 
interpreted as an additional source of light from a hot accretion disk 
<cit.>. A large hot star spot could also produce an 
asymmetric temperature distribution.  It is also possible that the 
intrinsic magnetic field associated with the secondary star could channel 
the pulsar wind and cause enhanced local heating and an apparent hot spot 
on the star <cit.>.
Dark star spots due to strong magnetic activity on the secondary stars 
could also be present on the surface of the star <cit.>. Indeed 
Roche tomography has revealed observational evidence of star spots in 
binaries <cit.>.
Therefore, to account for the possible sources of extra light we include an 
additional flux components in each wave-band and to simulate the effects of a 
dark star spot or off-centre heating (hot spot) which produces an 
asymmetric temperature distribution on the star, we add a Gaussian 
temperature function to the elements of area on the star. The 
normalization is positive or negative for a hot or dark sport, 
respectively. The Gaussian position and width is constrained in latitude 
and is extended uniformly in longitude across the star, roughly simulating 
a spot.








§ XRBCURVE FITTING



Here we simultaneously fit the g',r',i'-band 
light curves and the absorption line radial velocity curve of  
presented in <cit.> with xrbcurve to determine the 
binary masses. The individual g' (380 data points), r' (629 data points), i' (382 data points) band data 
points, were phase folded according to the orbital ephemeris and averaged 
into 39, 32 and 39 orbital phase bins, respectively. Similarly, the 131 
radial velocity points were averaged into 30 phase bins.

The optical light curves of  clearly show unequal maxima. The fact 
that one sees a increase or decrease in light at phase 0.25 or 0.75, 
respectively, suggests that the extra source of light arises from the 
surface of the secondary star in terms of a dark star spot or a hot spot. 
To model these light curves with xrbcurve we assume either a hot 
or dark spot model. To reduce the number of free parameters in the model, 
we fix various model parameters. The value of the gravity-darkening 
exponent depends on the mean temperature and gravity of the secondary star 
<cit.>. The observed F6 late-type spectral type for the 
secondary star allows us to fix the gravity-darkening exponent to 0.072 
<cit.>. 
From the observed hydrogen column density we 
find the reddening to be E(B-V)=0.251±0.054 <cit.>.


The binary model parameters are f, q, cos i, T_ 2, log 
F_ X, D_ pc, K_ 2, F_ AV and E(B-V).
There is also the additional source 
of light in each filter E_ g', E_ r' and E_ i', the 
light curve phase shift δ_ g', δ_ r' and 
δ_ i', the radial velocity phase shift and systemic velocity 
δ_ AV, γ_ AV, and finally the Gaussian spot 
parameters, position C_ spot, width W_ spot, normalization 
T_ spot and extent E_ spot. 

Given q, cos i and K_ 
2 we calculate M_ 1 and M_ 2. Given q and f we 
determine R_ eq and which when combined with q and K_ 2 
gives  (see Section <ref>). However, note that the observed 
value for  depends on orbital phase, because it reflects the size 
of the tidally distorted star.
The value for R_ eq we determine represents the mean radius of the 
star observed and so the predicted value for   is a 
mean value.


To optimize the fitting procedure to 
first fit photometric light curve and the absorption-line radial 
velocity curve with xrbcurve assuming a dark spot 
 on the 
secondary star. We use the differential evolution algorithm described in 
<cit.> with a crossover probability of 0.9, 
a  mutation scaling factor of 0.9, and a 
 maximum number of generations of 2000, 
which corresponds to ∼ 100,000 calculations of the 
model. Given that there are a number of different types of data with 
different number of data points, to optimize the fitting procedure we 
assigned relative weights to the different data sets, determined after an 
initial fit. 
We rescaled the error bars on 
each data set so 
that the total reduced χ^2 of the fit was ∼ 1 for each data set 
separately.  

The simultaneous light and radial velocity curve fitting is a 
multi-dimensional nonlinear optimization problem. In order to obtain 
a robust error analysis 
we used Markov chain Monte Carlo (MCMC) method convolved with a differential 
evolution fitting algorithm; differential 
evolution adaptive metropolis (dream) 
<cit.>, which works well in problems with 
a high number of  dimensions. 
dream simultaneously runs multiple 
different chains uses differential evolution <cit.> as 
the genetic algorithm for population evolution with a Metropolis selection 
rule to decide whether parents are replaced by 
candidate points not. 
In order to ensure that the MCMC has converged 
we also performed the Gelman-Ruben test, 
which analyzes the difference between multiple 
Markov chains and we also visually inspected 
the trace of the parameters. 
The convergence is assessed by comparing the estimated 
between-chains and within-chain variances for each model parameter 
<cit.>.
We use a Bayesian framework to make statistical inferences on our binary 
model parameters. 
Our fitting makes use of flat prior probability 
distributions for all the model parameters, except E(B-V), where we use 
a Gaussian prior
[see <cit.> for an 
introduction to Bayesian analysis]. 
We use 20 individual chains to explore 
the parameter space and followed each for 40,000 iterations.  We reject 
the first 500 iterations ("burn-in") and only include every 10th point 
("thinning"), based on the observed auto correlation lengths for the 
individual parameters. 












§ RESULTS





We simultaneously fit g', r' and i'-band light curves and an 
absorption line radial velocity curve of  with xrbcurve 
with a constant source of light in each band to account for a possible 
disk and/or intrabinary shock and two model scenarios (a) a dark star spot or 
(b) a hot spot due to off-centre heating  to account for the unequal maxima.
Our preliminary fits show that the dark or hot spot model F_ AV is not 
constrained, which means that the model requires all the inner face of the 
star contributes to the radial velocity curve, so we fix it to a large value 
accordingly.
Furthermore, we find that the  phase shifts for each filter are the same, 
so in our final fits we set  
δ_ i'=δ_ r'=δ_ g'.

Also, for the hot spot model
we find that no heating is required, so we fix F_ X to zero. 
The dark and hot spot model have 18 and 17 free model parameters, respectively.
For our final fits, the dark and hot spot model best fits have a χ^2 
of 182 and 193, with 122 and 123 degrees of freedom,  respectively. 
Using the F-test to test the null hypothesis that the χ^2 of two 
models are the same, we find that we can reject the null hypothesis at the 
1.1 per cent significance level, 
Hence, statistically the dark spot 
model is only better than the 
hot spot model at the 98.9 per cent  (∼2.29σ) confidence level.


The best fit models for the hot and dark spot are 
shown in Figures <ref> and <ref>, respectively. The plots 
shows the best fit to the light and radial velocity curves as well as 
the light curve of the extra source of light in each band. We also show 
the observed maps of the continuum flux and Hαabsorption line strength 
on the secondary star at different orbital phases. In 
Table <ref> we give the mean and 1-σ limits of the 
posterior probability distribution functions on each parameter for both 
models. The result of the MCMC for the hot and dark spot model are shown 
in Figs. <ref> and <ref>, respectively.  As one can see, most 
of the two-dimensional  probability distribution functions are relatively 
well determined. 
Changing the gravity darkening coefficient only changes the χ^2 by  
3, resulting in a change in M_1 of 10 per cent.


The fits with the dark and hot 
spot model give similar values for f, i and K_ 2 (see 
Table <ref>). 
However, the value for q and T_ 2 are different.  The dark spot 
model best fit suggests a system at an inclination angle of 65^∘ and mass 
ratio q=0.18 with some effects of X-ray heating; 24 percent 
of the secondary star has a change in temperature more than 10 degrees 
larger than the non-irradiated level. The dark spot has a minumum 
temperature of T_ spot=332 K and covers 2
percent of the star; see Fig. <ref>). In contrast the best fit with a 
hot spot suggest a system at slightly higher inclination angle of 69^∘ and 
a less extreme mass ratio q=0.28 with no heating. 
The hot spot has a temperature of T_ spot=538 K 
and covers 3 percent of the star.
For the dark spot model we find that the extra source of light contributes 
13, 11 and 8  per cent to the observed flux in the g', r' and i' 
bands, respectively. Similarly for the hot spot model we find that the 
extra source of light contributes 31, 25 and 23  per cent to the observed 
flux in the g', r' and i' bands, respectively. In 
Fig. <ref> we show the spectrum of the extra flux component for 
the dark and hot spot model, which can be represented with a spectral index 
(F_ν∝ν^α) of 1.1± 0.3 and 0.8± 0.2, respectively.
     
In principle the shape of the ellipsoidal light curve can provide 
important clues to the binary inclination angle, because the peak to peak 
amplitude and difference between the two minima depend primarily on the 
binary inclination angle and mass ratio. Irradiation not only fills the 
minimum at phase 0.5, but it also a produces a larger amplitude 
modulation. 
However, for a Roche lobe under-filling star, this dependence 
with mass ratio disappears, because the star is not as tidally distorted 
compared to a fully Roche lobe filling star. For our values of f the 
dependence of the light curve amplitude with mass ratio still exists.
<cit.>. For a similar inclination angle, a system with a less extreme
mass ratio and no heating will produce a similar light curve to a 
system with a more extreme mass ratio and  with  heating. Our hot and dark 
spot model fits suggest either a non-irradiated system at an extreme mass ratio  
or weakly irradiated system at a less extreme mass ratio, respectively. 

The lack of heating of the secondary star 
required by the hot spot model by the pulsar wind and/or radiation is not 
surprising, given the wide orbit. The only other MSP which also has a wide 
orbit with a similar lack of irradiation is PSR J1740–5340 
<cit.>. 
Both the hot and dark spot model give similar values for 
f, K_ 2, i and , but different values for q and  
T_ 2. T_ 2 determined from the hot spot model 
is consistent with the observed 
spectral type F6±2 <cit.>, corresponding to a temperature 
range of 6170–6640 K <cit.>.









§ DISCUSSION




 §.§ The extra light source


<cit.> determined the spectral energy distribution of 
PSR J1023+0038 from near-IR to X-rays, when the system was in an 
accretion-powered phase. They found that the spectral energy distribution 
is well modelled by contributions from the secondary star, the accretion 
disc and an intrabinary shock. The neutron star spin-down luminosity 
irradiates the accretion disc and the secondary star, accounting for the 
UV and optical emission. X-rays and gamma-rays are produced in an 
intrabinary shock and the shock emission is powered by the neutron star 
spin-down luminosity which extends to lower energies with a photon index 
of Γ∼1.5 which corresponds to a spectral index of α∼2.5 (α = 1+Γ; EN(E) ∝ f_ν; N(E) ∝ 
E^Γ; F_ν∝ν^α, where N(E) is the X-ray (0.5–10 keV) 
photon spectrum 
and f_ν is the source spectrum). Indeed, the photon index measured for 
redbacks in the accretion or pulsar state is in the range 0.9–1.8 
<cit.> and  is no exception with a photon index of 
Γ∼ 1.3, corresponding to a spectral index of α∼ 2.3.

When determining  <cit.> found that their best match 
template star required a non-stellar light veiling of 10–30 per cent in 
the ∼ g'-band. For the dark spot model we find that the extra source 
of light contributes 13, 11 and 8 per cent to the observed flux in the g', 
r' and i' bands, respectively. Similarly for the hot spot model we 
find that the extra source of light contributes 31, 25 and 23 per cent to 
the observed flux in the g', r' and i' bands, respectively. Both models 
agree well with the veiling measured from the spectra <cit.>.
The extra light source contribution required to fit our light curves can be 
represented with a spectral index; F_ν∝ν^α. We find a spectral 
index of 1.1 and 0.8, for the dark and hot spot model, respectively, which does not
 agree with the X-ray spectral index, suggesting that the extra optical emission 
may not be produced by an intrabinary shock. 




 §.§ Energy flux


<cit.> and <cit.> have determined the spectral energy 
distribution of  which shows the optical band dominated by the 
secondary star, the gamma-rays from the MSP and the shock between the MSP 
and secondary star's wind (intrabinary shock), which most likely powers 
the X-ray emission. The unabsorbed 0.1–100 GeV 
energy flux is 1.71×10^-11, whereas the unabsorbed 
0.5–10 keV luminosity is 1.8×10^-12<cit.>. We can compare these values with the bolometric 
flux  determined from our dark spot model; no X-ray heating is 
required in the hot spot model.
For the  dark spot model we find a flux of  6.6×10^-11which 
is a factor of ∼4 more than what is observed in the γ-rays. 
The extra source of heating arise from the pulsar which 
emits prompt particles to heat the secondary star <cit.>




 §.§ The secondary stars rotational broadening



Synchronization via tidal forces will suppress differential rotation 
and make the angular velocity of the secondary star constant.
However, since the star is distorted, the linear rotational 
velocity will vary with longitude around the star and so  will vary 
across the orbit. The variations show two maxima and two minima per 
orbital cycle and have an amplitude of ∼10, which depends on q 
and i <cit.>. The procedure normally used to measure the 
secondary star's rotational broadening is to compare it with the spectrum 
of a slowly rotating template star convolved with a limb-darkened standard 
rotation profile <cit.>. 
The width of the 
limb-darkened standard rotation 
profile (adopting the continuum value for the limb-darkening coefficient) is varied until it matches the width of the target spectrum 
<cit.>. 
However, because of the 
assumed spherical shape of the rotation profile. 
there are assumptions inherent in this method.
The secondary stars in 
binary MSPs substantially fill their Roche lobes, and so they will have distorted 
line profiles which depend on the exact Roche binary geometry <cit.>. 
Also, because both temperature and gravity vary over the star's photosphere 
due to the Roche lobe shape, the spectrum of the secondary 
cannot be described by a single-star spectrum. 
Finally, we use the 
continuum value for the limb-darkening coefficient.
Firstly, because the line flux arises from higher regions in the 
atmosphere than the continuum flux, absorption lines have core 
limb-darkening coefficients much less than the continuum value <cit.>. 
Hence using the standard rotation profile with zero and the continuum value 
for the line 
limb-darkening coefficient gives a value for q that brackets the value 
found using the full geometrical treatment <cit.>. Using the 
extreme cases for the limb-darkening coefficient of zero or the continuum 
value introduces a systematic uncertainty of about 14 per cent in the 
determination  <cit.>.


<cit.> used the standard method described above to estimate the 
projected rotational velocity of the secondary star, which they find to be 
=73.2±1.6from spectra around orbital phase 0.2. They then 
use this  with their value for K_ 2, determined from a 
sinusoidal fit to their radial velocity curve and obtain q=0.26 (see 
equation  <ref>), assuming a Roche lobe filling, tidally locked and 
spherically symmetric companion star (see Section <ref>). 
However, the variation of  with orbital phase (10–15 per cent) 
introduces coupled with the uncertainties in the line limb-darkening 
coefficient (∼14 per cent) and the fact that the star does not fill 
its Roche lobe, contributes significantly to the accuracy to which one can 
determine q using this method.




 §.§ Favoured model



 is thought to be a redback based on its optical, 
X-ray and γ-ray properties, 
which are similar to other redbacks <cit.>. Most redbacks 
have secondary stars with minimum secondary star masses of <0.8<cit.> and so the secondary star masses determined from  
the dark and hot spot model are consistent with the system being a redback. 
 The reduced χ^2 of the best-fit for the dark and hot spot models 
are significantly different, with the dark spot model being only statistically 
better at the 98.9 (=2.29σ) confidence level (see Section <ref>), 
The value determined for  for both models agrees (within ∼3-σ) with
the observations (73.2±1.6  ), however, the hot
spot model with T_ 2=6416 K more constent with observations. 
(6170–6640 K <cit.>). We therefore, favour the
hot spot model which gives  a neutron star mass of 
M_ 1=1.85 ^+0.32_-0.26  and a secondary star mass of 
M_ 2=0.50 ^+0.22_-0.19 .







 §.§ Binary masses


<cit.> obtained R and g' band light curves and two absorption  line
radial velocity points. They used the ELC code  <cit.> to fit
the data, assuming T_ 2=5750 K, M_ 2∼0.4 M_ with no
irradiation effects. They performed a grid  search in i=60–90 degrees and
found β_V=0.70 to 0.64 (where  β_V is the ratio of the equivalent volume
radius of the companion star  to the Roche lobe radius) corresponding to M_
1=1.5 to  2.2 M_. Note that their definition for the Roche lobe
filling  factor is different to that in our model. Their estimated range for 
β_V corresponds to f=0.50 to 0.54 (determined using xrbcurve), 
which is very different to the value we obtain (see  Table <ref>).
However, as noted by the authors, given their  poor fits and large systematics,
the binary parameters they obtain are  merely indicative. In contrast we simultaneously
fit the absorption line radial  velocity and g', r', i' light curve, which
allow us to constrain various model parameters such as the filling factor of the
companion and its surface temperature. It is evident that this redback companion
does not fill its Roche lobe, thus confirming the previous finding  by
<cit.> that redbacks and black widows may have smaller radii than
originally expected (due to the presence of radio eclipses in some systems).

The presence of temperature asymmetry also indicates that a hot or cold spot is
present on the companion's surface. The physical reason for this remains
unclear  (see Section <ref> for possible models), but long-term
monitoring could provide further insights to shed light on this mystery. A cold
or a hot spot due to an atmospheric phenomenon might move on a time scale
related to the magnetic activity on the star for instance.  Star spots are 
believed to be common features on many late-type stars
and typically vary on timescales of months to a few years <cit.>. 
Star spots have been observed in the surface map of the secondary star in the
LMXB Cen X–4  <cit.>. Indeed, in cataclysmic variables, studies of
the surface maps of the secondary stars show star spots moving on timescale of
few days  <cit.>.

There are about 32 reliable neutron star mass measurements in X-ray binaries, 
double neutron star systems and millisecond pulsars, most of which are 
are at least partly based on the radio timing techniques <cit.>,
Studies statistically distinguish  between different types of neutron stars and 
 between those believed to be close to 
their birth masses, with a mean mass of 1.35±0.05  which reflects a highly tuned 
formation channel <cit.>, and the ones that have undergone  
long-term accretion episodes.
Most recycled pulsars are accompanied by mid-F to late-K type low-mass 
<0.5 companions but for these systems to form in the age of the 
Universe, the companion must originally have been a star of ∼1. 
Thus, the companion star must have lost a large fraction of 
its mass, accreted onto the neutron star, leading to a increase in the neutron 
star's mass. The most precise massive neutron star measurements have been 
obtained through the measurement of the Shapiro delay, 1.94±0.04 and 
1.667±0.021 in PSR J1614-2230 and PSR J1903+0327, respectively 
<cit.>. These masses has allowed us to place fundamental constraints 
on the equation of state of nuclear matter at high densities, excluding 
many of the soft equations of state <cit.>. Dynamical photometric 
and spectroscopic studies of binary MSPs have largely revealed that 
neutron star masses in these systems are generally heavier than the 
1.4 M_⊙ canonical value, 
with some particularly heavy ones in PSR B1957+20 (M_ 1=2.4, 
) and PSR J1816+4510 (M_ 1≥1.84, 
). 
Our mass measurement for the neutron star in  using  the 
hot spot model suggests a massive 1.85 ^+0.32_-0.26 neutron star. 




§ ACKNOWLEDGEMENTS


This research has been supported by the Spanish Ministry of Economy and 
Competitiveness (MINECO) under the grant AYA2013-42627.
M.L. is supported by EU's Horizon 2020 programme through a 
Marie Sklodowska-Curie Fellowship (grant nr. 702638). R.P.B. received
funding  from  the  European  Union  Seventh  Framework
Programme under grant agreement PIIF-GA-2012-332393.
This paper makes use of the IAC's Supercomputing facility
condor.



 

@urlcharsothermakeother $&#_%
@doi@urlcharsother ifnextchar [ @doi@
  @doi@[]
@doi@[#1]#2tempa#1tempaempty http://dx.doi.org/#2 doi:#2http://dx.doi.org/#2 #1
@eprint#1#2@eprint@#1:#2::nil
@eprint@arXiv#1http://arxiv.org/abs/#1 arXiv:#1
@eprint@dblp#1http://dblp.uni-trier.de/rec/bibtex/#1.xml
  dblp:#1
@eprint@#1:#2:#3:#4niltempa #1tempb #2tempc
  #3tempc empty tempc tempb tempb tempa tempb empty tempb arXivifundefined
  mn@eprint@tempbtempb:tempcmn@eprint@tempbtempc

[AbdoAbdo2009]Abdo09
Abdo A. A. e. a.,  2009, @doi [Science] 10.1126/science.1176113, http://adsabs.harvard.edu/abs/2009Sci...325..848A 325, 848

[AceroAcero2015]Acero15
Acero F. e. a.,  2015, @doi [] 10.1088/0067-0049/218/2/23, http://adsabs.harvard.edu/abs/2015ApJS..218...23A 218, 23

[Alpar, Cheng, Ruderman  &
  ShahamAlpar et al.1982]Alpar82
Alpar M. A.,  Cheng A. F.,  Ruderman M. A.,   Shaham J.,  1982, @doi
  [] 10.1038/300728a0, http://adsabs.harvard.edu/abs/1982Natur.300..728A 300, 728

[Archibald et al.,Archibald
  et al.2009]Archibald09
Archibald A. M.,  et al., 2009, @doi [Science] 10.1126/science.1172740,
  http://adsabs.harvard.edu/abs/2009Sci...324.1411A 324, 1411

[Bagnulo, Jehin, Ledoux, Cabanac,
  Melo, Gilmozzi  & ESO Paranal Science Operations TeamBagnulo
  et al.2003]Bagnulo03
Bagnulo S.,  Jehin E.,  Ledoux C.,  Cabanac R.,  Melo C.,  Gilmozzi
  R.,   ESO Paranal Science Operations Team 2003, The Messenger, http://adsabs.harvard.edu/abs/2003Msngr.114...10B 114, 10

[Bassa et al.,Bassa
  et al.2014]Bassa14
Bassa C. G.,  et al., 2014, @doi [] 10.1093/mnras/stu708, http://adsabs.harvard.edu/abs/2014MNRAS.441.1825B 441, 1825

[Benvenuto, De Vito  &
  HorvathBenvenuto et al.2014]Benvenuto14
Benvenuto O. G.,  De Vito M. A.,   Horvath J. E.,  2014, @doi []
  10.1088/2041-8205/786/1/L7, http://adsabs.harvard.edu/abs/2014ApJ...786L...7B 786, L7

[Bhattacharya & van den
  HeuvelBhattacharya & van den Heuvel1991]Bhattacharya91
Bhattacharya D.,  van den Heuvel E. P. J.,  1991, @doi []
  10.1016/0370-1573(91)90064-S, http://adsabs.harvard.edu/abs/1991PhR...203....1B 203, 1

[Bouvier & BertoutBouvier &
  Bertout1989]Bouvier89
Bouvier J.,  Bertout C.,  1989, , http://adsabs.harvard.edu/abs/1989A

[Breton et al.,Breton
  et al.2013]Breton13
Breton R. P.,  et al., 2013, @doi [] 10.1088/0004-637X/769/2/108,
  http://adsabs.harvard.edu/abs/2013ApJ...769..108B 769, 108

[Chen, Chen, Tauris  & HanChen
  et al.2013]Chen13
Chen H.-L.,  Chen X.,  Tauris T. M.,   Han Z.,  2013, @doi []
  10.1088/0004-637X/775/1/27, http://adsabs.harvard.edu/abs/2013ApJ...775...27C 775, 27

[Claret & BloemenClaret &
  Bloemen2011]Claret11
Claret A.,  Bloemen S.,  2011, @doi []
  10.1051/0004-6361/201116451, http://cdsads.u-strasbg.fr/abs/2011A

[Collins & TruaxCollins &
  Truax1995]Collins95
Collins II G. W.,  Truax R. J.,  1995, @doi [] 10.1086/175225,
  http://adsabs.harvard.edu/abs/1995ApJ...439..860C 439, 860

[Coti Zelati et al.,Coti Zelati
  et al.2014]Zelati14
Coti Zelati F.,  et al., 2014, @doi [] 10.1093/mnras/stu1552,
  http://adsabs.harvard.edu/abs/2014MNRAS.444.1783C 444, 1783

[Crawford et al.,Crawford
  et al.2013]Crawford13
Crawford F.,  et al., 2013, @doi [] 10.1088/0004-637X/776/1/20,
  http://adsabs.harvard.edu/abs/2013ApJ...776...20C 776, 20

[Demorest, Pennucci, Ransom, Roberts
  & HesselsDemorest et al.2010]Demorest10
Demorest P. B.,  Pennucci T.,  Ransom S. M.,  Roberts M. S. E.,
  Hessels J. W. T.,  2010, @doi [] 10.1038/nature09466, http://adsabs.harvard.edu/abs/2010Natur.467.1081D 467, 1081

[Deneva et al.,Deneva
  et al.2016]Deneva16
Deneva J. S.,  et al., 2016, @doi [] 10.3847/0004-637X/823/2/105,
  http://adsabs.harvard.edu/abs/2016ApJ...823..105D 823, 105

[EggletonEggleton1983]Eggleton83
Eggleton P. P.,  1983, @doi [] 10.1086/160960, http://adsabs.harvard.edu/abs/1983ApJ...268..368E 268, 368

[Ekşİ & AlparEkşİ &
  Alpar2005]Eksi05
Ekşİ K. Y.,  Alpar M. A.,  2005, @doi [] 10.1086/425959,
  http://adsabs.harvard.edu/abs/2005ApJ...620..390E 620, 390

[Freire et al.,Freire
  et al.2011]Freire11
Freire P. C. C.,  et al., 2011, @doi []
  10.1111/j.1365-2966.2010.18109.x, http://adsabs.harvard.edu/abs/2011MNRAS.412.2763F 412, 2763

[Gelman & RubinGelman &
  Rubin1992]Gelman92
Gelman A.,  Rubin D.,  1992, Statistical Science, 7, 457

[Gentile et al.,Gentile
  et al.2014]Gentile14
Gentile P. A.,  et al., 2014, @doi [] 10.1088/0004-637X/783/2/69,
  http://adsabs.harvard.edu/abs/2014ApJ...783...69G 783, 69

[GrayGray2005]Gray05
Gray D. F.,  2005, The Observation and Analysis of Stellar Photospheres:, 3
  edn.
Cambridge University Press, Cambridge, @doi10.1017/CBO9781316036570

[GregoryGregory2005]Gregory05
Gregory P. C.,  2005, Bayesian Logical Data Analysis for the Physical
  Sciences: A Comparative Approach with `Mathematica' Support.
Cambridge University Press

[Hauschildt, Allard  &
  BaronHauschildt et al.1999]Hauschildt99
Hauschildt P. H.,  Allard F.,   Baron E.,  1999, @doi []
  10.1086/306745, http://adsabs.harvard.edu/abs/1999ApJ...512..377H
  512, 377

[Hessels et al.,Hessels
  et al.2011]Hessels11
Hessels J. W. T.,  et al., 2011, in Burgay M.,  D'Amico N.,  Esposito
  P.,  Pellizzoni A.,   Possenti A.,  eds,  American Institute of Physics
  Conference Series Vol. 1357, American Institute of Physics Conference Series.
  pp 40–43 (@eprint arXiv 1101.1742), @doi10.1063/1.3615072

[Hill, Watson, Shahbaz, Steeghs  &
  DhillonHill et al.2014]Hill14
Hill C. A.,  Watson C. A.,  Shahbaz T.,  Steeghs D.,   Dhillon V. S.,
   2014, @doi [] 10.1093/mnras/stu1460, http://adsabs.harvard.edu/abs/2014MNRAS.444..192H 444, 192

[Kaplan, Bhalerao, van Kerkwijk,
  Koester, Kulkarni  & StovallKaplan et al.2013]Kaplan13
Kaplan D. L.,  Bhalerao V. B.,  van Kerkwijk M. H.,  Koester D.,
  Kulkarni S. R.,   Stovall K.,  2013, @doi []
  10.1088/0004-637X/765/2/158, http://adsabs.harvard.edu/abs/2013ApJ...765..158K 765, 158

[Kiziltan, Kottas, De Yoreo  &
  ThorsettKiziltan et al.2013]Kiziltan13
Kiziltan B.,  Kottas A.,  De Yoreo M.,   Thorsett S. E.,  2013, @doi
  [] 10.1088/0004-637X/778/1/66, http://adsabs.harvard.edu/abs/2013ApJ...778...66K 778, 66

[Lattimer & PrakashLattimer &
  Prakash2004]Lattimer04
Lattimer J. M.,  Prakash M.,  2004, @doi [Science]
  10.1126/science.1090720, http://adsabs.harvard.edu/abs/2004Sci...304..536L 304, 536

[Li, Halpern  & ThorstensenLi
  et al.2014]Li14
Li M.,  Halpern J. P.,   Thorstensen J. R.,  2014, @doi []
  10.1088/0004-637X/795/2/115, http://adsabs.harvard.edu/abs/2014ApJ...795..115L 795, 115

[Li, Kong, Hou, Mao, Strader,
  Chomiuk  & TremouLi et al.2016]Li16
Li K.-L.,  Kong A. K. H.,  Hou X.,  Mao J.,  Strader J.,  Chomiuk
  L.,   Tremou E.,  2016, @doi [] 10.3847/1538-4357/833/2/143, http://adsabs.harvard.edu/abs/2016ApJ...833..143L 833, 143

[LinaresLinares2014]Linares14
Linares M.,  2014, @doi [] 10.1088/0004-637X/795/1/72, http://adsabs.harvard.edu/abs/2014ApJ...795...72L 795, 72

[Linares, Miles-Páez,
  Rodríguez-Gil, Shahbaz, Casares, Fariña  &
  KarjalainenLinares et al.2016]Linares16
Linares M.,  Miles-Páez P.,  Rodríguez-Gil P.,  Shahbaz T.,
  Casares J.,  Fariña C.,   Karjalainen R.,  2016, preprint, http://adsabs.harvard.edu/abs/2016arXiv160902232L  (@eprint arXiv
  1609.02232)

[LucyLucy1967]Lucy67
Lucy L. B.,  1967, , http://adsabs.harvard.edu/abs/1967ZA.....65...89L 65, 89

[Marsh, Robinson  & WoodMarsh
  et al.1994]Marsh94
Marsh T. R.,  Robinson E. L.,   Wood J. H.,  1994, @doi []
  10.1093/mnras/266.1.137, http://adsabs.harvard.edu/abs/1994MNRAS.266..137M 266, 137

[Orosz & HauschildtOrosz &
  Hauschildt2000]Orosz00
Orosz J. A.,  Hauschildt P. H.,  2000, , http://adsabs.harvard.edu/abs/2000A

[Orosz & van KerkwijkOrosz & van
  Kerkwijk2003]Orosz03
Orosz J. A.,  van Kerkwijk M. H.,  2003, @doi []
  10.1051/0004-6361:20021468, http://adsabs.harvard.edu/abs/2003A

[Özel & FreireÖzel &
  Freire2016]Ozel16
Özel F.,  Freire P.,  2016, @doi []
  10.1146/annurev-astro-081915-023322, http://adsabs.harvard.edu/abs/2016ARA

[Özel, Psaltis, Narayan  & Santos
  VillarrealÖzel et al.2012]Ozel12
Özel F.,  Psaltis D.,  Narayan R.,   Santos Villarreal A.,  2012,
  @doi [] 10.1088/0004-637X/757/1/55, http://adsabs.harvard.edu/abs/2012ApJ...757...55O 757, 55

[Papitto et al.,Papitto
  et al.2013]Papitto13
Papitto A.,  et al., 2013, @doi [] 10.1038/nature12470, http://adsabs.harvard.edu/abs/2013Natur.501..517P 501, 517

[Pecaut & MamajekPecaut &
  Mamajek2013]Pecaut13
Pecaut M. J.,  Mamajek E. E.,  2013, @doi []
  10.1088/0067-0049/208/1/9, http://adsabs.harvard.edu/abs/2013ApJS..208....9P 208, 9

[Phillips, Shahbaz  &
  PodsiadlowskiPhillips et al.1999]Phillips99
Phillips S. N.,  Shahbaz T.,   Podsiadlowski P.,  1999, @doi []
  10.1046/j.1365-8711.1999.02357.x, http://adsabs.harvard.edu/abs/1999MNRAS.304..839P 304, 839

[RobertsRoberts2013]Roberts13
Roberts M. S. E.,  2013, in van Leeuwen J.,  ed.,  IAU Symposium Vol. 291,
  Neutron Stars and Pulsars: Challenges and Opportunities after 80 years. pp
  127–132 (@eprint arXiv 1210.6903), @doi10.1017/S174392131202337X

[Romani & SanchezRomani &
  Sanchez2016]Romani16
Romani R. W.,  Sanchez N.,  2016, @doi []
  10.3847/0004-637X/828/1/7, http://adsabs.harvard.edu/abs/2016ApJ...828....7R 828, 7

[Romani, Filippenko  & CenkoRomani
  et al.2015]Romani15
Romani R. W.,  Filippenko A. V.,   Cenko S. B.,  2015, @doi []
  10.1088/0004-637X/804/2/115, http://adsabs.harvard.edu/abs/2015ApJ...804..115R 804, 115

[Schroeder & HalpernSchroeder &
  Halpern2014]Schroeder14
Schroeder J.,  Halpern J.,  2014, @doi []
  10.1088/0004-637X/793/2/78, http://adsabs.harvard.edu/abs/2014ApJ...793...78S 793, 78

[ShahbazShahbaz1998]Shahbaz98
Shahbaz T.,  1998, @doi [] 10.1046/j.1365-8711.1998.01618.x, http://adsabs.harvard.edu/abs/1998MNRAS.298..153S 298, 153

[ShahbazShahbaz2003]Shahbaz03b
Shahbaz T.,  2003, @doi [] 10.1046/j.1365-8711.2003.06258.x, http://adsabs.harvard.edu/abs/2003MNRAS.339.1031S 339, 1031

[Shahbaz, Groot, Phillips, Casares,
  Charles  & van ParadijsShahbaz et al.2000]Shahbaz00
Shahbaz T.,  Groot P.,  Phillips S. N.,  Casares J.,  Charles P. A.,
   van Paradijs J.,  2000, @doi []
  10.1046/j.1365-8711.2000.03341.x, http://adsabs.harvard.edu/abs/2000MNRAS.314..747S 314, 747

[Shahbaz, Zurita, Casares, Dubus,
  Charles, Wagner  & RyanShahbaz et al.2003]Shahbaz03a
Shahbaz T.,  Zurita C.,  Casares J.,  Dubus G.,  Charles P. A.,
  Wagner R. M.,   Ryan E.,  2003, @doi [] 10.1086/346001, http://adsabs.harvard.edu/abs/2003ApJ...585..443S 585, 443

[Shahbaz, Casares, Watson, Charles,
  Hynes, Shih  & SteeghsShahbaz et al.2004]Shahbaz04
Shahbaz T.,  Casares J.,  Watson C. A.,  Charles P. A.,  Hynes R. I.,
   Shih S. C.,   Steeghs D.,  2004, @doi [] 10.1086/426504, http://adsabs.harvard.edu/abs/2004ApJ...616L.123S 616, L123

[Shahbaz, Watson  & DhillonShahbaz
  et al.2014]Shahbaz14
Shahbaz T.,  Watson C. A.,   Dhillon V. S.,  2014, @doi []
  10.1093/mnras/stu267, http://adsabs.harvard.edu/abs/2014MNRAS.440..504S 440, 504

[Stappers, van Kerkwijk, Bell  &
  KulkarniStappers et al.2001]Stappers01
Stappers B. W.,  van Kerkwijk M. H.,  Bell J. F.,   Kulkarni S. R.,
  2001, @doi [] 10.1086/319106, http://adsabs.harvard.edu/abs/2001ApJ...548L.183S 548, L183

[Storn & PriceStorn &
  Price1995]Storn95
Storn R.,  Price K.,  1995, Differential Evolution - A simple and efficient
  adaptive scheme for global optimization over continuous spaces

[Strader, Chomiuk, Sonbas, Sokolovsky,
  Sand, Moskvitin  & CheungStrader et al.2014]Strader14
Strader J.,  Chomiuk L.,  Sonbas E.,  Sokolovsky K.,  Sand D. J.,
  Moskvitin A. S.,   Cheung C. C.,  2014, @doi []
  10.1088/2041-8205/788/2/L27, http://adsabs.harvard.edu/abs/2014ApJ...788L..27S 788, L27

[Tang et al.,Tang et al.2014]Tang14
Tang S.,  et al., 2014, @doi [] 10.1088/2041-8205/791/1/L5, http://adsabs.harvard.edu/abs/2014ApJ...791L...5T 791, L5

[VrugtVrugt2016]Vrugt2016
Vrugt J. A.,  2016, @doi [Environmental Modelling     Software]
  http://dx.doi.org/10.1016/j.envsoft.2015.08.013, 75, 273

[Wang, Archibald, Thorstensen, Kaspi,
  Lorimer, Stairs  & RansomWang et al.2009]Wang09
Wang Z.,  Archibald A. M.,  Thorstensen J. R.,  Kaspi V. M.,  Lorimer
  D. R.,  Stairs I.,   Ransom S. M.,  2009, @doi []
  10.1088/0004-637X/703/2/2017, http://adsabs.harvard.edu/abs/2009ApJ...703.2017W 703, 2017

[Welsh, Horne  & GomerWelsh
  et al.1995]Welsh95
Welsh W. F.,  Horne K.,   Gomer R.,  1995, @doi []
  10.1093/mnras/275.3.649, http://adsabs.harvard.edu/abs/1995MNRAS.275..649W 275, 649

[van Kerkwijk, Breton  & Kulkarnivan
  Kerkwijk et al.2011]Kerkwijk11
van Kerkwijk M. H.,  Breton R. P.,   Kulkarni S. R.,  2011, @doi
  [] 10.1088/0004-637X/728/2/95, http://adsabs.harvard.edu/abs/2011ApJ...728...95V 728, 95

[van Staden & Antoniadisvan Staden &
  Antoniadis2016]Staden16
van Staden A. D.,  Antoniadis J.,  2016, @doi []
  10.3847/2041-8213/833/1/L12, http://adsabs.harvard.edu/abs/2016ApJ...833L..12V 833, L12






