# The Progenitor of the Vela Pulsar

C. S. Kochanek1,2
1 Department of Astronomy, The Ohio State University, 140 West 18th Avenue,
Columbus OH 43210
2 Center for Cosmology and AstroParticle Physics, The Ohio State University,
191 W. Woodruff Avenue, Columbus OH 43210

###### Abstract

With Gaia parallaxes it is possible to study the stellar populations
associated with individual Galactic supernova remnants (SNR) to estimate the
mass of the exploding star. Here we analyze the luminous stars near the Vela
pulsar and SNR to find that its progenitor was probably ($\mathrel{\raise
1.29167pt\hbox{$>$}\mkern-14.0mu\lower 2.58334pt\hbox{$\sim$}}90\%$) low mass
($8.1$-$10.3M_{\odot}$). The presence of the O star $\gamma^{2}$ Vel a little
over 100 pc from Vela is the primary ambiguity, as including it in the
analysis volume significantly increases the probability (to 5%) of higher mass
($>20M_{\odot}$) progenitors. However, to be a high mass star associated with
$\gamma^{2}$ Vel’s star cluster at birth, the progenitor would have to a
runaway star from an unbound binary with an unusually high velocity. The
primary impediment to analyzing large numbers of Galactic SNRs in this manner
is the lack of accurate distances. This can likely be solved by searching for
absorption lines from the SNR in stars as a function of distance, a method
which yielded a distance to Vela in agreement with the direct pulsar parallax.
If Vela was a $10M_{\odot}$ supernova in an external galaxy, the $50$ pc
search region used in extragalactic studies would contain only $\simeq 10\%$
of the stars formed in a 50 pc region around the progenitor at birth and
$\simeq 90\%$ of the stars in the search region would have been born
elsewhere.

###### keywords:

stars: massive – supernovae: general – supernovae

##  1 Introduction

We would like to understand which massive stars explode as supernovae and the
nature of the resulting compact objects. Modern theoretical models (e.g.,
O’Connor & Ott 2011, Pejcha & Thompson 2015, Ertl et al. 2016, Sukhbold et al.
2016, Ghosh et al. 2021) find a complex mapping between progenitor mass and
explosion driven by changes in core structure related to the balance between
radiative and convective carbon burning (Sukhbold et al. 2018, Sukhbold &
Adams 2020). In these modern models, essentially all stars either explode and
produce neutron stars or fail to explode and produce black holes – “fall back”
supernovae where the star explodes but a significant amount of mass falls back
onto the proto-neutron star to form a black hole are extremely rare.111Despite
this, fall back models for remnant masses such as Fryer et al. (2012) remain
in common use for binary population synthesis models (e.g., Belczynski et al.
2020, Breivik et al. 2020, Chawla et al. 2021, Eldridge et al. 2017, Riley et
al. 2021). In the absence of fall back, lower mass ($\mathrel{\raise
1.29167pt\hbox{$<$}\mkern-14.0mu\lower 2.58334pt\hbox{$\sim$}}10M_{\odot}$)
black holes are produced either from explosions of stars stripped by pre-
supernova mass loss/transfer or by the failed explosions of the more massive
red supergiants. These produce black holes with the mass of the helium core
(Kochanek 2014) because the Nadezhin (1980) mechanism (also see, e.g.,
Lovegrove & Woosley 2013, Fernández et al. 2018) ejects the weakly bound
hydrogen envelope of the supergiant. None of these surveys of outcomes are
based on full ab initio core collapse simulations but instead use “calibrated”
explosion models to explore outcomes. The outcomes in true core collapse
simulations remain an open problem because of the complexity of the physics
and the need for high resolution three-dimensional simulations (e.g., Bollig
et al. 2021, Burrows et al. 2020, Pan et al. 2021).

To test these theoretical predictions, we need to observationally determine
the mapping between progenitors and outcomes. The cleanest approach for the
stars which explode is simply to measure the properties of the progenitor
star, as first done for SN 1987A (Gilmozzi et al. 1987). Since the stellar
luminosity is determined by the mass of the core, and the mass of the core is
determined by the initial mass of the progenitor this is a fairly robust
approach if there is sufficient data to well-determine the luminosity
(although the mass at the time of explosion cannot be well-constrained, see
Farrell et al. 2020). The challenge is that this must be done in distant
($\sim 1$-$10$ Mpc) galaxies and largely depends on the existence of multi-
band archival Hubble Space Telescope data for robust results. The progenitors
of Type IIP supernovae are red supergiants with an upper mass limit that is
consistent with the theoretical explosion studies (Smartt et al. 2009, Smartt
2009, Smartt 2015, although there are rebutted (e.g., Kochanek et al. 2012,
Kochanek 2020, Beasor et al. 2020) counterarguments (e.g., Walmswell &
Eldridge 2012, Groh et al. 2013, Davies & Beasor 2020). Less is known about
the progenitors of Type Ibc supernovae because they are generally undetected,
which likely implies that in most cases the envelopes are stripped by binary
processes (e.g., Eldridge et al. 2013, Folatelli et al. 2016, Johnson et al.
2017, Kilpatrick et al. 2021). Unfortunately, the nature of the compact
remnant formed in these systems will likely always be unknown.

Failed supernovae can be found by searching for stars which disappear
independent of the nature of any associated transient (Kochanek et al. 2008).
A search for failed supernovae using the Large Binocular Telescope (Gerke et
al. 2015, Adams et al. 2017b, Neustadt et al. 2021) has identified one
excellent candidate (Gerke et al. 2015, Adams et al. 2017a, Basinger et al.
2020) and a second, weaker candidate (Neustadt et al. 2021). As expected, the
progenitor of the strong candidate was a massive RSG, and finding one
candidate implies a fraction of core collapses leading to failed supernovae
consistent with the current theoretical predictions.

The primary alternative to searching for individual progenitors is to use the
local stellar population to infer the probable mass of the progenitor. This
has been done for supernova remnants in the Magellanic Clouds (Badenes et al.
2009, Auchettl et al. 2019) and in nearby galaxies (e.g., Jennings et al.
2012, Jennings et al. 2014, Díaz-Rodríguez et al. 2018, Williams et al. 2018,
Williams et al. 2019, Díaz-Rodríguez et al. 2021, Koplitz et al. 2021). Except
for Williams et al. (2018) and Díaz-Rodríguez et al. (2021), which examine
stellar populations near known supernovae, these are studies of the stellar
populations near supernova remnants. The primary advantage of this method is
that it can be applied to large numbers of supernova remnants or historical
supernovae compared to the numbers of directly observed progenitors. For the
Magellanic Clouds, several of the remnants are associated with neutron stars
(N49, N158A and possibly SN 1987A in the LMC, and 1KT6 and 1E 0102.2$-$7219 in
the SMC, see Badenes et al. 2009, Auchettl et al. 2019), so the outcome of the
explosion is also known.

There are also several disadvantages. First, unless the local stellar
population has a single well-defined starburst, the method does not provide
individual well-constrained masses and so is best suited for making
statistical models of the progenitor distribution. Second, the nature of the
supernova is unknown for SNRs, so, for example, one cannot separately
investigate the progenitors of Type IIP and Type Ibc supernovae. No SNR
studies can directly address the deficit of more massive RSG progenitors to
Type IIP supernovae because the supernova types are unknown. Third, the life
time and detectability of the SNR has some dependence on the nature of the
explosion (e.g., Sarbadhicary et al. 2017, Jacovich et al. 2021) which will
introduce some bias into the results. Most, but not all, of these studies have
supported a deficit of higher mass SN progenitors. The progenitor mass
distribution models used to date in these studies are simple functional forms
that do not resemble current theoretical expectations, so it is difficult to
evaluate the degree to which they agree or disagree with these expectations.

With the advent of Gaia (Gaia Collaboration et al. 2016, Gaia Collaboration et
al. 2021) it is now possible to apply this second method to Galactic supernova
remnants if the distance to the remnant is known and the extinction is not
severe. The present number of such systems is small, primarily because so few
supernova remnants have well-constrained distances. As with extragalactic
SNRs, the supernova type will generally be unknown, but, unlike extragalactic
SNRs, the compact object outcome of the explosion frequently is known. Because
the Galaxy is so well-surveyed from the mid-IR into the ultraviolet, it will
generally be possible to characterize the individual stars extremely well,
with well-determined individual stellar luminosities, temperatures and
extinctions. In many cases, the brighter stars will also have spectroscopic
classifications. Here we demonstrate this for the Vela pulsar (PSR
J0835$-$4510, PSR B0833$-$45, Large et al. 1968), which is so nearby (280 pc,
Dodson et al. 2003) that many of the massive stars near it are visible to the
naked eye. We describe the selection of the stars and spectral energy
distribution (SED) models of the more luminous stars in §2. We analyze these
stars to estimate the likely mass of Vela’s progenitor in §3, and we discuss
the results and future prospects in §4.

##  2 Selecting the Stars

We select the stars from the Gaia EDR3 catalog (Gaia Collaboration et al.
2016, Gaia Collaboration et al. 2021), requiring them to have parallaxes,
proper motions and all three Gaia magnitudes. We apply no restriction on the
RUWE statistic for the quality of the parallax. We adopt the position (J2000
08:35:20.61149$\pm$0.00002, $-$45:10:34.8751$\pm$0.0003), parallax
($\varpi=3.5\pm 0.2$ mas) and proper motions ($\mu_{\alpha}=-49.68\pm 0.06$
and $\mu_{\delta}=29.9\pm 0.1$ mas/year) of the pulsar from Dodson et al.
(2003). The parallax is sufficiently accurate that it is essentially
unaffected by Lutz-Kelker bias (see Verbiest et al. 2012). The spin down age
of Vela is $11.3$ kyr, although Lyne et al. (1996) argue that the braking
index implies a larger age. The proper motions combined with the offset from
the center of the remnant imply an age of $18000\pm 9000$ years based on the
$25\pm 5$ arcmin offset of the pulsar from the center of the remnant
(Aschenbach et al. 1995). Unfortunately, Aschenbach et al. (1995) do not
report their estimated position for the center, but if we adopt an age of
$T=20000$ years we can estimate the position of the explosion from the
position shift of ($+17$,$-10$) arcmin. The remnant itself has a diameter of
approximately $8.3$ degrees, so the shift is small compared to the scale of
the remnant and the size of our search region. We select the stars in a region
centered on the shifted equatorial position of ($128.84^{\circ}$,
$-45.18^{\circ}$) where we have reduced the coordinate precision due to the
uncertain age.

The next question is the size of the region to extract around Vela to obtain a
representative stellar population. The studies of the regions around
supernovae and SNRs in nearby galaxies have almost uniformly used a region 50
pc in radius based on the argument that stars form in compact clusters ($\sim
1$ pc, Lada & Lada 2003), drift apart at low relative velocities ($\sim$
km/s), and that stars which explode after their parent binary was disrupted by
a supernova have only modest velocities (Eldridge et al. 2011, Renzo et al.
2019). With Gaia, we can simply consider this question empirically. In this
section we consider the three dimensional geometry we use for Vela, and in §4
we consider the two-dimensional projected geometry of extragalactic analyses.

We started by extracting all stars within $R=250$ pc of Vela. If the pulsar is
at a distance of $D=286$ pc, a sphere of radius $R$ subtends a maximum angle
of $\theta=\sin^{-1}R/D\simeq 62$ deg relative to its center. We select stars
with $\theta<62$ deg, $1.87<\varpi<28.6$ mas, which corresponds to $\pm 250$
pc around the pulsar distance, and $G<9.5$ mag. This will include all stars
with $M_{G}\mathrel{\raise 1.29167pt\hbox{$<$}\mkern-14.0mu\lower
2.58334pt\hbox{$\sim$}}0$ mag since the extinctions are small. At this
magnitude limit we are including all $M\mathrel{\raise
1.29167pt\hbox{$>$}\mkern-14.0mu\lower 2.58334pt\hbox{$\sim$}}5M_{\odot}$
stars and the more luminous, evolved, lower mass stars. The initial search
yields $\sim 55,000$ stars, of which $\sim 37000$ are within $R=250$ pc of the
pulsar, and $3160$ are brighter than $M_{G}<0$ mag.

Figure 1:  Color-magnitude diagram without extinction corrections for the
stars with $M_{G}>0$ mag and within $R<125$ pc of Vela. The background
contours are the maximum likelihood model for the stellar density including
extinction with solid (dashed) contours at higher (lower) densities that are
spaced by factors of $3$. The red curves are Solar metallicity PARSEC
isochrones with ages of $10^{8.5}$ (bottom, dashed), $10^{8.0}$ (dashed)
$10^{7.5}$ (solid), $10^{7.0}$ (solid) and $10^{6.6}$ (top, solid) years. The
maximum mass for the $10^{7.5}$ year isochrone is $9.1M_{\odot}$, so only the
solid isochrones still have stars which will explode as supernovae. The
isochrones are shifted by the mean extinction of $E(B-V)=0.057$ mag.

The bright magnitude limit of Gaia at $G\simeq 3$ mag can be a problem because
the most luminous stars may not be present in the Gaia catalog. We searched
the Hipparcos (Perryman et al. 1997, using the updated astrometric solution
from van Leeuwen 2007) and the Bright Star (Hoffleit & Warren 1995) catalogs
for any such stars in the search volume, finding nine Hipparcos stars: HIP
30324 ($\beta$ CMa) HIP 33579 ($\epsilon$ CMa), HIP 39953 ($\gamma^{2}$ Vel),
HIP 41037 ($\epsilon$ Car) HIP 44816 ($\lambda$ Vel), HIP 66657 ($\epsilon$
Cen) HIP 68702 ($\zeta$ Cen), and HIP 81173 ($\alpha$ TrA). Since $\gamma^{2}$
Vel is in a wide binary (Tokovinin 2018) with $\gamma^{1}$ Vel (HD 68243), we
use the Gaia parallax of $\gamma^{1}$ Vel for $\gamma^{2}$ Vel. For the other
stars we adopt the updated Hipparcos parallaxes. Including the bright
Hipparcos stars, there are $N=3169$ stars with $M_{G}<0$ mag within $R=250$ pc
of Vela. Figure 1 shows the distribution of the $R<125$ pc stars in absolute
magnitude $M_{G}$ and color $B_{P}-R_{P}$ (with no extinction corrections) as
compared to Solar metallicity PARSEC (Bressan et al. 2012, Marigo et al. 2013,
Pastorelli et al. 2020) isochrones with ages of $10^{6.6}$, $10^{7.0}$,
$10^{7.5}$, $10^{8.0}$ and $10^{8.5}$ years. For the bright Hipparcos stars we
used the SED fits described below to synthesize the Gaia magnitudes.

To explore optimizing the region from which we select stars, we only want to
consider the younger, massive stars. The main contaminants are lower mass red
giants, which we can largely remove by considering only stars with $M_{G}<0$
for $B_{P}-R_{P}<0.6$, $M_{G}<-2.0$ for $0.6<B_{P}-R_{P}<1.8$, and
$M_{G}<-3.5$ for $B_{P}-R_{P}>1.8$. Only a small fraction of the stars have
Gaia DR2 radial velocities, so we assigned a random line of sight velocity
with a dispersion of $3.0$ km/s based on the velocity dispersion derived from
the proper motions of the stars within 50 pc of Vela (see below).

In our present analysis we are selecting the stars in three dimensions, so we
want to explore the effect of changing the radius of the selection sphere. We
can estimate the completeness and contamination of the selection by examining
the relative positions of the stars today and in the past. We start by taking
the stars near the pulsar and determining where they were $10^{7}$ years and
$10^{7.5}$ years in the past – these roughly correspond to the lifetimes of
$20M_{\odot}$ and $10M_{\odot}$ stars, respectively. Because the pulsar
received a kick, we cannot do the same for the past location of the progenitor
star. For simplicity, we assume that the progenitor had the mean velocity of
the stars. The completeness will be lower if we include a random motion for
the progenitor.

Following the extragalactic studies, we first considered a 50 pc sphere. This
sphere around Vela today contains 19 of these young, high mass stars. After
subtracting their median tangential velocities, they have a two-dimensional
(2D) velocity dispersion of 4.2 km/s defined by half the width of the velocity
range encompassing 68% of the stars. As we expand the radius of the sphere,
the number of stars rises rapidly to 152, 411, 757 and 1160 within 100, 150,
200 and 250 pc, and the 2D velocity dispersion rises slowly to $5.4$, $6.7$,
$6.7$ and $7.9$ km/s, respectively. We used the 4.2 km/s dispersion derived
from the proper motions to set the $4.2/\sqrt{2}=3.0$ km/s dispersion of the
randomly assigned line of sight velocity. The three dimensional velocity
dispersion is then $5.1$ km/s, so the typical typical massive has a random
motion of $\simeq 50$ pc in the lifetime of a $\sim 20M_{\odot}$ star ($\sim
10^{7}$ years) and $\simeq 150$ pc in the lifetime of a $\sim 10M_{\odot}$
star ($\sim 10^{7.5}$ years).

If we place a sphere at the median position of these stars $10^{7}$ years ago,
we can estimate the completeness and contamination. The completeness is the
fraction of stars inside a sphere of some radius today that are inside a
sphere of some other radius centered at the median position of the stars in
the past. The contamination is the fraction of the stars in the “past sphere”
that were not in the “present sphere”. For a 50 pc sphere today, 4, 12, 16 and
18 of the 19 stars are in a past sphere of radius $50$, $100$, $150$, and
$200$ pc, along with 10, 61, 179 and 356 stars which were not in the 50 pc
sphere today. So the completenesses are 21%, 63%, 84% and 95%, while the
contamination rates are 71%, 84%, 92% and 95%. Basically, since the typical
random motion corresponds to moving $50$ pc in $10^{7}$ years, a sphere of
comparable size must suffer from poor completeness and significant
contamination. If we increase the size of the present sphere to $100$ pc, the
completenesses in past spheres of radius $100$, $150$ and $200$ pc are 36%,
68% and 84%, respectively, and the contamination rates are 21%, 47% and 67%.
Finally, if we use $150$ pc, the completeness/contamination are 40%/14% and
70%/25% for past spheres of radius $150$ and $200$ pc. This obviously gets
much worse if we use the $10^{7.5}$ year lifetime of a $\simeq 10M_{\odot}$
star. For a $50$ pc sphere today, and 50, 100, 150 and 200 pc spheres
$10^{7.5}$ years ago, the completeness/contamination rates are 5%/75%,
16%/79%, 26%/88% and 47%/89%, respectively (3%/61%, 9%/61%, and 22%/58% for
$100$ pc today; 7%/15% and 15%/25% for $150$ pc today).

Completeness and contamination are, or course, a trade off – there is no right
answer. We will present all of our results for a radius of $125$ pc. We were
originally going to use $100$ pc but the most massive star near Vela,
$\gamma^{2}$ Vel, lies just outside a $100$ pc sphere, so we increased the
radius to $125$ pc to double the volume and include $\gamma^{2}$ Vel. For the
final results on the probable age and mass of Vela’s progenitor, we will
present the results for a range of radii.

We fit the spectral energy distributions (SEDs) of the more luminous stars in
the $125$ pc sphere to estimate their luminosity, temperature and extinctions.
We also fit the bright Hipparcos stars to synthesize their Gaia magnitudes.
The SED fits are moderately labor intensive in terms of collating the data, so
we did not model all the lower luminosity stars. We initially selected all
stars with $M_{G}<-2$ ($M_{G}<-3.5$) for $B_{P}-R_{P}<1.8$ ($>1.8$), where the
higher luminosity limit for redder stars eliminates lower mass red giants.
While the SEDs showed no need to include circumstellar dust, we used the same
methods as in Adams et al. (2017a), running DUSTY (Elitzur & Ivezić 2001)
inside a Markov Chain Monte Carlo (MCMC) driver to both optimize the fits and
then estimate the uncertainties. We used Castelli & Kurucz (2003) model
atmospheres for all but the coolest stars where we used MARCS (Gustafsson et
al. 2008) model atmospheres. We used near-IR data from 2MASS (Cutri et al.
2003). The optical data primarily came from Johnson et al. (1966) and Tonry et
al. (2018), supplemented by Cousins (1971), Ducati (2002), Mermillod (1977),
Morel & Magnenat (1978), NOMAD (Zacharias et al. 2005), and Neckel & Klare
(1980). Almost all of the hot stars had UV data extending to $\sim 1500$ Å
from either Thompson et al. (1978) or Wesselius et al. (1982).

We used temperature priors based on the reported spectral types and weak
extinction priors. The temperature prior widths were roughly one spectral type
(so B1 to B3 for a B2 star) and $\pm 0.1$ mag for $E(B-V)$. For extinction
priors, we use the three dimensional combined19 mwdust models (Bovy et al.
2016), which combine the Drimmel et al. (2003), Marshall et al. (2006) and
Green et al. (2019) models to provide estimates for any sky position. We
extracted the V band extinction, and then used an $R_{V}=3.1$ extinction law
to convert the V band extinction to those for the $G$, $B_{P}$ and $R_{P}$
bands. For most stars, the SEDs could only modestly improve the spectral
temperature estimates but strongly constrained the extinction. The agreement
with the mwdust estimates was generally good, except when the mwdust estimate
was high ($E(B-V)\mathrel{\raise 1.29167pt\hbox{$>$}\mkern-14.0mu\lower
2.58334pt\hbox{$\sim$}}0.1$). In these cases, the SED models generally
required much less dust. With this caveat, we will use the mwdust extinctions
to model the effects of extinction for the full sample.

Table 2 gives the goodness of fits, temperature, luminosity, mass, age,
separation from estimated explosion center and some comments for the stars
with $L_{*}>10^{3.5}L_{\odot}$. The individual masses and ages are simply the
range of PARSEC isochrone ages and masses where the luminosity and temperature
are within twice the estimated uncertainties of the fitted values with a
minimum uncertainty of $0.02$ dex. The ages and masses are strongly correlated
- the maximum masses correspond to the minimum ages and vice versa. Except for
$\gamma^{2}$ Vel, they are all less massive than $\simeq 15M_{\odot}$ and
older than $\simeq 10^{7}$ years even for the upper (lower) limits on the
masses (ages).

While many of the stars are binaries (see Table 2), the extra light from the
secondary seems to have little consequence for the SED fits. This is not very
surprising because B stars are still in the regime where the mass-luminosity
relation is fairly steep and a modestly lower mass companion will not greatly
perturb the SED. For example, consider the most luminous star, $\gamma^{2}$
Vel, which is an O star plus Wolf-Rayet (WR) star binary. $\gamma^{2}$ Vel is
an interferometrically resolved double lined spectroscopic binary with present
day masses of $(28.5\pm 1.1)M_{\odot}$ and $(9.0\pm 0.6)M_{\odot}$ for the O
star and the companion WR star, respectively (North et al. 2007). While our
SED fit finds a higher total luminosity, most of the difference is due to both
the Gaia EDR3 distance (where we used the distance to the wide binary
companion $\gamma^{1}$ Vel) and our extinction estimate being larger than
those used by North et al. (2007). Still, our rough individual mass estimate
of $25.3$-$27.3M_{\odot}$ agrees reasonably well with their dynamical
measurement. North et al. (2007) estimate that only $22\%$ of the V-band light
comes from the WR star, so the SED fit itself is not strongly biased by the WR
companion.

Table 1: Luminous Stars Near the Vela Pulsar

Star | $\chi^{2}/N_{dof}$ | $\log(T_{*})$ | $\log(L_{*})$ | $M_{*}$ | $\log t$ | Sep | Comments
---|---|---|---|---|---|---|---
|  | (K) | ($L_{\odot}$) | ($M_{\odot}$) | (yr) | (pc) |
HD 68273 | $0.79$ | $4.622\pm 0.020$ | $6.029\pm 0.070$ |  $25.3$-$27.3$ |  $6.61$-$6.63$ | 102 |  $\gamma^{2}$ Vel, O7.5III+WC8 binary
HD 63462 | $1.82$ | $4.383\pm 0.032$ | $4.347\pm 0.094$ |  $9.8$-$13.9$ |  $6.96$-$7.40$ | 124 | O Pup, B1IVe
HD 68243 | $1.00$ | $4.358\pm 0.033$ | $4.319\pm 0.084$ |  $9.5$-$13.2$ |  $7.07$-$7.43$ | 102 |  $\gamma^{1}$ Vel, B2III
HD 65818 | $1.74$ | $4.390\pm 0.029$ | $4.297\pm 0.075$ |  $9.9$-$13.4$ |  $6.92$-$7.39$ | 92 | V Pup, B1Vp+B2 ecl. binary
HD 74575 | $0.27$ | $4.345\pm 0.018$ | $4.149\pm 0.053$ |  $8.9$-$11.3$ |  $7.19$-$7.49$ | 70 |  $\alpha$ Pyx, B1.5III
HR 3307 | $2.29$ | $3.624\pm 0.012$ | $4.137\pm 0.029$ |  $7.9$-$13.1$ |  $7.19$-$7.62$ | 116 |  $\epsilon$ Car, K3III+B2 binary
HD 57150 | $1.03$ | $4.325\pm 0.040$ | $4.076\pm 0.121$ |  $7.7$-$11.5$ |  $7.10$-$7.61$ | 83 |  $\upsilon^{1}$ Pup, B2V+B3IV binary
HD 78647 | $2.03$ | $3.618\pm 0.010$ | $4.050\pm 0.044$ |  $2.0$-$12.1$ |  $7.25$-$8.73$ | 121 |  $\lambda$ Vel, K4I
HD 63032 | $1.53$ | $3.619\pm 0.010$ | $4.008\pm 0.036$ |  $2.0$-$11.7$ |  $7.27$-$8.74$ | 84 | c Pup, K4III
HD 65551 | $1.61$ | $4.407\pm 0.014$ | $3.988\pm 0.034$ |  $9.8$-$11.9$ | $<7.27$ | 51 | N Pup, B2II/IV
HD 83058 | $0.57$ | $4.345\pm 0.035$ | $3.940\pm 0.098$ |  $7.6$-$10.8$ |  $6.84$-$7.62$ | 90 | L Vel, B2IV
HD 56139 | $2.75$ | $4.246\pm 0.028$ | $3.931\pm 0.087$ |  $7.3$-$9.3$ |  $7.42$-$7.66$ | 120 |  $\omega$ CMa, B2.5Ve
HD 51799 | $2.58$ | $3.597\pm 0.008$ | $3.862\pm 0.034$ |  $1.9$-$10.0$ |  $7.38$-$8.75$ | 122 | M1III
HD 63465 | $0.93$ | $4.261\pm 0.032$ | $3.843\pm 0.108$ |  $6.6$-$9.3$ |  $7.41$-$7.75$ | 105 | B2IV/V
HD 68324 | $0.31$ | $4.353\pm 0.028$ | $3.835\pm 0.075$ |  $8.2$-$10.3$ |  $6.62$-$7.47$ | 65 | IS Vel, B2V
HD 64740 | $0.84$ | $4.369\pm 0.024$ | $3.800\pm 0.059$ |  $8.4$-$10.4$ | $<7.42$ | 56 | B2V
HD 73155 | $1.78$ | $3.671\pm 0.011$ | $3.799\pm 0.040$ |  $6.6$-$9.5$ |  $7.43$-$7.78$ | 90 | C Vel, K1.5I
HD 68217 | $0.34$ | $4.324\pm 0.034$ | $3.772\pm 0.094$ |  $6.9$-$9.7$ |  $6.79$-$7.70$ | 48 | B2IV
HD 64503 | $0.44$ | $4.307\pm 0.029$ | $3.727\pm 0.081$ |  $6.8$-$9.0$ |  $7.17$-$7.72$ | 64 | b Pup, B2V binary
HD 65575 | $0.73$ | $4.232\pm 0.029$ | $3.661\pm 0.099$ |  $6.0$-$8.1$ |  $7.53$-$7.84$ | 117 |  $\chi$ Car, B3IV
HD 57219 | $0.25$ | $4.263\pm 0.030$ | $3.643\pm 0.087$ |  $6.1$-$8.2$ |  $7.42$-$7.83$ | 102 |  $\upsilon^{2}$ Pup, B3Ve
HD 89388 | $1.63$ | $3.639\pm 0.008$ | $3.642\pm 0.035$ |  $5.2$-$8.7$ |  $7.50$-$8.03$ | 112 | q Car, K2.5II
HD 85123 | $0.75$ | $3.879\pm 0.018$ | $3.635\pm 0.041$ |  $6.0$-$7.4$ |  $7.64$-$7.86$ | 108 |  $\upsilon$ Car, A9
HD 54893 | $0.32$ | $4.282\pm 0.028$ | $3.602\pm 0.082$ |  $6.2$-$8.2$ |  $7.30$-$7.80$ | 83 | A Pup, B3IV/V
HD 63949 | $0.58$ | $4.286\pm 0.029$ | $3.563\pm 0.091$ |  $6.2$-$8.1$ |  $7.21$-$7.81$ | 112 | QS Pup, B2III $\beta$ Ceph
HD 79275 | $0.91$ | $4.320\pm 0.025$ | $3.509\pm 0.073$ |  $6.8$-$8.5$ | $<7.60$ | 63 | B2III/IV
HD 56779 | $0.61$ | $4.247\pm 0.024$ | $3.490\pm 0.074$ |  $5.7$-$7.3$ |  $7.53$-$7.89$ | 85 | B3V
HD 69144 | $0.79$ | $4.195\pm 0.025$ | $3.487\pm 0.084$ |  $5.4$-$7.1$ |  $7.66$-$7.94$ | 82 | NO Vel, B3III ecl. binary
HD 84816 | $1.06$ | $4.268\pm 0.026$ | $3.465\pm 0.082$ |  $5.9$-$7.5$ |  $7.36$-$7.86$ | 82 | B2V
HD 85622 | $2.08$ | $3.697\pm 0.012$ | $3.438\pm 0.046$ |  $5.3$-$7.2$ |  $7.67$-$8.01$ | 66 | m Vel, G6IIa binary
HD 63922 | $0.42$ | $4.152\pm 0.020$ | $3.429\pm 0.053$ |  $5.6$-$6.4$ |  $7.77$-$7.92$ | 40 | P Pup, B7/8
HD 61641 | $0.78$ | $4.246\pm 0.023$ | $3.402\pm 0.073$ |  $5.6$-$7.0$ |  $7.51$-$7.92$ | 106 | e Pup, B3III
HD 59550 | $0.89$ | $4.269\pm 0.020$ | $3.401\pm 0.060$ |  $6.2$-$7.2$ |  $7.36$-$7.74$ | 118 | B2IV
HD 71510 | $0.41$ | $4.298\pm 0.047$ | $3.394\pm 0.132$ |  $5.6$-$8.6$ | $<7.91$ | 62 | B3IV, Be star
HD 68895 | $0.70$ | $4.203\pm 0.026$ | $3.295\pm 0.086$ |  $5.0$-$6.5$ |  $7.67$-$8.04$ | 104 | B5V binary
HD 76566 | $1.02$ | $4.244\pm 0.025$ | $3.251\pm 0.080$ |  $5.5$-$6.6$ |  $7.36$-$7.84$ | 89 | B3IV
HD 75710 | $1.67$ | $3.959\pm 0.005$ | $3.251\pm 0.025$ |  $5.0$-$5.6$ |  $7.91$-$8.03$ | 101 | g Vel, A2III
HD 82150 | $4.19$ | $3.630\pm 0.008$ | $3.173\pm 0.042$ |  $2.5$-$6.2$ |  $7.81$-$8.93$ | 92 |  $\epsilon$ Ant, K3III
HD 91504 | $0.28$ | $3.675\pm 0.011$ | $3.156\pm 0.029$ |  $4.5$-$6.2$ |  $7.81$-$8.18$ | 121 | t Vel, K1/2III
HD 72227 | $0.77$ | $3.655\pm 0.012$ | $3.129\pm 0.038$ |  $3.4$-$6.1$ |  $7.83$-$8.53$ | 85 | K3III
HD 92449 | $0.49$ | $3.697\pm 0.013$ | $3.105\pm 0.049$ |  $4.5$-$5.9$ |  $7.86$-$8.18$ | 110 | x Vel, G5II
HD 50235 | $0.16$ | $3.646\pm 0.010$ | $3.087\pm 0.029$ |  $2.5$-$5.9$ |  $7.86$-$8.91$ | 111 | K2/3III
HD 85355 | $0.95$ | $4.100\pm 0.013$ | $3.057\pm 0.040$ |  $4.5$-$5.1$ |  $7.99$-$8.14$ | 66 | n Vel, B7III
HD 75630 | $1.42$ | $3.958\pm 0.006$ | $3.027\pm 0.029$ |  $4.4$-$4.8$ |  $8.06$-$8.17$ | 94 | h Vel, A2/EIV
HD 49689 | $0.48$ | $3.659\pm 0.015$ | $3.026\pm 0.051$ |  $3.2$-$5.7$ |  $7.88$-$8.59$ | 94 | K3III
HD 65699 | $0.22$ | $3.694\pm 0.013$ | $2.933\pm 0.043$ |  $4.2$-$5.3$ |  $7.96$-$8.27$ | 112 | 12 Pup, G9IIa pec
HD 68512 | $0.50$ | $3.668\pm 0.013$ | $2.922\pm 0.044$ |  $3.2$-$5.3$ |  $7.96$-$8.60$ | 80 | K3III
HD 85656 | $0.55$ | $3.678\pm 0.012$ | $2.863\pm 0.029$ |  $3.4$-$5.1$ |  $8.00$-$8.53$ | 101 | K1III
HD 76304 | $0.90$ | $3.692\pm 0.009$ | $2.708\pm 0.017$ |  $3.4$-$4.7$ |  $8.10$-$8.53$ | 119 | K0II/III
HD 80126 | $0.23$ | $3.726\pm 0.012$ | $2.696\pm 0.040$ |  $3.8$-$4.5$ |  $8.13$-$8.37$ | 118 | G8II

Figure 2:  Age distribution $N_{i}F_{i}$ of the modeled stars. The points show
the median number of stars associated with each age bin and the 90% confidence
range. The horizontal error bars span the bin widths, and the mass range
corresponding to the more massive age bins is listed. The solid red curves
shows the median integral distribution and the dashed curves show its 90%
confidence range. All uncertainties are highly correlated. The integral
distribution converges exactly to the number of stars because we are
determining how to distribute the modeled stars over the age bins.  Figure 3:
Inferred numbers $N_{i}$ of $M>M_{min}=M_{\odot}$ stars formed for each age
bin. This corresponds to the number of modeled stars shown in Figure 2 divided
by the fraction $F_{i}$ of $M>M_{min}$ stars which have $M_{G}>0$ mag. The
points show the median number of $M>M_{min}$ stars associated with each age
bin and the 90% confidence range. The horizontal error bars span the bin
widths, and the mass range corresponding to the more massive age bins is
listed. All uncertainties are highly correlated.  Figure 4:  Rate of forming
$M>M_{min}=M_{\odot}$ stars per million years for each age bin. This
corresponds to the number of $M>M_{min}$ stars $N_{i}$ shown in Figure 3
divided by the temporal bin width $\Delta t_{i}$. The points show the median
star formation rate associated with each age bin and the 90% confidence range.
The horizontal error bars span the bin widths, and the mass range
corresponding to the more massive age bins is listed. All uncertainties are
highly correlated.

Figure 5:  Mass distribution of the $M_{G}>0$ mag stars. The black solid shows
the maximum likelihood model and the dashed curves show the mean and 90%
confidence range for the integral distributions from the MCMC chains. This
distribution should be nearly complete for masses above $5M_{\odot}$ (the
vertical dashed line) and increasingly incomplete for lower masses. The red
curves with points show the integral distributions for the minimum and maximum
masses of the luminous stars with SED models from Table 2. The maximum mass
star is $\gamma^{2}$ Vel. The maximum (minimum) masses correspond to younger
(older) ages, and the CMD models strongly favor an older population and so
match the minimum mass curve well.  Figure 6:  The number of stars dying in
the last $10^{5}$ years for each age bin. The points show the probability for
each bin and its 90% confidence range. The horizontal error bars span the bin
widths, and the mass range corresponding to each age bin is listed. The solid
red curves shows the median integral distribution and the dashed curves show
its 90% confidence range. The arbitrarily normalized dotted cyan curve shows
what the age distribution would be for a constant star formation rate. All
uncertainties are highly correlated.  Figure 7:  Integral probability
distributions for the age of the progenitor of Vela as a function of the
radius $R$ used to select the stars. The dashed curve is for the $R=125$ pc
radius used until now. The large difference between this distribution and the
distribution for $R=100$ pc is due to the inclusion of $\gamma^{2}$ Vel in the
larger sphere.

##  3 The Progenitor of Vela

The striking property of Figures 1 and the models in Table 2 is the marked
absence of very young stars or massive stars with the exception of
$\gamma^{2}$ Vel. All of the other stars are only consistent with populations
older than ten million years, corresponding to maximum progenitor masses
$\mathrel{\raise 1.29167pt\hbox{$<$}\mkern-14.0mu\lower
2.58334pt\hbox{$\sim$}}15M_{\odot}$. The next step is to use these stars to
estimate the mass of the progenitor of the Vela pulsar. Here we will model the
Gaia color-magnitude distribution to estimate the numbers of stars as a
function of age. As discussed in §2, we present all of the results for a $125$
pc sphere around Vela, and the estimated age/mass of the progenitor for
regions from $100$ to $200$ pc in radius.

We will assume a Salpeter (1955) initial mass function (IMF) with a minimum
mass of $M_{min}=1M_{\odot}$ since we are not interested in the lower mass
stars. If the star formation rate ($M_{\odot}$/year) for $M>M_{min}$ stars is
$SFR$, then the rate of forming stars is

${dN\over dtdM}={(x-2)SFR\over M_{min}^{2}}\left({M\over M_{min}}\right)^{-x}$
(1)

with $x=2.35$ and a mean mass of $\langle M\rangle=(x-1)M_{min}/(x-2)$. We
divide the star formation history into logarithmic time intervals
$t_{min,i}<t<t_{max,i}$ with $\Delta t_{i}=t_{max,i}-t_{min,i}$. Assuming the
star formation rate $SFR_{i}$ in the interval is constant, the number of
$M>M_{min}$ stars formed is

$N_{i}={SFR_{i}\Delta t_{i}\over\langle M\rangle}.$ (2)

The number of stars formed in this period which die in a short time interval
$\delta t$ today is

$N_{i}{\delta t\over\Delta t_{i}}\left[\left({M(t_{min,i})\over
M_{min}}\right)^{1-x}-\left({M(t_{max,i})\over
M_{min}}\right)^{1-x}\right]=N_{i}S_{i}\delta t$ (3)

where $M(t)$ is the most massive surviving star on the isochrone, and
$S_{i}\delta t$ is the fraction of $M>M_{min}$ stars that died in the last
$\delta t$ years. This expression is explicitly assuming single star
evolution. We used $8$ temporal bins spanning $6.4<\log t<8.0$, bin widths of
$0.2$ dex and use the number $N_{i}$ of $M>M_{min}$ stars formed as the
variable to be determined.

If the global IMF is Salpeter (1955) down to $0.5M_{\odot}$ and then flattens
to $M^{-1.3}$ from $0.08M_{\odot}$ to $0.5M_{\odot}$ (Kroupa 2001), then the
$M>M_{\odot}$ stars represent 9.1% of the stars formed. The global mean mass
is $0.61M_{\odot}$, so for each $M>M_{\odot}$ star formed, the total mass of
new stars is $6.7M_{\odot}$. Local estimates of the mean stellar mass density
are $0.04M_{\odot}$ pc-3 (e.g., Flynn et al. 2006), so we should find $\sum
N_{i}\simeq 6\times 10^{4}$. If we use $2\pi R_{d}^{2}H$ with $R_{d}=3$ kpc
and $H=0.1$ kpc as the volume of the Galaxy, forming ten $M>M_{\odot}$ stars
per megayear in a 125 pc sphere corresponds to a global star formation rate of
$\sim 0.04M_{\odot}$/year.

We start by building density maps $\rho_{i}^{jk}$ in absolute magnitude
$M_{G}$ (index $j$) and color $B_{P}-R_{P}$ (index $k$) for each time interval
$i$. We use Solar metallicity PARSEC isochrones sampled at $\Delta\log t=0.01$
dex. For each density map $i$ we carry out $N_{trial}$ trials. We uniformly
select a time between $t_{min,i}$ and $t_{max,i}$, which corresponds to
assuming a constant star formation rate for each bin, randomly draw an initial
stellar mass from the IMF, and use the isochrones to determine the absolute
magnitude and color if a star of this mass still exists.

For the magnitude range we consider, the uncertainties in the Gaia magnitudes
are negligible. Similarly, the median parallax uncertainty of $0.027$ mas is
also unimportant given that the smallest parallax we consider is $2.43$ mas.
The median mwdust extinctions increases as $E(B-V)=0.19(d/\hbox{kpc})$ with
distance $d$ from the Sun. The width of the distribution is approximately
$0.13(d/\hbox{kpc})$. Since the mwdust estimates agreed reasonably well with
the extinction estimates from the SED fits, we randomly assigned each trial
the mwdust extinction of one of the $M_{G}<0$ mag stars.

If the extincted trial star has $M_{G}<0.0$ mag, we add $N_{trial}^{-1}$ to
the appropriate cell of $\rho_{i}^{jk}$. The maps span a finite range of color
($-0.75<B_{P}-R_{P}<3.5$) and absolute magnitude ($0.0>M_{G}>-8.0$) so trial
stars that fall outside the color range or are brighter than $-8.0$ mag are
placed at the bin edge. This preserves the total probability for $M_{G}<0.0$
mag stars. As we build the maps, we also build a binned mass distribution of
the $M_{G}<0.0$ mag stars, $D_{i}^{j}$ for mass bin $j$. For $M_{G}<0.0$ mag
we are including almost all $M>5M_{\odot}$ stars, and then only mass-dependent
portions of the post-main sequence lifetimes of the lower mass stars.

With these definitions, the number distribution of stars in magnitude and
color brighter than the selection limit from time period $i$ is
$N_{i}\rho^{i}_{jk}$ and $N_{i}\sum_{jk}\rho_{i}^{jk}=N_{i}F_{i}$ is the
number of stars still living and passing the selection criterion. $F_{i}$ is
the fraction of $M>M_{min}$ stars which are more luminous than $M_{g}>0$ mag.
The expected number of stars in a color/absolute magnitude bin is

$e_{jk}=\sum_{i}N_{i}\rho_{i}^{jk}.$ (4)

If we now distribute the $N$ actual stars over the grid, the observed number
of stars in a cell is $n_{jk}$ with a Poisson probability of finding that
number of stars of $e_{jk}^{n_{jk}}\exp(-e_{jk})/n_{jk}!$. The logarithm of
the likelihood for all $N$ stars is

$\ln L=\sum_{stars}\ln\left({re_{ij}^{n_{ij}}\over
n_{ij}!}\right)-\sum_{all}re_{ij}$ (5)

where the first term is the sum over bins containing stars and the second is
the sum over all bins. Note that the factorial $n_{ij}!$ can be discarded
since the calculation depends only on likelihood differences and not the
absolute likelihood. Putting trial stars falling off the grid when
constructing the distributions on the bin edges ensures that the second sum is
correct.

In Equation 5 we have also introduced a “renormalization” factor $r$. Equation
5 with $r=1$ still includes the Poisson uncertainties from the total number N
of stars being modeled. For the problem of estimating the progenitor mass,
however, we need the relative probabilities (i.e. ratios) of the numbers of
stars in the age bins, not their absolute values, and the ratios are
unaffected by Poisson fluctuations from the finite number of stars. We solve
this problem by optimizing the likelihood with respect to the renormalization
factor to find that

$r=N\left[\sum_{all}e_{jk}\right]^{-1},$ (6)

which we then use to renormalize $N_{i}\rightarrow rN_{i}$. With this
renormalization, $\sum_{all}e_{jk}\equiv N$ and we have effectively converted
the Poisson likelihood into a multinomial likelihood for how to divide the $N$
stars over the age bins.

Operationally, we optimize the likelihood and estimate the uncertainties using
Markov Chain Monte Carlo (MCMC) methods with the $\log N_{i}$ as the
variables. Once trial values of $\log N_{i}$ are selected, they are
renormalized before computing the likelihood. Some age bins were susceptible
to $\log N_{i}\rightarrow-\infty$ (i.e., $N_{i}$ arbitrarily close to zero),
so we added a weak prior that the star formation rates of adjacent temporal
bins should be similar by adding

$\lambda^{-2}\sum_{i}\left[\ln\left({N_{i}\Delta t_{i+1}\over N_{i+1}\Delta
t_{i}}\right)\right]^{2}$ (7)

with $\lambda=6.91$ to the likelihood. This adds a penalty of unity to the
likelihood if adjacent bins have star formation rates ($SFR_{i}\propto
N_{i}/\Delta t_{i}$) that differ by a factor of 1000. This is just to prevent
numerical divergences and has no significant effect on the results.

Figure 1 shows contours for the numbers of stars as a function of luminosity
and temperature for the maximum likelihood model. Not surprisingly the maximum
density lies along the main sequence of B stars. The density contours follow
the stellar distribution closely. There are a few red giants in very low
density regions. This may mean that there some stars with higher extinctions
than our model, probably cases which required some circumstellar dust. They
are, however, old lower mass stars and should have no effect on the inferences
about young, high mass stars.

Figures 2, 3 and 4 show three different ways of viewing the distribution of
stars in age. Figure 2 shows the distribution of the modeled stars,
$N_{i}F_{i}$, in age. Because of the renormalization procedure, the numbers
add up exactly to the $N=607$ stars used. The five youngest age bins,
corresponding to progenitors more massive than $10.3M_{\odot}$ contain a
median of $N_{i}F_{i}=5$ modeled stars and fewer than 19 at 95% confidence.
For these young bins, we are basically counting all $M\mathrel{\raise
1.29167pt\hbox{$>$}\mkern-14.0mu\lower 2.58334pt\hbox{$\sim$}}5M_{\odot}$
stars, so the number of higher mass stars is lower. For ten $M>5M_{\odot}$
stars and a Salpeter (1955) IMF, there are $3.9$, $2.3$, $1.5$ and $0.9$ stars
more massive than $10$, $15$, $20$ and $30M_{\odot}$, respectively. The field
contains very few high mass stars, consistent with Figure 1. For comparison,
the $10^{7.4}$-$10^{7.6}$ year bin, which spans a final mass range of $8.1$ to
$10.3M_{\odot}$ that corresponds to the minimum mass range expected to produce
supernovae for single stars, has a median of $N_{i}F_{i}=95$ modeled stars.

Only a small fraction of the $M>M_{min}$ stars are luminous enough to be
selected with $M_{G}>0$ mag, so the actual number of stars associated with the
older bins is much larger. Figure 3 show the number $N_{i}$ of $M>M_{min}$
stars as a function of age. As expected, the total number of $M>M_{\odot}$
stars ($(6.8\pm 2.3)\times 10^{4}$) in the volume is much larger than the
$607$ we modeled. and the total number of $M>0.08M_{\odot}$ stars would be
another $\sim 11$ times larger. The implied stellar mass density of $(0.056\pm
0.019)M_{\odot}$ pc-3 is remarkably close to the local estimate of
$0.05M_{\odot}$ pc-3 (Flynn et al. 2006) given the huge correction required to
go from the number of $M_{G}>0$ mag stars to the total mass of stars - for the
oldest bin, the $\sim 50$ modeled $M_{G}>0$ mag stars represent some $\sim
20,000$ $M>M_{\odot}$ stars which then correspond to some $\sim 200,000$ stars
in total.

Finally, Figure 4 shows the star formation rates for each temporal bin as the
number of $M>M_{\odot}$ stars formed per $10^{6}$ years. As discussed above, a
formation rate of 10 $M>M_{\odot}$ stars per million years corresponds to a
global Galactic star formation rate of order $0.04M_{\odot}$/year, so the peak
star formation rate of $\sim 0.3M_{\odot}$/year on a global basis is not
impressive. The peak is, however, again found in the $10^{7.4}$-$10^{7.6}$
year bin. The star formation rates in the younger bins are at least an order
of magnitude lower, with even fewer new stars because the associated times
$\Delta t_{i}$ become shorter (i.e. the drop off in Figure 3 is sharper than
in Figure 4).

Figure 5 shows the model for the integral distribution of the selected stars
in mass

$D^{j}=\sum_{i}N_{i}D_{i}^{j}$ (8)

where the $D_{i}^{j}$ are the mass distributions associated with each temporal
bin. As expected from the distributions in Figures 1, the region contains few
higher mass stars. We also show the integral distribution of the mass
estimates from the SED models in Figure 5. Broadly speaking they are in good
agreement, although the results from the CMD models track the lower mass
limits much more closely than the upper mass limits. As noted earlier, at
roughly constant luminosity and temperature, there is a strong correlation
between age and mass for these main sequence B stars in the sense that the
minimum mass estimates correspond to the maximum ages and vice versa. The CMD
models strongly disfavor having significant numbers of stars with the young
ages associated with the maximum mass estimates and strongly favor ages
corresponding to the minimum mass estimates and so track the minimum mass
estimates. The CMD models constrain the age distribution better than the SED
models of individual stars.

The outlier in Figure 5 is $\gamma^{2}$ Vel. The model mass distribution based
on the CMD predicts only $\simeq 0.3$ ($\simeq 0.1$) stars with
$M>20M_{\odot}$ ($>30M_{\odot}$). The existence of $\gamma^{2}$ Vel is not a
huge statistical anomaly, as the likelihood of having one or more such stars
is $26\%$ ($10\%$). However, in the CMD (Figure 1), there are simply no other
stars between where the $10^{7.5}$ year isochrone (maximum mass
$9.1M_{\odot}$) starts to turn off the main sequence and $\gamma^{2}$ Vel. The
absence of these stars drives the model to make $\gamma^{2}$ Vel moderately
unlikely.

Figure 6 shows the number of stars $N_{i}S_{i}\delta t$ (Equation 3) predicted
to have died in the last $\delta t=10^{5}$ years for each age bin. That the
integral probability sums to near unity is happenstance – the distribution is
not normalized. The probability for the ages leading to supernovae is
dominated by the age range corresponding to initial masses of
$8.1$-$10.3M_{\odot}$. Note that for the range of ages producing supernovae,
the expected number of deaths in the last $10^{5}$ years is only $\sim 0.1$
and so the time scale for this volume to produce a supernovae is $\delta
t\simeq 10^{6}$ years, far longer than the life time of an SNR. This is not a
statistical problem – we chose this volume because it contained an SNR. To
estimate $\delta t$ we would have to analyze a much larger volume chosen
without using any prior knowledge of the number of enclosed remnants.

Figure 6 also shows what the prediction would be if the star formation rate
was constant across all age bins, roughly normalized to the intermediate age
bins where the estimated star formation rates are nearly constant (Figure 4).
The four youngest and two oldest bins have lower estimated SFRs, and so fewer
stars are predicted to have died in the last $10^{5}$ years, while the
$10^{7.4}$-$10^{7.6}$ year bin has a higher SFR and more predicted deaths. The
difference between the two curves represents the information added by using
the local stellar population to infer the star formation history.

The probability of a star dying in the last $\delta t=10^{5}$ years is
dominated by the older bins and the volume is likely to contain a number of
young white dwarfs. Out of curiosity we did a cursory search for candidates.
Based on the MIST (Dotter 2016, Choi et al. 2016) evolutionary tracks, a
$<10^{5}$ year old white dwarf should be a very hot star with a luminosity
$\sim 10^{2}L_{\odot}$. The Gaia EDR3 source ID#5319832121597913984 is the
bluest and most luminous source near the tip of the usual white dwarf cooling
sequence. It is flagged as a very high probability white dwarf candidate by
Gentile Fusillo et al. (2019) but with no estimates of its physical
properties. There are no UV or spectroscopic observations of it, but for an
assumed temperature of $40,000$ K it would have a luminosity of $\sim
2L_{\odot}$. Since the luminosity scales linearly with the temperature on the
Rayleigh-Jeans tail of the SED, this star has to be significantly older than
$10^{5}$ years even at twice the temperature. There are five similarly blue
but $\Delta M_{G}\sim-3$ mag more luminous stars (#5592257426113535104,
#5411002594979483392, #5442024044243183232, #5512856125193586176, and
#5539780553619327360) whose luminosities could be in the right range. Three of
them are spectroscopically classified as hot sub-dwarfs (#5592257426113535104,
Garrison & Hiltner (1973); #5442024044243183232, Barlow et al. 2013;
#5539780553619327360, Kilkenny et al. 1988, one (#5411002594979483392) is
spectroscopically classified as an A0 star (Nesterov et al. 1995, but this
difficult to reconcile with its location in the CMD), and #5512856125193586176
has no spectroscopic classifications. Possibly one of these latter two sources
is a very young white dwarf masquerading as a hot subdwarf.

To produce a final constraint on the mass of Vela’s progenitor we must
consider three remaining issues. First, we have to put a minimum mass in by
hand. When analyzing a single region around one target, we cannot determine a
minimum mass for explosion. This requires analyzing many such regions both
with and without SNRs. We must impose a minimum mass by simply dropping the
older age bins. The obvious choice is to keep only the $10^{7.4}$-$10^{7.6}$
year, with its maximum masses at death of $8.1$-$10.3M_{\odot}$, and younger
bins. For single star evolution, this is the correct mass range for the
cutoff. With binary evolution the next older bin can contribute through
mergers of longer lived, lower mass stars (see, e.g., Zapartas et al. 2017).
Retaining this next older bin would only strengthen our final conclusion that
the progenitor was significantly less massive than $20M_{\odot}$. The change
would be modest because many fewer stellar deaths are expected from this age
bin than from the $10^{7.4}$-$10^{7.6}$ year bin (see Figure 6).

In our formalism, the number of stellar deaths is $N_{i}S_{i}\delta t$ where
$S_{i}$ is independent of $\delta t$. For our final result we want the
relative probabilities of the age bins with no dependence on $\delta t$. The
probability that the progenitor came from bin $i$ and no other bin $j$ is

$\displaystyle P(i\&!j)$ $\displaystyle=$ $\displaystyle N_{i}S_{i}\delta
t\exp(-N_{i}S_{i}\delta t)\Pi_{i\neq j}\exp(-N_{j}S_{j}\delta t)$
$\displaystyle=$ $\displaystyle N_{i}S_{i}\delta
t\Pi_{all}\exp(-N_{J}s_{j}\delta_{j}).$

The total probability summed over all the bins is

$P_{tot}=\sum_{i}P(i\&!j)=\delta
t\left[\sum_{all}N_{i}S_{i}\right]\exp\left(-\delta
t\sum_{all}N_{i}S_{i}\right)$ (10)

and thus the normalized probability for each bin of

$P_{tot}^{-1}P(i\&!j)=N_{i}S_{i}\left[\sum_{all}N_{i}S_{i}\right]^{-1}$ (11)

is independent of $\delta t$ as desired. Basically, this just corresponds to
normalizing the integral probability distribution so that the total
probability is unity.

Finally, we have a maximum likelihood solution and all the MCMC samples, each
of which represents a realization of the integral probability distribution. We
could go to each age bin, sort these distributions and report a median and
some range, say 90% confidence, but it is unclear how to interpret this. What
does it mean to say there is a 5% chance that there is a 10% chance that the
progenitor was younger than the age of some bin? There really should be only
one distribution which incorporates this information. We are really combining
the $N_{MCMC}$ results of the MCMC chains, each of which has a probability of
$N_{MCMC}^{-1}$, and each trial predicts a probability for the age of the
progenitor. The way to combine them is to sum the probability the trial
predicts for the age bin weighted by the probability of the trial – in short,
the final probability distribution is simply the average of the MCMC samples.

Figure 7 shows these final probability distributions as a function of the
radius $R$ of the sphere used to select the stars. Recall that the results so
far have all been for $R<125$ pc. The maximum likelihood age distributions
modestly favor older, lower mass progenitors compared to these distributions.
For $R<100$ pc, a low mass progenitor is very strongly favored with a 94%
probability of it coming from the lowest mass bin ($8.1$-$10.3M_{\odot}$).
However, as we discussed in §2, this sphere probably only contains $\sim 2/3$
of the stars born within a similar radius. If we increase the radius to
$R<125$ pc, the structure of the distribution changes considerably, and this
is entirely driven by the inclusion of $\gamma^{2}$ Vel. The probability of
the oldest age bin is still high (90%), but the probability of an age bin
corresponding to a progenitor more massive than $>20M_{\odot}$ increases by
almost an order of magnitude. None the less, the probability of a progenitor
less massive than $20M_{\odot}$ is still 95%. Expanding the sphere still
further primarily increases the probability of the $13.6$-$19.3M_{\odot}$ and
$10.3$-$13.6M_{\odot}$ progenitor bins relative to the $8.1$-$10.3M_{\odot}$
bin. The probability of being less massive than $20M_{\odot}$ is still $90\%$.

In fact, $\gamma^{2}$ Vel is associated with a concentrated cluster of pre-
main sequence stars with a very low velocity dispersion (e.g., Franciosini et
al. 2018). This means that associating the progenitor of Vela with the
formation of $\gamma^{2}$ Vel and its cluster requires the progenitor to have
been a runaway star. To explore this we used the fiducial runaway model from
Renzo et al. (2019) and estimated the distance the surviving stars from
disrupted binaries could travel in the time left for their current
evolutionary phase. Since we now have an isolated star, this will be close to
the remaining life time of the star. Because of mass transfer, we examined the
fraction of stars which could travel more than $100$ pc (the distance from
$\gamma^{2}$ Vel to the pulsar) in bins of either the zero age main sequence
(ZAMS) mass or the mass after the binary is disrupted.

For the ZAMS mass and no restriction on the mass of the star which died to
disrupt the binary, the chances of traveling 100 pc were significant only for
$10$-$15M_{\odot}$ (20%) and $15$-$20M_{\odot}$ (5%) stars. Stars with higher
post-disruption masses could do so because they are initially lower mass stars
that were mass gainers when the primary began to evolve. Even so, only 7.6%,
2.7% and 2.2% of stars with post disruption masses of $20$-$25$, $25$-$30$ and
$30$-$35M_{\odot}$ could travel 100 pc. Demanding a high ($M>30M_{\odot}$)
ZAMS mass primary helps only modestly. Thus, a volume including $\gamma^{2}$
Vel does significantly increase the probability of a higher mass
($\mathrel{\raise 1.29167pt\hbox{$>$}\mkern-14.0mu\lower
2.58334pt\hbox{$\sim$}}20M_{\odot}$) progenitor, but a Vela progenitor formed
in the cluster associated with $\gamma^{2}$ Vel would (a) have to be a runaway
star from a disrupted binary, and (b) would have to have a statistically
unlikely runaway velocity.

##  4 Discussion

The environment of the Vela pulsar is dominated by B stars with only one
nearby O star, $\gamma^{2}$ Vel. If we consider only stars within $R<100$ pc
of Vela, the most likely (95%) age range for its progenitor is
$10^{7.4}$-$10^{7.6}$ years corresponding to a mass range of
$8.1$-$10.3M_{\odot}$. There is clearly a local burst of star formation
associated with this age bin. For these ages and the observed velocity
dispersions of the stars, a radius of $R<100$ pc will not encompass all of the
stars formed within 50-100 pc of the progenitor. Any larger sphere encompasses
the most massive nearby star, the O star plus WR star binary $\gamma^{2}$ Vel.
So for $R<125$ pc, there is still a 90% probability of associating the
progenitor with the $10^{7.4}$-$10^{7.6}$ year age bin, but there is now a 5%
chance of an age allowing masses $>20M_{\odot}$. However, $\gamma^{2}$ Vel is
associated with a very low velocity dispersion cluster, so for the progenitor
to be that massive it would also need to be a runaway star from a disrupted
binary with an unusually high velocity based on the models of Renzo et al.
(2019).

We have assumed single star evolution models in this analysis. We know that
Vela itself was not a binary or triple at the time of its death (Kochanek
2021, Fraser & Boubert 2019) down to companion mass limits $\mathrel{\raise
1.29167pt\hbox{$<$}\mkern-14.0mu\lower 2.58334pt\hbox{$\sim$}}M_{\odot}$. This
does not rule out the progenitor as a merger remnant or as an unbound
secondary from a previous explosion. As a merger remnant or an unbound
secondary that gained significant mass through mass transfer, we have
underestimated the probability of a low mass progenitor because we imposed the
minimum mass/maximum age limit by hand since it cannot be determined by
analyzing a single system. The changes from including the next lower mass bin
would not be huge, because the star formation rate estimated for this next bin
is much lower than for the $10^{7.4}$-$10^{7.6}$ year age bin.

As noted in the introduction, the two fundamental limitations to applying this
method in the Galaxy are estimating the distance to the supernova remnant and
extinction. There are three additional systems which can be analyzed easily.
G180.0$-$01.7 is associated with the radio pulsar PSR J0538$+$2817, which has
a VLBI parallax (Ng et al. 2007, Chatterjee et al. 2009). In Kochanek (2021)
we examined the properties of the nearby luminous stars but did not carry out
a formal analysis as done here for Vela. G205.5$+$0.00.5 (Monoceros Loop) and
G284.3$-$01.8 both contain neutron star high mass X-ray binaries (Hinton et
al. 2009, Corbet et al. 2011) and the stellar companions have Gaia parallaxes
(see Kochanek 2021). Unfortunately, the interacting (probably) black hole
binary SS 433 (see the review by Margon 1984) in G039.7$-$02.0 is likely too
distant ($\varpi^{-1}=8.5$ kpc) to use Gaia parallaxes to trivially select
stars in a $\sim 100$ pc sphere around the binary. The Gaia EDR3 parallax of
SS 433 is also in strong ($\sim 5\sigma$) disagreement with the distances
estimated from kinematic models of the relativistic jets (Blundell & Bowler
2004, Marshall et al. 2013). As discussed in Kochanek (2021), the success of
Cha et al. (1999) in determining the distance to the Vela remnant based on the
distance at which stars began to show absorption features from the supernova
remnant provides a simple observational approach to better determining the
distances to other remnants.

If we consider the 165 pulsars with parallaxes in the ATNF Pulsar Catalog
(Manchester et al. 2005), the vast majority cannot be analyzed using this
method because the formation region cannot be well-localized. In particularly,
the only available age estimate is the spin down age which at best estimates
the time of explosion to within a factor of two. If we require that the
projected distance traveled in the spin down age is less than $100$ pc, only
10 pulsars are left after excluding Vela. Because the proper motions represent
two components of the kick velocity while the motion along the line of sight
contains only one component, the line of sight motions should be less of a
problem. If we assume a typical neutron star kick velocity of $265$ km/s
(Hobbs et al. 2005), the typical line-of-sight distance traveled in the spin
down time for the pulsars with transverse motions less than $100$ pc ranges
from 24 to 152 pc. If we require this line-of-sight motion to be $<100$ kpc
and restrict the parallax distance to be $<2$ kpc, we are left with four
systems other than Vela (J0157$+$6212, J0633$+$1746, J0659$+$1414 and
J2032$+$4127). This would be a sample biased towards supernovae with low kick
velocities.

In the absence of good parallaxes, the alternative is to simply look at all
the stars projected within a fixed projected separation from the SNR over some
broad line of sight distance range consistent with estimates for the distance
to the SNR. This is what is done for all the extragalactic analyses since
there is no possibility of using parallaxes to remove the foreground and
background contamination. However, while some of the external galaxies that
have been examined are highly inclined (e.g., Andromeda at $\sim 13^{\circ}$),
none have as unfavorable geometry for minimizing contamination as we face for
examining Galactic SNRs. Which Galactic SNRs are suitable for such an analysis
will depend on their Galactic coordinates.

Related to this is the question of completeness and contamination in analyses
of stellar populations near extragalactic supernovae and remnants. We can use
our nearly complete knowledge of the environment around Vela to evaluate the
50 pc projected search region used by Jennings et al. (2012) and subsequent
papers. We again restrict ourselves to the more massive stars using the
magnitude and color cuts described in §2. We transform the stellar positions
to axes aligned with Galactic coordinates, select stars in a sphere around
Vela, and then count stars using their positions $10^{7}$ years ago in circles
centered on their median position as if we were looking down on the plane of
the Galaxy (i.e. we ignore the distance of the star from the Galactic plane).
As before, there are 19 stars in a 50 pc sphere around Vela, while a 50 pc
circle centered on their median position $10^{7}$ years ago contains only 9 of
them along with 22 other stars, so the completeness is 47% and the
contamination is 79%. A 100 pc circle contains 16 of the stars along with 140
other stars leading to a completeness of 84% and a contamination rate of 90%.
If instead we consider a 100 pc sphere, which contains 152 stars, a 100 pc
circle contains 73 of them with 73 additional stars for a completeness of 58%
and a contamination rate of 50%. For a 150 pc circle, the completeness and
contamination increase to 76% and 63%, respectively. This age corresponds to
the life time of a $\sim 20M_{\odot}$ star – a $\sim 10M_{\odot}$ star lives
three times longer so the completeness will be far lower and the contamination
will be far higher. For the 50 pc sphere, the completeness/contamination for
the 50 pc and 100 pc circles are 11%/89% and 26%/90%, respectively. For the
100 pc sphere, the completeness/contamination for the 100 pc and 150 pc
circles are 14%/56% and 28%/64%, respectively.

If the environment of Vela is typical, a 50 pc projected search radius
includes almost none of the $\mathrel{\raise
1.29167pt\hbox{$>$}\mkern-14.0mu\lower 2.58334pt\hbox{$\sim$}}5M_{\odot}$
stars that were within 50 pc of the progenitor at birth – the stars within
this radius are overwhelmingly stars which were more distant from the
progenitor at birth. A larger, 100 pc projected search radius does capture a
reasonable fraction of the stars born within 100 pc of the progenitor, but
they are still only about half of the stars within that projected radius. This
strongly suggests that this method of analysis cannot estimate the progenitor
properties for individual supernova remnants or supernovae in external
galaxies– it can only statistically estimate the progenitor properties for
large ensembles of targets. Even with complete three dimensional information,
as we have for Vela, there are significant problems with completeness and
contamination in such small regions on the 30 million year time scale before
$\simeq 10M_{\odot}$ stars explode.

## Acknowledgments

The author thanks K. Auchettl, J. Johnson, M. Pinsonneault, K. Stanek, and T.
Thompson for discussions. CSK is supported by NSF grants AST-1908570 and
AST-1814440. This research has made use of the VizieR catalogue access tool,
CDS, Strasbourg, France (DOI : 10.26093/cds/vizier). The original description
of the VizieR service was published in Ochsenbein et al. (2000). This research
has made use of the SIMBAD (Wenger et al. 2000) database, operated at CDS,
Strasbourg, France.

## Data Availability Statement

All data used in this paper are publicly available.

## References

  * Adams et al. (2017a) Adams, S. M., Kochanek, C. S., Gerke, J. R., et al. 2017a, MNRAS, 468, 4968.
  * Adams et al. (2017b) Adams, S. M., Kochanek, C. S., Gerke, J. R., et al. 2017b, MNRAS, 469, 1445.
  * Aschenbach et al. (1995) Aschenbach, B., Egger, R., & Trümper, J. 1995, Nature, 373, 587
  * Auchettl et al. (2019) Auchettl, K., Lopez, L. A., Badenes, C., et al. 2019, ApJ, 871,
  * Basinger et al. (2020) Basinger, C. M., Kochanek, C. S., Adams, S. M., et al. 2020, arXiv:2007.15658
  * Badenes et al. (2009) Badenes, C., Harris, J., Zaritsky, D., et al. 2009, ApJ, 700, 727.
  * Barlow et al. (2013) Barlow, B. N., Kilkenny, D., Drechsel, H., et al. 2013, MNRAS, 430, 22.
  * Beasor et al. (2020) Beasor, E. R., Davies, B., Smith, N., et al. 2020, MNRAS, 492, 5994.
  * Belczynski et al. (2020) Belczynski, K., Klencki, J., Fields, C. E., et al. 2020, AAP, 636,
  * Blundell & Bowler (2004) Blundell, K. M., & Bowler, M. G. 2004, ApJL, 616, L159
  * Bollig et al. (2021) Bollig, R., Yadav, N., Kresse, D., et al. 2021, ApJ, 915, 28.
  * Bovy et al. (2016) Bovy, J., Rix, H.-W., Green, G. M., et al. 2016, ApJ, 818, 130.
  * Breivik et al. (2020) Breivik, K., Coughlin, S., Zevin, M., et al. 2020, ApJ, 898, 71.
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Burrows et al. (2020) Burrows, A., Radice, D., Vartanyan, D., et al. 2020, MNRAS, 491, 2715.
  * Castelli & Kurucz (2003) Castelli, F., & Kurucz, R. L. 2003, Modelling of Stellar Atmospheres, 210, 20P
  * Cha et al. (1999) Cha, A. N., Sembach, K. R., & Danks, A. C. 1999, ApJL, 515, L25.
  * Chatterjee et al. (2009) Chatterjee, S., Brisken, W. F., Vlemmings, W. H. T., et al. 2009, ApJ, 698, 250
  * Chawla et al. (2021) Chawla, C., Chatterjee, S., Breivik, K., et al. 2021, arXiv:2110.05979
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
  * Corbet et al. (2011) Corbet, R. H. D., Cheung, C. C., Kerr, M., et al. 2011, The Astronomer’s Telegram, 3221,
  * Cousins (1971) Cousins, A. W. J. 1971, Royal Observatory Annals, 7
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, “The IRSA 2MASS All-Sky Point Source Catalog, NASA/IPAC Infrared Science Archive.
  * Davies & Beasor (2020) Davies, B. & Beasor, E. R. 2020, MNRAS, 493, 468.
  * Díaz-Rodríguez et al. (2018) Díaz-Rodríguez, M., Murphy, J. W., Rubin, D. A., et al. 2018, ApJ, 861, 92.
  * Díaz-Rodríguez et al. (2021) Díaz-Rodríguez, M., Murphy, J. W., Williams, B. F., et al. 2021, MNRAS, 506, 781.
  * Dodson et al. (2003) Dodson, R., Legge, D., Reynolds, J. E., et al. 2003, ApJ, 596, 1137
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8
  * Drimmel et al. (2003) Drimmel, R., Cabrera-Lavers, A., & López-Corredoira, M. 2003, AAP, 409, 205
  * Ducati (2002) Ducati, J. R. 2002, VizieR Online Data Catalog
  * Eldridge et al. (2011) Eldridge, J. J., Langer, N., & Tout, C. A. 2011, MNRAS, 414, 3501.
  * Eldridge et al. (2013) Eldridge, J. J., Fraser, M., Smartt, S. J., et al. 2013, MNRAS, 436, 774.
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058.
  * Elitzur & Ivezić (2001) Elitzur, M. & Ivezić, Ž. 2001, MNRAS, 327, 403.
  * Ertl et al. (2016) Ertl, T., Janka, H.-T., Woosley, S. E., et al. 2016, ApJ, 818, 124.
  * Farrell et al. (2020) Farrell, E. J., Groh, J. H., Meynet, G., et al. 2020, MNRAS, 494, L53.
  * Fernández et al. (2018) Fernández, R., Quataert, E., Kashiyama, K., et al. 2018, MNRAS, 476, 2366.
  * Flynn et al. (2006) Flynn, C., Holmberg, J., Portinari, L., et al. 2006, MNRAS, 372, 1149.
  * Folatelli et al. (2016) Folatelli, G., Van Dyk, S. D., Kuncarayakti, H., et al. 2016, ApJL, 825, L22.
  * Franciosini et al. (2018) Franciosini, E., Sacco, G. G., Jeffries, R. D., et al. 2018, AAP, 616, L12.
  * Fraser & Boubert (2019) Fraser, M. & Boubert, D. 2019, ApJ, 871, 92
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91.
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, AAP, 595, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, AAP, 649, A1
  * Garrison & Hiltner (1973) Garrison, R. F. & Hiltner, W. A. 1973, ApJL, 179, L117.
  * Gentile Fusillo et al. (2019) Gentile Fusillo, N. P., Tremblay, P.-E., Gänsicke, B. T., et al. 2019, MNRAS, 482, 4570.
  * Gerke et al. (2015) Gerke, J. R., Kochanek, C. S., & Stanek, K. Z. 2015, MNRAS, 450, 3289.
  * Ghosh et al. (2021) Ghosh, S., Wolfe, N., & Fröhlich, C. 2021, arXiv:2107.13016
  * Gilmozzi et al. (1987) Gilmozzi, R., Cassatella, A., Clavel, J., et al. 1987, Nature, 328, 318.
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., et al. 2019, ApJ, 887, 93.
  * Groh et al. (2013) Groh, J. H., Meynet, G., Georgy, C., et al. 2013, AAP, 558, A131.
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, AAP, 486, 951
  * Hinton et al. (2009) Hinton, J. A., Skilton, J. L., Funk, S., et al. 2009, ApJL, 690, L101
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974
  * Hoffleit & Warren (1995) Hoffleit, D. & Warren, W. H. 1995, VizieR Online Data Catalog, V/50x
  * Jennings et al. (2012) Jennings, Z. G., Williams, B. F., Murphy, J. W., et al. 2012, ApJ, 761, 26.
  * Jennings et al. (2014) Jennings, Z. G., Williams, B. F., Murphy, J. W., et al. 2014, ApJ, 795, 170.
  * Johnson et al. (1966) Johnson, H. L., Mitchell, R. I., Iriarte, B., et al. 1966, Communications of the Lunar and Planetary Laboratory, 4, 99
  * Johnson et al. (2017) Johnson, S. A., Kochanek, C. S., & Adams, S. M. 2017, MNRAS, 472, 3115.
  * Kilkenny et al. (1988) Kilkenny, D., Heber, U., & Drilling, J. S. 1988, South African Astronomical Observatory Circular, 12
  * Kilpatrick et al. (2021) Kilpatrick, C. D., Drout, M. R., Auchettl, K., et al. 2021, MNRAS, 504, 2073.
  * Kochanek et al. (2008) Kochanek, C. S., Beacom, J. F., Kistler, M. D., et al. 2008, ApJ, 684, 1336.
  * Kochanek et al. (2012) Kochanek, C. S., Khan, R., & Dai, X. 2012, ApJ, 759, 20
  * Kochanek (2014) Kochanek, C. S. 2014, ApJ, 785, 28.
  * Kochanek (2020) Kochanek, C. S. 2020, MNRAS, 493, 4945.
  * Kochanek (2021) Kochanek, C. S. 2021, arXiv:2107.13559
  * Koplitz et al. (2021) Koplitz, B., Johnson, J., Williams, B. F., et al. 2021, ApJ, 916, 58.
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231.
  * Lada & Lada (2003) Lada, C. J. & Lada, E. A. 2003, ARA&A, 41, 57
  * Large et al. (1968) Large, M. I., Vaughan, A. E., & Mills, B. Y. 1968, Nature, 220, 340.
  * Lovegrove & Woosley (2013) Lovegrove, E. & Woosley, S. E. 2013, ApJ, 769, 109.
  * Lyne et al. (1996) Lyne, A. G., Pritchard, R. S., Graham-Smith, F., et al. 1996, Nature, 381, 497
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., et al. 2005, VizieR Online Data Catalog, VII/245
  * Margon (1984) Margon, B. 1984, ARA&A, 22, 507
  * Marigo et al. (2013) Marigo, P., Bressan, A., Nanni, A., et al. 2013, MNRAS, 434, 488
  * Marshall et al. (2006) Marshall, D. J., Robin, A. C., Reylé, C., et al. 2006, AAP, 453, 635.
  * Marshall et al. (2013) Marshall, H. L., Canizares, C. R., Hillwig, T., et al. 2013, ApJ, 775, 75
  * Mermillod (1977) Mermillod, J. C. 1977, Geneva: Observatory, 1977
  * Morel & Magnenat (1978) Morel, M. & Magnenat, P. 1978, A&ASup, 34, 477
  * Nadezhin (1980) Nadezhin, D. K. 1980, Ap&SS, 69, 115.
  * Neckel & Klare (1980) Neckel, T. & Klare, G. 1980, A&ASup, 42, 251
  * Nesterov et al. (1995) Nesterov, V. V., Kuzmin, A. V., Ashimbaeva, N. T., et al. 1995, A&AS
  * Neustadt et al. (2021) Neustadt, J. M. M., Kochanek, C. S., Stanek, K. Z., et al. 2021, arXiv:2104.03318
  * Ng et al. (2007) Ng, C.-Y., Romani, R. W., Brisken, W. F., Chatterjee, S., & Kramer, M. 2007, ApJ, 654, 487
  * North et al. (2007) North, J. R., Tuthill, P. G., Tango, W. J., et al. 2007, MNRAS, 377, 415
  * O’Connor & Ott (2011) O’Connor, E., & Ott, C. D. 2011, ApJ, 730, 70
  * Ochsenbein et al. (2000) Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23.
  * Pan et al. (2021) Pan, K.-C., Liebendörfer, M., Couch, S. M., et al. 2021, ApJ, 914, 140.
  * Pastorelli et al. (2020) Pastorelli, G., Marigo, P., Girardi, L., et al. 2020, MNRAS, 498, 3283
  * Jacovich et al. (2021) Jacovich, T., Patnaude, D., Slane, P., et al. 2021, ApJ, 914, 41.
  * Pejcha & Thompson (2015) Pejcha, O., & Thompson, T. A. 2015, ApJ, 801, 90
  * Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, AAP, 500, 501
  * Renzo et al. (2019) Renzo, M., Zapartas, E., de Mink, S. E., et al. 2019, AAP, 624, A66
  * Riley et al. (2021) Riley, J., et al. 2021, arXiv:2109.10352
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161.
  * Sarbadhicary et al. (2017) Sarbadhicary, S. K., Badenes, C., Chomiuk, L., et al. 2017, MNRAS, 464, 2326.
  * Smartt et al. (2009) Smartt, S. J., Eldridge, J. J., Crockett, R. M., et al. 2009a, MNRAS, 395, 1409.
  * Smartt (2009) Smartt, S. J. 2009b, ARA&A, 47, 63.
  * Sukhbold & Adams (2020) Sukhbold, T. & Adams, S. 2020, MNRAS, 492, 2578
  * Smartt (2015) Smartt, S. J. 2015, PASA, 32, e016.
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H.-T. 2016, ApJ, 821, 38
  * Sukhbold et al. (2018) Sukhbold, T., Woosley, S. E., & Heger, A. 2018, ApJ, 860, 93.
  * Thompson et al. (1978) Thompson, G. I., Nandy, K., Jamar, C., et al. 1978, “Catalogue of stellar ultraviolet fluxes: a compilation of absolute stellar fluxes measured by the Sky Survey Telescope (S2/68) aboard the ESRO satellite TD-1”
  * Tokovinin (2018) Tokovinin, A. 2018, VizieR Online Data Catalog,
  * Tonry et al. (2018) Tonry, J. L., Denneau, L., Flewelling, H., et al. 2018, ApJ, 867, 105.
  * van Leeuwen (2007) van Leeuwen, F. 2007, AAP, 474, 653
  * Verbiest et al. (2012) Verbiest, J. P. W., Weisberg, J. M., Chael, A. A., et al. 2012, ApJ, 755, 39
  * Walmswell & Eldridge (2012) Walmswell, J. J. & Eldridge, J. J. 2012, MNRAS, 419, 2054
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9.
  * Wesselius et al. (1982) Wesselius, P. R., van Duinen, R. J., de Jonge, A. R. W., et al. 1982, A&ASup, 49, 427
  * Williams et al. (2018) Williams, B. F., Hillis, T. J., Murphy, J. W., et al. 2018, ApJ, 860, 39.
  * Williams et al. (2019) Williams, B. F., Hillis, T. J., Blair, W. P., et al. 2019, ApJ, 881, 54.
  * Zacharias et al. (2005) Zacharias, N., Monet, D. G., Levine, S. E., et al. 2005, VizieR Online Data Catalog, I/297
  * Zapartas et al. (2017) Zapartas, E., de Mink, S. E., Izzard, R. G., et al. 2017, AAP, 601, A29.