or down of neutron stars has long been attributed to the coupling between the
accreted material and the magnetic field of the neutron star. As described in
Section 2.3, the modalities of the coupling differ between the quasi-spherical
and disk geometries, which both lead to different accretion-induced torques
onto the neutron star magnetosphere. For disk accretion and given the
parameters of Vela X-1, the torque applied to the neutron star computed by
Ghosh & Lamb (1979b) spins it up and produces a spin period derivative of
approximately (Ho et al., 2014; Malacaria et al., 2020)

$\begin{split}|\dot{P}_{\text{disk,+}}|\sim&3.5\text{s}\text{yr}^{-1}\left(\frac{P}{283\,\text{s}}\right)^{2}\left(\frac{M}{2\,M_{\odot}}\right)^{-4/7}...\\\
&...\left(\frac{\dot{M}_{acc}}{2\times
10^{-10}\text{M}_{\odot}\text{yr}^{-1}}\right)^{6/7}\left(\frac{B}{2.6\times
10^{12}\,G}\right)^{2/7}\left(\frac{R_{\text{NS}}}{12\,\text{km}}\right)^{-8/7}\end{split}.$
(12)

For a quasi-spherical geometry, the evolution of the neutron star spin period
depends on the accretion regime. The transition from the direct accretion
regime to the subsonic propeller regime occurs when the mass accretion rate
decreases below $4\times 10^{36}$ erg s-1, which corresponds to the median
X-ray luminosity of Vela X-1 (Shakura et al., 2012). Then, both regimes should
occur in this system. Shakura et al. (2012) set the theoretical framework for
computing the accretion-induced torques applied by a quasi-spherical flow onto
the magnetosphere. Shakura et al. (2018) interpreted the absence of net
spinning-up or -down of the neutron star over the past decades as an
indication that the system must be in spin equilibrium. Based on this
assumption and on the observed spinning-up and -down rates in Vela X-1, they
estimated the value of the dimensionless parameters of their model ,which
encapsulate the physical mechanisms at stake at the outer rim of the
magnetosphere. With the parameters of Vela X-1, it yields the following
spinning-up and -down torques for Vela X-1:

$\cases{&}|\dot{P}_{\text{sph,+}}|\sim
0.1\text{s}\text{yr}^{-1}\left(\displaystyle\frac{v_{r}}{700\,\text{km
s${}^{-1}$}}\right)^{-4}\left(\displaystyle\frac{R_{\text{sh}}}{R_{\text{acc}}}\right)^{2}\\\
&|\dot{P}_{\text{sph,$-$}}|\sim 1\text{s}\text{yr}^{-1}{},$ (13)

respectively, where $R_{\text{sh}}$ is the distance from the bow shock to the
neutron star, which typically is on the order of $R_{\text{acc}}/5$ (El Mellah
& Casse, 2015), $\dot{P}_{sph,-}$ weakly depends on the system parameters, and
only the two strongest dependences are shown for the spin-up rate. These
values are to be compared to the spin period derivatives measured in Vela X-1,
which are about 0.25–0.5 s yr-1 (spin up) and 0.05 s yr-1 (spin down,
Malacaria et al., 2020). In contrast, the quasi-spherical subsonic settling
model predicts higher spin-down than spin-up rates: While the predicted spin-
up rate is on the order of the observed rate, the spin-down rate is
overestimated by more than an order of magnitude. Furthermore, in the moderate
coupling regime, Vela X-1 is expected to accrete. In this regime, a
correlation of the value of the observed spin period derivatives and the X-ray
luminosity is expected but has not been reported. These inconsistencies might
be due both to the numerous unknowns of this model and to the strong
dependence of the spin-up rate on the location of the bow shock with respect
to the accretion radius, which like the relative wind speed is unknown to
within at least a factor of 2. Finally, in the supersonic accretion regime of
the quasi-spherical settling model, the neutron star is generally spun up by
accretion-induced torques. Spinning down can only occur if the neutron star
spin axis is upside down with respect to the orbital angular momentum axis, or
if clumps carry enough angular momentum to momentarily cause an inversion of
the sign of the accreted net angular momentum.

###  4.3 Emission from the accretion column or polar cap

Most of the hard X-rays we observe are produced very close to the neutron star
within the accretion column. In this column, the hot in-falling plasma can
upscatter through Compton-interaction seed photons to very high energies that
roughly resemble the observed continuum spectrum (Davidson, 1973; Davidson &
Ostriker, 1973; Meszaros & Nagel, 1985a; Rebetzky et al., 1988). However,
because the problem is complex, these early models had to make many
simplifying assumptions, such as reducing the problem to 1D or neglecting the
accretion flow dynamics. When they escape the accretion column, the photons
cool the plasma and prevent a full thermalization of the flow (Arons et al.,
1987). The source of the seed photons is still debated: A large part very
likely is due to cyclotron radiation (Nagel, 1981; Bussard et al., 1985; Arons
et al., 1987), but thermal photons from the accretion mound on the neutron
star surface will also play a role (Becker & Wolff, 2005a).

While this basic idea of Comptonized photons agrees roughly with the observed
spectra, detailed modeling of the spectrum produced in the accretion column
continues to be difficult even with modern computers. These difficulties
mainly arise because for a correct description of the particles and their
interaction within the column, a full quantum-mechanical as well as
relativistic treatment is necessary. In particular, the electron cross-
sections are strongly polarization dependent and are increased at the
cyclotron energy by orders of magnitude (e.g., Daugherty & Harding, 1986, see
also Sect. 4.4), which is computationally expensive to model. Additionally,
the open question remains of where and how the in-falling plasma is
decelerated (e.g., Becker & Wolff, 2005a).

In the past decade, Becker & Wolff (2005a, b, 2007); Becker et al. (2012) have
presented a more detailed model of the spectral formation in the accretion
column, taking bulk (first-order Fermi acceleration) and thermal
Comptonization into account and injecting seed photons from bremsstrahlung,
cyclotron radiation, and thermal contribution. This model is able to reproduce
the spectral shape of a few selected X-ray binaries, but requires a simplified
description of the velocity in the accretion column and does not account for
the variability of theses sources (Wolff et al., 2016). It assumes a high
luminosity, in which a radiation-dominated shock is formed above the neutron
star surface and decelerates the in-falling plasma to subsonic speeds.

Farinelli et al. (2012, 2016) expanded on the model by Becker & Wolff (2007),
allowing for different velocity profiles and a magnetic field gradient within
the accretion column. Their model allowed them to report a satisfying
description of the X-ray spectra of bright sources such as Her X-1 and 4U
0115+63.

Vela X-1 only rarely exhibits luminosities above $10^{37}$ erg s-1 . It is
therefore unclear whether a radiation dominated shock forms. This
intermediate-luminosity range is in particular hard to model because it is
unclear which approximations can be made. Vela X-1 also shows a complex pulse
profile, which might indicate that both fan and pencil beam components are
present. The most recent update of the Becker & Wolff models was presented by
(West et al., 2017a, b), but the model has not yet been applied to Vela X-1.

The observed pulse profile strongly depends on the assumed emission profile of
the accretion column and of the polar cap, and also on light-bending effects
around the neutron star. For the emission geometry, two simple approximations
are often made: a slab, or a column geometry. In the slab geometry, the
accretion mainly escapes from the sides of the accretion column (and
perpendicular to the magnetic field) close to the neutron star surface,
creating a so-called fan-beam pattern (i.e., in a similar configuration as a
polar cap). In a column geometry, the emission escapes along the column and
parallel to the magnetic field, forming a pencil beam, with a certain (narrow)
opening angle (Meszaros & Nagel, 1985b; Meszaros & Riffert, 1988).

While these two different geometries predict significantly different pulse
profiles (Meszaros & Nagel, 1985b; Wang & Welter, 1981), they prove difficult
to compare with observations. The observed profiles are highly complex and
strongly energy dependent (see Sect. 3.1.9), and this makes them far more
complicated than simple models predict.

Over the years, various efforts were made to offer physical descriptions of
the pulse profile shapes. Early studies were concerned with the physical
production region of the X-rays and their emission geometry. They found that
the models were able to match the data for Vela X-1 approximately (e.g., Wang
& Welter, 1981; Leahy, 1990; Sturner & Dermer, 1994). However, the omission of
gravitational light-bending led to unrealistic estimates of the emission
geometry (Riffert et al., 1993; Leahy & Li, 1995; Kraus et al., 1995).

More detailed work led to descriptions of profiles for single sources, using a
large number of assumptions (Kraus et al., 1995; Caballero et al., 2011;
Falkner, 2018). Recently, models have been presented that allow the user to be
very free in the choice of geometry and emission pattern, and they no longer
require many of the previous assumptions (Cappallo et al., 2017; Falkner,
2018). However, the sheer number of free parameters and large degeneracy
between them has so far made it impossible to identify unique solutions. Some
of the most relevant physical parameters are the amount of relativistic
boosting of the in-falling plasma, the size and mass of the neutron star, and
the configuration of the magnetic field. For Vela X-1, no recent physical
description of the pulse profile has been put forward at the time of this
review.

###  4.4 Cyclotron resonant scattering features

As explained in Section 2.4, the CRSF features that are visible in broadband
X-ray spectra of Vela X-1 (see Sect. 3.1.7) are caused by resonant scattering
of X-ray photons at the corresponding energies on electrons with quantized
motion perpendicular to the magnetic field. Calculating the scattering cross-
sections requires detailed fully relativistic QED-based calculations. Over the
decades, such calculations have been undertaken by various authors and have
included increasingly complex details (e.g., Harding & Daugherty, 1991; Sina,
1996; Isenberg et al., 1998a, b; Araya & Harding, 1999; Schönherr et al.,
2007). The most recent treatment based on elaborate Monte Carlo calculations
is presented in Schwarm et al. (2017b, a).

The scattering cross-sections are very strongly peaked close to the Landau
level energies, up to several orders of magnitude compared to Thomson
scattering. This makes the plasma optically very thick for photons at these
energies. Because the the energy of the electrons is only quantized in the
direction perpendicular to the magnetic field, however, the angle between the
magnetic field and the photon becomes relevant. For large angles, the cross-
sections become thermally broadened and shift to higher energies for the
harmonic levels. An open question in CRSF research is the fact that model
calculations (Araya & Harding, 1999; Schwarm et al., 2017b, a) tend to predict
asymmetrical lines, frequently showing emission wings at energies below and
above the central energy, while observed features tend to be broad, without a
marked asymmetry, and without a clear sign of the predicted emission features.

In order for discrete CRSFs to be observable, the sample magnetic field has to
be confined to a very narrow range, indicating a closely confined region
within the accretion column, possibly a shock region in the column or close to
the poles (Becker & Wolff, 2007, and references therein). The often observed
relatively broad and shallow CRSFs might arise because multiple line-forming
regions contribute (Nishimura, 2008). Poutanen et al. (2013) proposed another
explanation according to which, CRSFs are formed due to reflection of the
downward-beamed radiation from the accretion column on the neutron star
surface around the poles.

The strong angular dependence of the line features predicted in theoretical
calculations would in principle lead to another diagnostic for possible
emission geometries from the pulse phase, and thus geometrically, resolved
spectroscopy. This is limited by the fact that because of gravitational light-
bending, we can expect that multiple zones of the accretion column or neutron
star surface contribute at least most of the time (Falkner, 2018). We would
therefore also need to have a good grasp on separating these mixed
contributions.

###  4.5 Variability of the observed X-ray emission

In variability studies, the two privileged simulation outputs to be compared
to observations are occurrence diagrams and power spectral distributions.
While the first show the fraction of time that a variable spends at a given
value (e.g., the mass or angular momentum accretion rates, or the column
density), the second retains the information and reveals coherence timescales
and possible periodicities.

####  4.5.1 Variable absorption along the line of sight

The mass-loss rate from the donor star is high, therefore absorption along the
line of sight from the neutron star to the observer is high but is also
variable (Section 3.1.3). Because we observe Vela X-1 close to edge-on, we
expect a periodic modulation of the column density as the neutron star orbits
the blue supergiant. The average orbital $N_{\textrm{H}}$ profile grants
access to the orbital parameters and to the stellar mass-loss rate, provided
estimates are available for the wind speed and the stellar radius. The origin
of the stochastic variability of the measured absorbing column density at a
given orbital period is threefold. It can be due either to changes in the
amount of material that is integrated along the LOS because of clumps in the
stellar wind or because of a change in the flow structure near the neutron
star magnetosphere, or modifications in the opacity of the medium can cause it
through changes in the wind ionization structure.

Regardless of whether they participate in the accretion process, clumps
intercepting the LOS will produce a variability in column density on short
timescales. Oskinova et al. (2012) designed a 2D model of massive ($3\times
10^{21}$g) and intermediate-size clumps (2% of the stellar radius at one
stellar radius above the stellar photosphere) that propagate radially from the
donor star. They computed the extinction coefficient as a function of time and
showed how clumps that are compressed in the radial direction tend to produce
a higher level of variability than their spherical counterparts (see also
Oskinova et al., 2006). Grinberg et al. (2017) estimated the characteristic
amplitude and timescale of absorption episodes with a simple model of clumps,
to be compared to column densities derived from hardness monitoring with time-
resolved spectroscopy at orbital phase $\phi\sim 0.25$. This motivated a
comprehensive exploration of the $N_{\textrm{H}}$ variability induced by a
clumpy wind: El Mellah et al. (2020b) designed a 3D model of clumps flowing
radially away from the star with different velocity profiles, orbital
parameters, and clump properties. Vela X-1 was one of the two HMXBs to which
the model was applied. The authors computed the median $N_{\textrm{H}}$
orbital profile over many orbital periods, along with the standard deviation
and coherence timescale at each orbital phase. The timescale was found to be
associated with the flyby time of the smallest clumps along the LOS, which
granted access to their spatial extent. Based on the porosity length, the
authors also showed that the mass of the clumps could be constrained based on
the standard deviation of the column density: Because they are rarer, higher-
mass clumps produce a higher spread of the measured column density at a given
orbital phase.

Manousakis et al. (2012) computed the structure of a smooth stellar wind in
hydrodynamics numerical simulations in the orbital plane. Although their model
was tailored to another high-mass X-ray binary and not to Vela X-1, they were
able to reproduce the excess of column density at orbital phases 0.4-0.8
reported in Manousakis & Walter (2011) for IGR J17252-3616 and traced its
origin back to an overdense accretion wake trailing the neutron star. In Vela
X-1, Malacaria et al. (2016) measured episodic enhanced absorption events near
inferior conjunction. They performed a detailed spectral analysis that showed
that the data could be attributed to an inhomogeneous accretion wake, in
agreement with the idea proposed by Doroshenko et al. (2013). The eccentricity
of the orbit also introduces a minor asymmetry between the evolution of the
column density when the neutron star moves toward and away from Earth.

Finally, intrinsic variations of the X-ray ionizing emission from the
immediate vicinity of the neutron star (see Section 4.5.2) alter the wind
ionization structure. The highly ionized fraction of the wind, near the
neutron star and in low-density regions, is transparent to X-rays and does not
contribute to the absorbing column density. It induces $N_{\textrm{H}}$
variations even if the amount of material integrated along the LOS remains
unchanged.

While off-states have usually been explained as due to a quenching of the
accretion itself (see Section 4.5.2), one model found that the lowered
luminosity during off-states might be due to scattering of hard X-rays along
the line of sight: Sidoli et al. (2015) reported different off-state
occurrence rates at orbital ingress and egress that they interpreted as an
indication that off-states were due to X-rays that were being scattered in the
photoionization wake identified by Blondin et al. (1990) in numerical
simulations.

####  4.5.2 Variability in mass accretion rate

The dominant mechanisms responsible for the intrinsic variability in mass
accretion rate in Vela X-1 have not yet been conclusively identified. Periodic
changes are expected because the eccentricity of the neutron star orbit leads
to non-negligible systematic differences in the mass density of the
environment in which the neutron star is embedded: Notwithstanding wind
acceleration and nonisotropic effects, the density at periastron would be 40%
lower than the density at apastron. This number is a lower limit because the
wind is still accelerating at the orbital separation and is more strongly
beamed toward the neutron star at periastron, when the filling factor is the
highest (see Section 4.2.1). Concerning the stochastic variability of the mass
accretion rate onto the neutron star surface, the very location of its origin
varies from one model to the next.

The X-ray ionizing feedback has long been suspected to modulate the mass
supply onto the neutron star in an unstable manner: An increase in X-ray
luminosity leads to a slower wind because the line acceleration is inhibited,
which increases the mass accretion rate (see Section 4.2.1). The cycle closes
when wind acceleration is so inhibited that stellar material no longer
penetrates the Roche lobe of the neutron star. In an attempt to quantify this
effect, Ho & Arons (1987b) designed a 1D model along the line that joins the
two bodies with a prescribed smooth and isotropic wind. They equated the
accretion X-ray luminosity associated with the BHL formula (4) with the X-ray
luminosity corresponding to the extent of the region around the neutron star
in which the ionization parameter is above $\xi_{\text{crit}}=10^{4}$ erg cm
s-1. Using the orbital parameters of Vela X-1, they obtained two solutions. In
the first, the wind was hardly affected by the X-ray ionizing feedback and
yielded an X-ray luminosity corresponding to the one observed in Vela X-1.
This solution was later found by Ho & Arons (1987a) to be prone to
fluctuations with realistic timescales (see, e.g., Kreykenbohm et al., 2008).
In the second solution derived by Ho & Arons (1987b), the X-ray luminosity was
about four orders of magnitude higher, with a wind speed significantly lowered
by line-acceleration quenching. In spite of episodic flares in which Vela X-1
reaches a few 1037erg$\cdot$s-1, the system has never been observed in such a
high state, however. Karino (2014) performed a similar computation, but
accounted for optical depth and a gradual effect of the X-rays on line
acceleration. Karino (2014) retrieved a bimodal behavior, but in contrast to
Ho & Arons (1987b), Vela X-1 was lying on the branch corresponding to a high
luminosity, associated with a slow wind. Manousakis & Walter (2015a) revived
the role of the X-ray ioninzing feedback that they found to be a possible
origin of the observed off-states. In 2D hydrodynamic simulations in the
orbital plane, they found that when wind acceleration was inhibited above a
critical ionization parameter of $10^{2.5}$erg$\cdot$cm$\cdot$s-1, a permanent
X-ray luminosity of $4\times 10^{36}$erg$\cdot$s-1 could produce an
oscillation of the position of the accretion front shock. This breathing mode
induces a subsequent modulation of the mass accretion rate onto the neutron
star, with $\sim$30-minutes-long off-states during which the accretion is
quenched, separated by a characteristic duration that agrees reasonably well
with the characteristic timescale of the X-ray emission from Vela X-1 reported
in Kreykenbohm et al. (2008). Recently, Bozzo et al. (2020) designed a
detailed semianalytical model of a smooth spherically symmetric wind and
showed that accounting for the eccentricity of the orbit in a simplified model
of accretion leads to an average mass accretion rate almost three times higher
near periastron than near apastron. Because the accretion wake intercepts the
line of sight at orbital phases 0.4-0.8, precisely when the neutron star is
near apoastron, it is difficult to separate between a decrease in intrinsic
emission and an increase in absorption, even using hard X-ray observations.
The authors also included the X-ray ionizing feedback and found that its
effect was limited to the vicinity of the neutron star. This challenges the
conclusions drawn by Manousakis & Walter (2015a) about the origin of off-
states, but agrees with Sander et al. (2018).

Because its structure is inherently clumpy, the line-driven wind from the
donor star provides material to the neutron star at variable rates. In a
seminal study, Ducci et al. (2009) designed a statistical model of clump
capture in which they accounted for this component by assuming that each clump
intersecting the disk of radius $R_{\text{acc}}$ around the neutron star is
instantaneously captured, with a correction factor if the overlap is only
partial. They found occurrence diagrams consistent with what is observed in
Vela X-1. However, these results were challenged by 3D hydrodynamic
simulations of the accretion of a clumpy wind. El Mellah et al. (2018) showed
that clumps were strongly compressed, distorted, and eventually mixed within
the shocked region. The bow shock that forms around the neutron star
significantly decreases the peak-to-peak variability compared to what is
expected in the ballistic BHL framework, in which clumps are instantaneously
accreted without accounting for their thermodynamical properties. These
results indicate that the intrinsic variability of the X-ray emission cannot
be considered a direct tracer of the wind microstructure and that other
mechanisms are likely to dominate the intrinsic variability of the emission.
We note that in simulations of accretion onto supermassive black holes
mediated by a centrifugal gating mechanism, Gaspari et al. (2015) found
similar properties for the variability of the mass accretion rate (e.g., a
log-normal occurrence diagram, indicative of self-organized criticality).

Although realistic 3D clumps per se do not produce a variability high enough
to account for the changes in intrinsic X-ray luminosity, they can play a role
in modifying the conditions of accretion at the neutron star magnetosphere.
Due to the proximity between the magnetospheric and circularization radii in
Vela X-1 (see Figure 2), the accretion flow between the shock and the
magnetosphere in Vela X-1 might transit between a disk-like and a spherical
geometry, depending on the instantaneous and local wind properties. In each
geometry, there are regimes of low and high mass accretion rates, depending on
the stellar material supplied at the orbital scale (see Section 4.2.1).
Transitions between regimes and/or geometries might contribute to the overall
variability. Bozzo et al. (2016) investigated this scenario in a model of
accretion of a clumpy wind accounting for the magnetic and centrifugal gating
mechanisms. They showed that for orbital and neutron star spin periods close
to the ones in Vela X-1, switches between regimes were occurring. Even if the
authors warned the readers about the limitations of their model, which is uni-
dimensional and does not solve the wind dynamics, their work highlighted the
importance of considering the coupling between the plasma and the
magnetosphere to fully appreciate the evolution of the intrinsic X-ray
variability.

###  4.6 Evolutionary scenario

High-mass X-ray binaries are an important evolutionary stage on the way to
form compact object binaries and thus an interesting benchmark to understand
the origin of gravitational wave sources. With a neutron star being the
compact object, Vela X-1 is a candidate for either a black hole-neutron star
merger or a double neutron star merger, assuming the system remains bound
after the current supergiant component has undergone core collapse. The
identification of the evolutionary scenario for an individual system such as
Vela X-1 is complicated. Given the high number of parameters in binary
evolution and the different possible channels, a snapshot of the present
status can thus only be used to rule out certain scenarios rather than
pinpointing at a particular evolutionary track.

Neglecting any scenarios that could involve higher-order systems, the standard
scenario to form an HXMB considers two massive stars with an uneven, but not
too extreme mass ratio (e.g. Tauris & van den Heuvel, 2006). The more massive
primary star evolves faster and eventually fills its Roche Lobe, which will
lead to mass transfer from the primary to the secondary (e.g. Paczyński, 1971;
van den Heuvel, 1976). Depending on the nuclear burning status of the primary,
this is referred to as either case A (during central H-burning), case B (after
central H burning, but before central He burning), or case C (after the start
of central He burning) mass transfer. For typical separations and assuming
close mass ratios, this first mass transfer period in the system is generally
assumed to be stable (e.g., Podsiadlowski et al., 1992; Ritter, 1999; de Mink
et al., 2008). The mass transfer will stop after a large fraction of the
hydrogen envelope has been removed from the primary. The secondary has now
gained hydrogen-rich matter. Depending on the precise parameters, the
secondary could now even evolve faster than the primary and reach core
collapse first (e.g. Pols, 1994). However, given that the mass transfer might
have happened in an already advanced stage of the primary, the primary as the
first object to undergo core collapse is the more common outcome. If the
system remains bound beyond this stage (Blaauw, 1961), we are left with a
compact object resulting from the primary and a secondary which has gotten
additional material and – depending on the mixing timescale of this “new”
material into the deeper layers – has been “rejuvenated” (e.g. Braun & Langer,
1995). Once the compact object is able to accrete matter from the secondary –
either by wind accretion or Roche Lobe overflow (RLOF) – the system will
appear as an HMXB.

In the case of Vela X-1, the compact object is a neutron star and the
secondary has a current mass on the order of $20\,M_{\odot}$. Consequently, we
would expect a primary more massive than that with suggestions going up to
$60\,M_{\odot}$ (Quast et al., 2019). Given that enough mass needed to be
removed to eventually yield a neutron star as the compact object, the latter
number might be a bit high, but illustrates the severe uncertainty when trying
to extrapolate the systems’ backstory. Even taking into account the
accompanying supernova during the formation of the neutron star, the core mass
of the primary had to be low enough to avoid the formation of a black hole,
which means that either the primary mass had to be low enough or the mass
transfer was efficient enough to alter the evolution before a more massive
core could be formed (e.g. Wellstein & Langer, 1999).

The near-circular orbit of Vela X-1 prompts the questions whether the first
mass transfer in the system might have lead to a Common Envelope (CE) phase
(Paczynski, 1976). Even though many uncertainties remain on the detailed
physical processes involved in CE evolution, their occurrence is typically
inferred from an evolutionary necessity (e.g. Ivanova et al., 2013). CE stages
are thought to occur for unstable RLOF and help to shrink and circularize the
orbit, at least if the system can eventually eject the envelope and avoid
merging. However, CE stages are unlikely to occur for the first mass transfer
stage in an HMXB (Tauris & van den Heuvel, 2006). The parameters for Vela X-1,
including its “runaway” nature (Kaper et al., 1997, see also Figure 14),
indeed point more towards this standard scenario (van den Heuvel, 2019).

The remaining secondary in Vela X-1 is classified as “overluminous” (e.g.
Wickramasinghe, 1975; Conti, 1978). This term generally refers to that fact
that the object is more luminous than expected for its mass, although the
details differ in whether this refers to a simple comparison with single star
evolution (e.g. Quast et al., 2019) or is also used for remaining
discrepancies in the context of binary evolution scenarios (e.g. Vanbeveren et
al., 1993). Regardless of terminology, it is important to account for the fact
that the secondary has been rejuvenated and thus will usually have a different
luminosity and temperature compared to a single star of the same spectral
type.

The recent atmosphere analysis of HD 77581, the donor of Vela X-1, yields
$\log L/L_{\odot}\approx 5.5$ (see Table 5). Relying on stellar structure
calculations by Gräfener et al. (2011), a single star on the zero age main
sequence with this luminosity would have a mass of approximate
$47\,M_{\odot}$. This is twice as large as measured (cf. Sect. 5.3). The other
extreme would be a He star with a thin hydrogen layer negligible in mass. In
this case, the donor star in Vela X-1 would be about $15.5\,M_{\odot}$. The
measured masses of about $20\,M_{\odot}$ are much closer to the second case
and thus illustrate that the now donor star has probably evolved considerably
before gaining mass from the primary. Vanbeveren et al. (1993) deduced that
the high ratio between luminosity and mass ($L/M$-ratio) could only be
explained if the central hydrogen abundance of the secondary was already down
to a mass fraction of $0.1$ when it started accreting from the primary.
Moreover, the star would need to have been fully convective during the
accretion stage and its surface should now be enriched in nitrogen, but
depleted in oxygen and carbon. The latter is indeed in line with the recent
analysis by Giménez-García et al. (2016), although their determined surface
hydrogen abundance of $X_{\text{S}}=0.65$ is slightly higher than the range
suggested by Vanbeveren et al. (1993). However, depending on the details of
the accretion into the secondary, the resulting values of $X_{\text{S}}$ can
vary quite a bit (Hellings, 1983; Braun & Langer, 1995) with existing
scenarios not excluding the measured abundance. In any case, the surface
hydrogen abundance of the donor in Vela X-1 is below the solar value, thus
giving rise to the label that the donor of Vela X-1 is “He enriched”.

The high $L/M$-ratio and the still rather high surface hydrogen abundance both
lead to a stronger stellar wind mass loss than assumed for an isolated star of
the same spectral type as stellar wind mass is significantly boosted for
higher values of the so-called “Eddington parameter” $\Gamma_{\text{e}}\propto
q_{\text{ion}}\,L/M$ (e.g. Gräfener & Hamann, 2008; Vink et al., 2011; Sander
et al., 2020; Sander & Vink, 2020). A higher rotation rate due to angular
momentum transfer in the earlier stages of the evolution would add to this.
Consequently, the wind mass-loss rate of the secondary is now higher than it
ever was for the secondary in the system before the onset of mass transfer.

In the standard picture, the HMXB type with a supergiant donor was seen as a
relatively short-lived stage (e.g. Tauris & van den Heuvel, 2006) due to the
unstable nature of the mass transfer arising from the – now – high mass ratios
between the donor and the compact object and due to the convective envelope of
the donor. This claim has recently been challenged by Quast et al. (2019), who
argue that despite the shrinking of the orbit also the donor radius is
decreasing as the He surface abundance increases due to the mass transfer.
This scenario would significantly prolong the HMXB stage and provide a better
explanation for the relatively large number of HMXBs hosting a blue supergiant
donor star in our own Milky Way. Depending on the $L/M$-ratio (Sander & Vink,
2020), the donor could eventually appear as Wolf-Rayet star (Quast et al.,
2019). Eventually, a CE stage with the compact object seems to be unavoidable
(Hjellming & Webbink, 1987; Tauris & van den Heuvel, 2006), which would either
lead to shrinking of the orbit or the eventual merger of the two stars after
the formation of a so-called “Thorne-Żytkow object” (TZO, Thorne & Zytkow,
1975, 1977). Recently, Oskinova et al. (2018) suggested that sgB[e] HMXBs
could represent those HMXBs currently in or shortly after a CE stage. A
successful ejection of the CE would remove further mass from the system,
leading to a He star even if the considerations by (Quast et al., 2019) were
not applicable. However, a system like Vela X-1 might not survive the CE phase
as van den Heuvel et al. (2017) argue that all systems with mass ratios larger
than approximately $3.5$ will lead to an inspiral of the compact object. In
contrast, recent simulations suggest that the parameter space leading to
successful CE ejection is more complex. Depending on the properties of the
donor star envelope, CE evolution can be immediately followed by stable mass
transfer (Klencki et al., 2020). So far, there are a lot of remaining
uncertainties on the scenarios and detailed evolution calculations are
missing, but given that the current mass ratio for Vela X-1 is about $10$, the
system would be a prime candidate for this fate and thus never form a compact
object binary. Yet, the observational confirmation of this scenario is
hampered by the severe uncertainties for the mass loss of lower-mass He stars
(e.g. Vink, 2017; Sander et al., 2020). While traditional considerations often
associated He stars with a WR-type spectral appearance, this is only true if
the stars have a sufficiently high $L/M$-ratio (Shenar et al., 2020; Sander &
Vink, 2020). For the case of Vela X-1, this would mean that the current donor
would have to lose about $5\,M_{\odot}$ – while at least keeping the current
luminosity – before reaching a sufficient $L/M$-ratio for a WR-type mass-loss
rate. Hence, the evolutionary fate of Vela X-1 is quite uncertain. However,
the most expected outcome according to current considerations and known system
parameters is a merger of two components in a CE stage (e.g., Belczynski et
al., 2012). The inspiral of the neutron star into the core of the now
secondary will eventually lead to the formation of a black hole. This could
give rise to a less common type of short gamma-ray bursts (Fryer et al.,
2013).

In addition to our knowledge of its evolutionary track, Vela X-1 was shown to
move at high speed ($\sim$90 km s-1) with respect to its local environment.
This is likely indicative of the natal kick produced by the asymmetric
supernova explosion when the neutron star formed, probably a few million years
ago (see Sect. 5.1.2). A bow shock associated with Vela X-1 was found by Kaper
et al. (1997) based on its H$\alpha$ emission. Although its main axis is
consistent with the relative motion of Vela X-1 through the interstellar
medium, Gvaramadze et al. (2018) identified nonsymmetric features in the
vicinity of the shock, such as an elongated H$\alpha$-emitting structure.
Based on magnetohydrodynamic simulations, they showed that the opening angle
of the bow shock and its H$\alpha$ brightness could be reproduced by
considering the interaction of a region of higher density in the interstellar
medium with the stellar wind from the donor star in Vela X-1. Interestingly
enough, they were able to reproduce the emission from the filamentary
structure only by accounting for a helical magnetic field in the stellar wind
that they tentatively attributed to the interaction between the stellar wind
and the magnetosphere of the neutron star.

###  4.7 Taxonomy of wind models

Relevant models to describe the physics at stake in Vela X-1 have been
published over the past decades. Most of them resorted to the numerical tool
to solve the flow dynamics , but a handful are essentially based on
(semi-)analytical approaches to study the X-ray ionizing feedback (Hatchett &
McCray, 1977), the flow structure (Foglizzo & Ruffert, 1997), the clump
capture rate (Ducci et al., 2009), the accretion mechanism (Bozzo et al.,
2016), the formation of wind-captured disks (Karino et al., 2019), or the
secular spin evolution of the neutron star (Karino, 2020), for instance.
Democratization of access to computing facilities and exploratory studies led
to an overall increase in the complexity of the models that were developed to
provide more detailed explanations of the physical mechanisms in Vela X-1. As
an unfortunate side effect, the numerous physical and, when applicable,
numerical assumptions required to capture the physics and solve the underlying
equations preclude any one-to-one comparison of results, even in the rare
cases when studies focus on similar aspects of the system: The price to pay
for more realistic models has been an accelerated division of the system
description into a constellation of models that each address more specific
aspects in more detail. When observational data are not enough to decide, the
robustness of the provided scenarios does not lie in the strict
reproducibility of the results, but instead in the repeated occurrence of
qualitatively similar features from one simulation to the next and in the
compatibility of an explanation with the models proposed by other teams.

As an illustration of this difficulty of comparing different frameworks
although they all solve the equations of radiative-hydrodynamics in 1D,
Watanabe et al. (2006), Krtička et al. (2012), and Krtička & Kubát (2017),
Krtička et al. (2018) and Sander et al. (2018) did not obtain the same wind
ionization and velocity profile (see Section 4.1.2). The models developed by
the last two teams account for more physics than in the more simplified
calculation by Watanabe et al. (2006), but they both chose to rely on
different assumptions concerning the frame in which the radiative transfer is
performed, the ions involved in the computation, the X-ray emission, and the
stellar properties or the wind clumping, among others.

Similarly, the 2D and 3D simulations ran by Čechura & Hadrava (2015), which
are primarily tailored to the properties of Cygnus X-1, another wind-fed high-
mass X-ray binary, Blondin et al. (1990), Manousakis et al. (2012), Manousakis
& Walter (2015a), Manousakis & Walter (2015b) and El Mellah et al. (2019a) did
not handle the effect of the X-rays on wind acceleration in the same way.
Because of their multidimensional aspect, none of these simulations can afford
to treat radiative transfer in a way as realistic as the studies above. For
instance, they cannot capture any spectral property, but instead have to focus
on the overall photoionizing flux, which is taken to be the X-ray flux from
the neutron star. In order to represent the wind acceleration quenching, they
relied on the critical ionization parameter (see equation (2)) beyond which
the ions contributing most to the resonant line absorption are assumed to be
so depleted that acceleration drops to zero. In this all-or-nothing approach,
it can only be guessed which critical ionization parameter is realistic
because the dominant ions are unknown and sharply depend on the stellar
properties and chemical composition. Even in this case, the sudden quenching
assumed by Manousakis & Walter (2015b) or the even steady decrease of the
acceleration assumed by Čechura & Hadrava (2015) has been shown to be a far
too simplified approximation (Sander et al., 2018). Instead, El Mellah et al.
(2019a) relied on the 1D line acceleration profile computed by Sander et al.
(2018) for the donor star in Vela X-1, but assumed that the nonisotropic 3D
features induced by the orbital motion in their simulations had no significant
effect on this profile, but were purely a function of the distance to the
donor star. This hypothesis remains to be confirmed. Below the critical
ionization parameter (100.5 erg$\cdot$cm$\cdot$s-1 in Čechura & Hadrava 2015,
102.5 erg$\cdot$cm$\cdot$s-1 in Manousakis & Walter 2015b), line acceleration
is computed differently in these studies, leading to different amounts of
specific kinetic energy that was deposited in the wind by the stellar photon
field when the flow reached the neutron star.

The wind microstructure has also been treated in different ways: While some
works assumed a smooth wind (Watanabe et al., 2006; El Mellah & Casse, 2017;
Blondin et al., 1991), others included wind clumpiness in their models or
simulations. The clumps are either indirectly represented by their effect on
the radiation field (Sander et al., 2018), are parameterized as spheres of a
given mass and radius (Oskinova et al., 2012; El Mellah et al., 2020b), or are
computed from 2D radiative-hydrodynamics simulations on pseudo-planar grids
before being extended to 3D (El Mellah et al., 2018). Because of their small
sizes, however, including more realistic clumps comes at the expense of a
smaller simulation space in radiative-hydrodynamics simulations, or of
simplified equations of motion of the flow.

A detailed taxonomy of the models of Vela X-1 in particular and of wind-fed
neutron stars in supergiant X-ray binaries in general is beyond the scope of
this paper. It would help to determine which model could prove fruitful in
order to interpret a given set of observations, however.

##  5 Properties of the Vela X-1 system

In this section we compile the existing knowledge on the essential parameters
of the Vela X-1 system and how they have been obtained. We also indicate where
results were strongly based on parameters from previous work. The section is
broken down into subsections for different essential parameters, but we note
that depending on the data and analysis methods that were used, various
parameters have frequently been derived jointly or cannot be determined in an
independent fashion from other parameters. A trivial example are the masses of
the binary partners and inclination (Figure 17).

###  5.1 Distance from Earth and relation to the Vela OB1 association

In order to derive the real luminosities of the supergiant and the X-ray
source, it has been the goal for many different authors to determine the
distance to Vela X-1. Figure 13 illustrates the evolution of different
distance estimates over time.

Figure 13: Upper panel: Evolution of distance estimates to the Vela X-1 system
over time derived with different approaches, as explained in the text. At the
bottom we also show the distance ranges derived for stars in the Vela OB1
association and our new distance determination based on the new Gaia-EDR3
data. Lower panel: Posterior distribution of our updated Gaia distance
estimate.

####  5.1.1 Evolution of the distance estimates

The earliest distance estimate to the Vela X-1 system, or more precisely, to
HD 77581, was reported by Hiltner et al. (1972), who derived a distance of 1.2
kpc (distance modulus 10.5) based on an assumed absolute magnitude of
$M_{V}=-6.0$ and a ratio of 3.0 for the total-to-selective absorption. This
distance estimate was still used as a baseline in the 1980s, for instance, by
Nagase et al. (1984b). Based on UV observations with the European satellite
Thor-Delta 1A (TD1), Nandy et al. (1975) arrived at a similar result of
$1.3\pm 0.2$ kpc.

In contrast, Wickramasinghe et al. (1974) derived a distance $>$2 kpc by
comparing the reddening of HD 77581 with that of stars in its vicinity.
Zuiderwijk et al. (1974) used the equivalent width of interstellar Ca ii K,
independent from the reddening diagnostics, to derive a distance of $2.2\pm
0.4$ kpc.

Based on the spectral type (Hiltner et al., 1972) and duration of the eclipse
(Forman et al., 1973), Conti (1978) derived a stellar radius of 35 $R_{\odot}$
and an absolute magnitude $M_{V}=-7.0$, which then led to an estimated
distance of 2.2 kpc (no uncertainty given) based on the photometry by Jones &
Liller (1973).

Sadakane et al. (1985) used ultraviolet spectra obtained with IUE to determine
an effective temperature of $25\,000\pm 1000$ K. They combined this
temperature with a radius of 31 $R_{\odot}$ (Rappaport & Joss, 1983) to arrive
at an absolute visual magnitude $M_{V}=-6.67$ and a derived distance of
$1.9\pm 0.2$ kpc.

From a spectroscopic analysis using NLTE models (see section 3.2.3),
Vanbeveren et al. (1993) derived a distance of 1.8–2 kpc.

Kaper et al. (1997) noted the possible connection to the Vela OB1 association
(see below) and introduced the distance of 1.82 kpc to the Vela OB1
(Humphreys, 1978) as a distance value. This was then used by several
subsequent X-ray studies without considering the implied uncertainties.

Coleiro & Chaty (2013) used a spectral energy distribution (SED) fitting
procedure in order to derive the distance and absorption of a sample of HMXBs.
For Vela X-1, they derived a value of 2.2$\pm$0.2 kpc.

This approach was further refined by Giménez-García et al. (2016). They
combined data from IR (2MASS) to UV (IUE) for the SED with a detailed analysis
of spectral lines using NLTE models created with the PoWR code (see section
3.2.3) in order to derive the extinction, luminosity, and distance to the
source, arriving at a distance of 2.0$\pm$0.2 kpc.

In the huge catalog of stellar distances based on Gaia Data Release 2 (Gaia
DR2, Gaia Collaboration, 2018) parallaxes, published by Bailer-Jones et al.
(2018), the distance to Vela X-1 is derived as 2.42$-$0.16$+$0.19 kpc from a
Bayesian analysis, with a basic parallax of 0.38$\pm$0.03 $\mu$as. Figure 13
shows that this distance is farther away than most previous estimates and is
formally consistent with only very few estimates.

The recently published third Gaia data release (EDR3, Gaia Collaboration et
al. 2020) gives us the opportunity to obtain an improved distance to Vela X-1.
Gaia EDR3 lists a parallax $\varpi$ = 496.2$\pm$15.2 $\mu$as. To derive a
distance, we first applied a parallax zero-point of $-$13.5 $\mu$as using the
Lindegren et al. (2020a) recipe for the relevant magnitude, color, and
ecliptic latitude. For this specific star, this zero-point does not differ
much from the median quasar value of $-$17 $\mu$as, but we note that a star as
bright as Vela X-1 ($G_{\rm EDR3}=6.7351$) is at the limit of the analysis in
subsection 4.4 of Lindegren et al. (2020a). Next, we considered transforming
from the internal parallax uncertainty of 15.2 $\mu$as to an external
uncertainty. Following Lindegren et al. (2018), astrometric uncertainties
should be transformed from their internal ($\sigma_{i}$) to their external
($\sigma_{e}$) uncertainties using

$\sigma_{e}=\sqrt{k^{2}\sigma_{i}^{2}+\sigma_{s}^{2}},$ (14)

where $k$ is a value to be determined (and likely a function of at least
magnitude), and $\sigma_{s}$ is the square root of the limit when the angular
covariance reaches zero. As of the time of this writing, the value of $k$ has
only been estimated for stars significantly fainter than Vela X-1 (Fabricius
et al., 2020; Maíz Apellániz et al., 2021). We used a value of $k=1$, which
may be a slight underestimation, but we compensate for this by using a very
conservative value of $26$ $\mu$as for $\sigma_{s}$, which is derived from
(faint) quasars by Lindegren et al. (2020b) and is likely to be an
overestimation because the checkered pattern seen in the EDR3 data for the LMC
has a smaller amplitude (Lindegren et al., 2020b). In this way, we arrive at a
corrected parallax of $\varpi$ = 509.7$\pm$30.1 $\mu$as, which may be revised
in the near future when the systematic uncertainties in Gaia EDR3 are better
characterized. Using this parallax as an input and the prior for runaway OB
stars from Maíz Apellániz (2001, 2005) with the updated parameters from Maíz
Apellániz et al. (2008), we arrive at a distance of $1.99^{+0.13}_{-0.11}$ kpc
for Vela X-1. We note from the analysis of Pantaleoni González et al. (2021)
with Gaia DR2 data (which should also apply to the EDR3 case) that when a
reasonable prior is used for the stellar spatial distribution in the Galaxy,
distances for stars with precise Gaia parallaxes are robust, that is, they are
nearly independent of the specific prior, especially for objects within a few
kiloparsec.

Figure 14: Comparison of the proposed tracks for the Vela X-1 system relative
to the Vela OB1 association following Figure 7 of van Rensbergen et al.
(1996), Figure 3 of Kaper et al. (1997), and Figure 7 of Gvaramadze et al.
(2018). The stars shown as belonging to the association (see Section 5.1.2)
are listed in Table 12. The numbers along the tracks indicate where the system
would have been one, two, or three million years ago, based on the assumptions
in these publications. The direction of proper motion, as listed in the
Hipparcos Input Catalog (HIC) (Turon et al., 1992) and as obtained by Sahu
(1992), is shown by gray vectors. The inset is adapted from Figure 6 of Maíz
Apellániz et al. (2018) and shows a $14^{\prime}\times 14^{\prime}$ WISE RGB
mosaic. Green and yellow arrows indicate the raw and corrected proper motions,
respectively (see the reference for details). The vectors shown in the inset
are not to scale with those of the main drawing.

####  5.1.2 Possible origin in the Vela OB1 association

In a wider study of six stars and four OB associations, van Rensbergen et al.
(1996) found that Vela X-1 presumably originated in the Vela OB1
association888Wright (2020) lists Vela OB1 under “OB associations for which
very little is known”. (Humphreys, 1978; Melnik & Dambis, 2020) of massive
stars. They found that HD 77581 was a runaway OB star, probably expelled from
the Vela OB1 association by the natal kick that occurred during the supernova
explosion associated with the formation of the neutron star.

The connection with Vela OB1 was again reported by Kaper et al. (1997), who
discovered a symmetric wind bow shock about 0.9 arcmin north of HD 77581 in a
narrow-band H$\alpha$ image. This bow shock is created by the balance between
the ram pressure of the stellar wind and the ambient interstellar medium, not
to be confused with the much smaller bow shock that formed between the blue
supergiant and the neutron star by the motion of the latter in the stellar
wind of the former (Sect. 2.3). Based on the symmetry axis of the bow shock,
they proposed a different direction of proper motion than was assumed by van
Rensbergen et al. (1996). They also noted that the distance of 1.82 kpc to
Vela OB1 given by Humphreys (1978) agreed with the distance to Vela X-1
reported by Sadakane et al. (1985) and the X-ray luminosity, thus popularizing
the use of this distance.

In a follow-up study, Huthoff & Kaper (2002) searched for bow shocks around
other HMXBs or single OB runaway stars, but found no other example in the
first group and only a minority in the second. They concluded that the
formation of a bow show strongly depends on the temperature and density of the
ambient medium and that apparently the majority of OB runaway stars is moving
through hot bubbles in the Galactic plane.

In their search for Galactic runaway stars based on the Gaia Data Release 1,
Maíz Apellániz et al. (2018) found the bow shock around HD 77581 (GP Vel in
that publication) also in images obtained by WISE999Wide-field Infrared Survey
Explorer. and demonstrated that the proper motion taken from the Tycho-Gaia
Astrometric Solution (TGAS, Michalik et al., 2015) aligned very well with the
symmetry axis of the shock, provided it was corrected for the mean proper
motion of OB stars in that Galactic direction (see inset in Figure 14).

Motivated by the results of Kaper et al. (1997), Gvaramadze et al. (2018)
searched for an H$\alpha$ counterpart in the Southern H-alpha Sky Survey Atlas
(SHASSA, Gaustad et al., 2001). They found extended filamentary structures
downstream of the bow shock within a wider area and presented the results of
optical spectroscopy of the bow shock, comparing the observational data with
3D magnetohydrodynamic simulations. Recalculating the space velocity of the
system (Tab. 2), Gvaramadze et al. (2018) confirmed the possible association
with the Vela OB1 association, but warned that this connection might be
spurious if the original binary system was ejected from its parent star
cluster before the supernova explosion. Even the lower estimate for the actual
space velocity of Vela X-1 places the system in the top few percentiles of the
velocity distribution of runaway stars derived from recent numerical
calculations of massive binary evolution (Renzo et al., 2019).

Table 2: Space velocity estimates of the Vela X-1 system. Reference | space velocity
---|---
Kaper et al. (1997) / Turon et al. (1992) |  $89\pm 40$ km s-1
Kaper et al. (1997) / Sahu (1992) |  $98\pm 26$ km s-1
Gvaramadze et al. (2018) |  $54.25\pm 0.55$ km s-1

As found in the previous section, the distance to Vela X-1 based on Gaia data
is not fully consistent with the commonly assumed mean distance to the Vela
OB1 association. Therefore we also used Gaia EDR3 astrometric data to revisit
Vela OB1 and obtain our own estimate of its distance.

Recently, Melnik & Dambis (2020) cross-matched the catalog of Blaha &
Humphreys (1989) with Gaia DR2 using a matching radius of 3 arcsec and a
magnitude tolerance of 3 mag. They identified 46 Vela OB1 sources with Gaia
DR2 counterparts. We cross-matched these 46 sources with Gaia EDR3 and
obtained a counterpart for all of them. However, we rejected 7 sources with
poor astrometric solution (renormalized unit weight error,
RUWE,101010https://www.cosmos.esa.int/web/gaia/dr2-known-issues$>$1.4).

We assumed that all Vela OB1 members have similar distances (parallaxes) and
kinematic (proper motions). We therefore applied an iterative 3$\sigma$
clipping in parallax, proper motion in RA, and proper motion in DEC. This
reduced the sample to 30 Gaia EDR3 sources with compatible parallaxes and
proper motions that are considered to belong to Vela OB1. They are shown in
Fig. 14 and listed in Table 12. This table contains the Gaia EDR3 ID, the
coordinates of the source, and the parallax and proper motion used in this
work.

Using this reduced sample, we estimated a mean parallax for Vela OB1 of
0.53$\pm$0.07 mas. This corresponds to a mean distance of 1.88$\pm$0.24 kpc
and is based on the inverse of the parallax and error propagation alone.

The distance to Vela OB1 is therefore compatible with a possible origin of
Vela X-1 in this stellar association. However, the estimated distances to Vela
OB1 have large uncertainties that for the moment prevent us from reaching a
firm conclusion about the origin of Vela X-1. It is expected that future Gaia
data releases will allow us a better determination of the distance to Vela OB1
and Vela X-1.

Table 3: Orbital parameters of the Vela X-1/HD 77581 system as derived by a selection of authors over the years. We also indicate where subsequent publications directly reused previous results. Different definitions of the time of phase zero have been used over time: Mid-eclipse time ($T_{\mathrm{ecl}}$) or time of mean longitude equal to $\pi/2$ ($T_{\pi/2}\equiv T_{\mathrm{90}}$) are the most common, but some authors also use periastron time ($T_{\mathrm{peri}}$ or $\tau$) or specific zero times (Hutchings, 1974). Reference | Time | Phase zero | Orbital period | semimajor axis | eccentricity | longitude of
---|---|---|---|---|---|---
| basis | [MJD] |  $P_{\mathrm{orb}}$ [d] |  $a_{\mathrm{X}}\sin i$ [lt-sec] | $e$ |  periastron $\omega$
Forman et al. (1973) | $T_{\mathrm{ecl}}$ |  41446.04$\pm$ | 0.07 |  8.95$\pm$ | 0.02 | – |  0.14$\pm$ | 0.05 | –
Hutchings (1974) (H74) | $T_{\mathrm{phot}}$ |  41593.40$\pm$ | 0.56 |  8.966$\pm$ | 0.001 | – | – | –
Rappaport et al. (1976) | $T_{\mathrm{ecl}}$ |  42611.501$\pm$ | 0.12 | (H74) |  111.4$\pm$ | 3.3 |  0.126$\pm$ | 0.041 |  146$\pm$ | 23
Ögelman et al. (1977) (Ö77) | $T_{\mathrm{ecl}}$ |  42446.03$\pm$ | 0.06 |  8.9643$\pm$ | 0.0005 |  113.0$\pm$ | 4.0 | – |  162$\pm$ | 25
van Paradijs et al. (1977b) | – | – |  8.9681$\pm$ | 0.0016 | – |  0.136$\pm$ | 0.046 | –
Watson & Griffiths (1977) | $T_{\mathrm{ecl}}$ |  42620.30$\pm$ | 0.05 |  8.9640$\pm$ | 0.0005 | – | – | –
Rappaport et al. (1980) | $T_{\mathrm{peri}}$ |  42822.90$\pm$ | 0.13 |  8.9649$\pm$ | 0.0002 |  113.0$\pm$ | 0.8 |  0.092$\pm$ | 0.005 |  154$\pm$ | 5
Nagase et al. (1983) | $T_{\mathrm{ecl}}$ |  41446.02$\pm$ | 0.05 |  8.9640$\pm$ | 0.0006 | – | – | –
Nagase et al. (1984b) | $T_{\mathrm{peri}}$ |  45329.05$\pm$ | 0.06 |  8.9642$\pm$ | 0.0006 |  114.10$\pm$ | 0.50 |  0.080$\pm$ | 0.006 |  157.3$\pm$ | 2.1
van der Klis & Bonnet-Bidaud (1984) | $T_{\pi/2}$ |  42727.750$\pm$ | 0.024 | (Ö77) | (Ö77) |  0.116$\pm$ | 0.022 |  166$\pm$ | 7
– ditto – | $T_{\pi/2}$ |  42966.628$\pm$ | 0.019 | (Ö77) |  113.4$\pm$ | 1.7 |  0.085$\pm$ | 0.022 |  121$\pm$ | 16
– ditto – | $T_{\pi/2}$ |  44170.937$\pm$ | 0.021 | (Ö77) |  113.20$\pm$ | 1.4 |  0.102$\pm$ | 0.024 |  177$\pm$ | 15
Boynton et al. (1986) | $T_{\pi/2}$ |  43821.3604$\pm$ | 0.0056 |  8.96443$\pm$ | 0.00022 |  112.70$\pm$ | 0.47 |  0.0881$\pm$ | 0.0036 |  152.8$\pm$ | 2.2
Deeter et al. (1987) | $T_{\pi/2}$ |  43955.8261$\pm$ | 0.0048 |  8.964353$\pm$ | 0.000063 |  112.66$\pm$ | 0.41 |  0.0896$\pm$ | 0.0031 |  152.8$\pm$ | 1.8
Deeter et al. (1987) | $T_{\pi/2}$ |  44278.5466$\pm$ | 0.0037 |  8.964416$\pm$ | 0.000049 |  112.98$\pm$ | 0.35 |  0.0885$\pm$ | 0.0025 |  150.6$\pm$ | 1.8
Bildsten et al. (1997) (B97) | $T_{\pi/2}$ |  48895.2186$\pm$ | 0.0012 |  8.964368$\pm$ | 0.00004 |  113.89$\pm$ | 0.13 |  0.0898$\pm$ | 0.0012 |  152.59$\pm$ | 0.92
Kreykenbohm et al. (2008) | $T_{\pi/2}$ |  52974.001$\pm$ | 0.012 |  8.964357$\pm$ | 0.000029 | (B97) | (B97) | (B97)
– ditto – | $T_{\mathrm{ecl}}$ |  52974.227$\pm$ | 0.007 | same as above |  |  |  |  |  |
Falanga et al. (2015) | $T_{\pi/2}$ |  42611.349$\pm$ | 0.013 |  8.964427$\pm$ | 0.000012 | (B97) | (B97) | (B97)
– ditto – | $T_{\mathrm{ecl}}$ |  42611.1693$\pm$ | 0.013 |  8.9644061$\pm$ | 0.0000064 |  |  |  |  |  |

###  5.2 Orbital parameters and ephemerides

Since the first clearly determined orbital period (Forman et al., 1973), the
orbital parameters have been further refined in the following years and
decades by various different groups based on eclipse timing, the analysis of
Doppler shifts in the X-ray pulsations, or RVs from optical spectroscopy, for
example. The values have essentially converged for the main parameters. Table
3 presents an overview of this evolution and gives the references. In the more
recent literature, the orbital parameters are frequently taken from Bildsten
et al. (1997) or from the updated values from Kreykenbohm et al. (2008), who
largely built on the previous work. Falanga et al. (2015) have undertaken a
systematic study of eclipsing HMXBs and derived new orbital parameters for
Vela X-1, which slightly contradict Bildsten et al. (1997) or Kreykenbohm et
al. (2008); see Table 3 and Figure 16.

Although the remaining differences are small, the orbital phases calculated
using ephemerides published in the past three decades (Deeter et al., 1987;
Bildsten et al., 1997; Kreykenbohm et al., 2008; Falanga et al., 2015) agree
very well, at a level of a few times $10^{-3}$ in orbital phase all the way
back to the detection of Vela X-1. The systematic phase shift induced by
choosing either $T_{\mathrm{ecl}}$ or $T_{\pi/2}$ is significantly larger,
about 0.025$\pm$0.01 (Kreykenbohm et al., 2008; Falanga et al., 2015). Figure
15 illustrates these facts for a recent date. Care should be taken in any
discussion of orbital variations that it is clearly indicated which ephemeris
and zero time is used.

Figure 15: Predicted orbital phase for a given date, demonstrating the small
differences between different ephemerides. In contrast, the difference in
phase between mid-eclipse time and $T_{\pi/2}$ as zero-point, as shown for
Kreykenbohm et al. (2008) and Falanga et al. (2015), is significant.

The orbit of Vela X-1 is very stable and shows no indication of orbital decay.
As the most recent study, Falanga et al. (2015), combining results from back
to Forman et al. (1973) with data up to 2011 derived a period decay consistent
with zero ($\dot{P}_{\mathrm{orb}}/P_{\mathrm{orb}}=0.1(3)\times
10^{-6}\,\mathrm{yr}^{-1}$). They reported indications for a slight apsidal
advance ($\dot{\omega}=0.41(27)\mathrm{yr}^{-1}$), in contrast with previous
studies (Deeter et al., 1987, and references therein).

Figure 16: Evolution of orbital period determinations since Hutchings (1974),
converging on very similar values by different groups.

###  5.3 Masses, radii, and inclination

Table 4: Estimated masses of the two binary components of the Vela X-1/HD 77581 system as given by various authors, together with the assumed inclination $i$. Where provided, the estimated radius of the optical star is also given. Errors or uncertainty ranges are given as in the original publications. See text for the multiple solutions provided by van Kerkwijk et al. (1995) or Stickland et al. (1997). The latter underestimated $K_{\textrm{opt}}$ according to Barziv et al. (2001) due to a data reduction error. Quaintrell et al. (2003) give a range of values, with the two extreme cases tabulated. Reference | $i$ |  $M_{\mathrm{NS}}$ [$M_{\sun}$] |  $M_{\mathrm{opt}}$ [$M_{\sun}$] |  $R_{\mathrm{opt}}$ [$R_{\sun}$] |  $K_{\mathrm{opt}}$ [km s-1]
---|---|---|---|---|---
Hutchings (1974) |  | 90 | 1.5–2.7 | – | – | $20.2\pm 1.1$
– ditto – |  | 70 | 2.1–2.9 | – | – | $20.2\pm 1.1$
Mikkelsen & Wallerstein (1974) |  | 90 | 1.0–3.0 | 10–30 | $\sim 30$ | –
– ditto – (plausible $M_{\mathrm{Opt}}$) |  | 90 | 1.5–3.0 | 20–25 | $\sim 30$ | –
Zuiderwijk et al. (1974) |  | 90 | $\geq 2.5\pm 0.3$ | $\geq 30\pm 5$ | – | $26\pm 0.7$
van Paradijs et al. (1976) |  | 74–90 | $1.61\pm 0.27$ | $21\pm 2.6$ | – | $20\pm 1$
Lamers et al. (1976) |  | – | 2.0 | 22 | 27 | –
Rappaport et al. (1976) |  |  $\geq$70 | $1.45\pm 0.16/(\sin i)^{3}$ | $21\pm 0.9/(\sin i)^{3}$ | – | –
Vidal (1976) |  | 90 | 1.56–2.1 | 18.1–23.0 | – | –
– ditto – |  | 80 | 1.64–2.2 | 19.0–23.7 | – | –
– ditto – |  | 70 | 1.9–2.55 | 20.7–27.7 | – | –
van Paradijs et al. (1977b) |  | $>74$ | $1.67\pm 0.12$ | $20.5\pm 0.9/(\sin i)^{3}$ | – | $21.75\pm 1.15$
Conti (1978) |  | 74 | 1.6 | 24 | 35 | –
Vanbeveren et al. (1993) |  | – | – | 21.6–26.5 | 31–35 | –
Joss & Rappaport (1984) | (JR84) | $>73$ | $1.85^{+0.35}_{-0.30}$ | $23.0^{+3.5}_{-1.5}$ | $31^{+4}_{-3}$ | –
Bulik et al. (1995) |  | – | $1.34\pm 0.16$ | – | – | –
van Kerkwijk et al. (1995) | #1 | $>74$ | $1.88^{+0.69}_{-0.47}$ | $23.5^{+2.2}_{-1.5}$ | $30.0^{+1.9}_{-1.9}$ | $17.0-29.7$
– ditto – | #2 | $>73$ | $1.75^{+0.34}_{-0.33}$ | $23.3^{+2.5}_{-1.1}$ | $30.2^{+1.7}_{-2.2}$ | $20.8\pm 1.7$
– ditto – | #3 | $>75$ | $2.08^{+0.46}_{-0.43}$ | $23.6^{+1.9}_{-2.2}$ | $29.9^{+1.4}_{-1.8}$ | $24.6\pm 2.3$
Stickland et al. (1997) | (a) | $>73$ | $1.39\pm 0.14/(\sin i)^{3}$ | $21.41\pm 0.22/(\sin i)^{3}$ | – | $17.8\pm 1.6$
– ditto – | (b) | $>73$ | $1.33\pm 0.13/(\sin i)^{3}$ | $21.42\pm 0.22/(\sin i)^{3}$ | – | $17.8\pm 1.6$
– ditto – | (c) | $>73$ | $1.45\pm 0.14/(\sin i)^{3}$ | $21.42\pm 0.22/(\sin i)^{3}$ | – | $17.8\pm 1.6$
Barziv et al. (2001) |  | $>73$ | $1.86\pm 0.32$ | $23.8^{+2.4}_{-1.0}$ | $30.4^{+1.6}_{-2.1}$ | $21.7\pm 1.6$
Quaintrell et al. (2003) |  | $70.1\pm 2.6$ | $2.26\pm 0.32$ | $27.9\pm 1.3$ | $32.1\pm 0.6$ | $21.2\pm 0.7$
– ditto – |  | 90 | $1.88\pm 0.13$ | $23.1\pm 0.2$ | $26.8\pm 0.9$ | $22.6\pm 1.5$
Rawls et al. (2011) (analytic) |  | $83.6\pm 3.1$ | $1.788\pm 0.157$ | – | – | –
– ditto – (numerical) |  | $78.8\pm 1.2$ | $1.770\pm 0.083$ | $24.00\pm 0.37$ | $31.82\pm 0.37$ | –
Krtička et al. (2015)† |  | – | – | 23.5 | 30 | –
Falanga et al. (2015) |  | $72.8\pm 0.4$ | $2.12\pm 0.16$ | $26\pm 1$ | $29\pm 1$ | –
Giménez-García et al. (2016) |  | – | – | $21.5\pm 4.0$ | (JR84) | –
† Refer to van Kerkwijk et al. (1995)

Figure 17: Upper panel: Selected estimates of the mass of the neutron star in
Vela X-1 taken from the literature. The estimates are for specific assumed
inclinations (see the original papers for how these where chosen) or plotted
for a range of values with $i=90^{\circ}$ as the lowest mass point and then
$5^{\circ}$ steps to $i=70^{\circ}$, to show the $\sin(i)$ dependence. See
text for the multiple solutions provided by van Kerkwijk et al. (1995). The
shaded areas on the right indicate possible maximum masses for neutron stars
based on multi-messenger observations of GW170817 (see Margalit & Metzger,
2017, and references therein). To set them into perspective, the lower panel
shows the distribution of derived masses of other neutron stars in the
literature taken from Alsing et al. (2018) for comparison.

Most commonly, the masses of the two components of the binary have been
determined from fits to RV curves (see Sect. 3.2.2). In some cases, the mass
and radius of HD 77581 have been estimated by other means. Lamers et al.
(1976), for example, used half the orbital separation as the stellar radius
estimate because this agreed with the mean value for these supergiants, and
matched the temperature and evolutionary tracks to determine the mass of the
donor star.

Early RV observations (Hiltner et al., 1972) were interpreted as indications
of a possible black hole in the binary system (Wickramasinghe et al., 1974),
but with more measurements, this interpretation was soon changed to a high-
mass neutron star system instead, with possible masses between 1.5 and 2.9
$M_{\sun}$ (Hutchings, 1974). In the following years, various other authors
have undertaken to determine the masses of the two partners, sometimes using
the same basic data sets, but arriving at different results.

As examples of how the results of the analysis may depend on assumptions, van
Kerkwijk et al. (1995) derived three different sets of solutions for different
ranges of the RV amplitude $K_{\mathrm{opt}}$ derived under different
assumptions: The first set is for a conservative range, set 2 is for the best
fit using all their data, and set 3 for $K_{\mathrm{opt}}$ derived excluding
orbital phase 0.65–0.85, in which they found the strongest indications for
systematic effects. As inclination constraints for these three cases, they
used $i>74$, $i>73$ and $i>75$ degrees, respectively. Later, Stickland et al.
(1997) used the same data set as van Kerkwijk et al. (1995), but a different
fitting procedure. They derived systematically lower estimates for the neutron
star mass and reported three different possible solutions for (a) an
unconstrained orbital fit, (b) using the orbital solution from Deeter et al.
(1987), and (c) with the spectra of the comparison star taken with a small
aperture.

An interesting complementary approach was used by Bulik et al. (1995), who fit
pulse-phase resolved X-ray spectra with a model of X-ray emission from polar
caps in an inhomogeneous highly magnetized neutron star atmosphere. They
included general relativistic effects to derive constraints on the neutron
star mass and radius. Their results depend on assumptions about the emission
region, however, and they reported systematic differences for fits made to
different energy ranges.

Table 4 compiles different mass estimates, and Figure 17 visualizes a
selection of results for the neutron star mass in comparison with the mass
distribution of other neutron stars with existing estimates.

In summary, while there are clear indications for a relatively high-mass
neutron star in Vela X-1, with a mass of about 1.8 $M_{\sun}$ or even higher,
the systematic uncertainties in the derivation of the mass still allow a
rather wide range of masses. The only moderately constrained inclination is a
significant factor that by itself is responsible for a difference of $\sim$15%
between the lowest and highest mass of a given study, but there is larger
scatter from the different approaches. The radius of the neutron star has not
been derived from observations. It will depend on the actual mass and the
equation of state of the neutron star. Discussing the many solutions for the
equation of state and the corresponding mass-radius relations is beyond the
scope of this article (see, e.g., Lattimer & Prakash, 2007, for a review), but
based on the tracks shown in Maselli et al. (2020) for models compatible with
current constraints, we arrive at possible radii between 11 and 12.5 km for a
plausible range in neutron star masses in the Vela X-1 system. Similarly, the
mass of HD 77581 is effectively constrained to a range 20–30 $M_{\sun}$ ,
while its radius is mostly found to be $\sim 30R_{\sun}$.

Finally, indirect theoretical arguments and numerical approaches can also be
used in favor of certain values for the mass of the neutron star. For
instance, Manousakis et al. (2012) proposed a mass estimate based on the
geometry of the accretion flow as deduced from the absorption column density
profile: All things being equal, the higher the mass of the neutron star, the
more distorted the wind they obtained in their hydrodynamical simulations,
leading to different column density profiles. Although they considered another
HMXB, IGR J17252$-$3616, as an example case, their method would be applicable
to Vela X-1 as well. We also note that for a mass ratio lower than $\sim$12,
the absence of RLOF (see Section 2.2) sets a more stringent lower limit on the
orbital inclination than the eclipse duration.

Table 5: Parameters of HD 77581 derived by different authors. See Sect. 3.2.3 for notes on the different approaches used. For Giménez-García et al. (2016) and Sander et al. (2018), the temperature we report here as $T_{\star}$ corresponds to the parameter $T_{2/3}$ as defined in the PoWR code used by these publications. For mass and radius estimates, see Table 4. Reference | $\log L$ | $E_{B-V}$ | $T_{\star}$ | $v_{\mathrm{rot}}\sin i$ | $\log g_{\star}$ |
---|---|---|---|---|---|---
|  [$L_{\sun}$] | [mag] | [kK] |  [km s-1] | [cgs] | Comments
Wickramasinghe et al. (1974) | – | 0.75 | 20–28 | $\sim 90$ | 2.5–3.0 | RV, spectrum and reddening
Mikkelsen & Wallerstein (1974) (MW74) | 5.2 | 0.75 | 22 | 130 | 2.64–2.80 |
Lamers et al. (1976) | – | – | 26.3 | – | – |
Ammann & Mauder (1978) | – | – | – | $89\pm 5$ | – | HeI line fitting
Conti (1978) (C78) | 5.7 | 0.73 | 26 | (MW74) | – |
Dupree et al. (1980) | (C78) | 0.7 | (C78) | – | – |
Sadakane et al. (1985) | 5.53 | – | $25\pm 1$ | – | – |
Vanbeveren et al. (1993) | 5.5–5.7 | – | 26 | 130 | 2.75–2.76 |
Zuiderwijk (1995) | – | – | – | $116\pm 6$ | – | HeI line profiles
Howarth et al. (1997) | – | – | – | 114 | – | cross-correlating IUE spectra
Fraser et al. (2010) (F10) | – | – | 26.5 | 56 | 2.90 |  Arrive at mass of 40$M_{\odot}$
Falanga et al. (2015) | – | – | – | $130.3\pm 0.2$ | – | From system geometry
Krtička et al. (2015) | 5.63 | – | 27 | – | – |
Giménez-García et al. (2016) | $5.5\pm 0.1$ | $0.77\pm 0.05$ | $24.4\pm 1.0$ | (F10) | $2.86\pm 0.10$ |  Uses $d=2.0\pm 0.2$ kpc
Sander et al. (2018), no X-rays | 5.485 | – | 23.5 | (F10) | 2.84 | same distance as above
– ditto –, moderate X-ray illum. | – ditto – | – | 23.5 | – ditto – | – |
– ditto –, strong X-ray illum. | – ditto – | – | 23.8 | – ditto – | – |
this work (Section 5.5) | 5.8–6.2 | $0.689\pm 0.018$ | $33.7\pm 5.2$ | – | $3.28\pm 0.19$ |  Uses $d=1.99$ kpc

###  5.4 Revising the spectral classification of HD 77581

SIMBAD lists different spectral classifications of HD 77581. They range in
spectral subtype from B0 to B0.5 and in luminosity class from Ib to Ia. This
means that the different authors agree that it is an early-type B supergiant,
but disagree on the precise spectral classification. This disagreement is
quite common for OB stars because many objects were classified in the 1950s to
1970s using photographic plates, and they have not been revisited with modern
digital data since then. For this reason, the Galactic O-Star Spectroscopic
Survey (GOSSS, Maíz Apellániz et al. 2011) was started a decade ago. It was
originally restricted to the modern spectral classification of O stars and was
later extended to B stars. A GOSSS spectrogram of HD 77581 was presented in
Maíz Apellániz et al. (2018) and given the classification of B0.5 Ia. Since
that paper, the GOSSS team (of which the PI is a coauthor here) has revisited
the spectral classification criteria for B supergiants and is currently
working on a new publication on this issue. The two main criteria for the
horizontal classification of early-B stars are the ratio of Si iii to Si iv
lines and the progressive disappearance of He ii lines in the spectrum with
evolution from B0 to B1. Originally, He ii lines were only seen in O stars,
but with improvements in S/N and the advent of digital data, the current
criteria establish that for supergiants, He ii $\lambda$4542 disappears at
B0.5 (and is just visible at B0.2) and that He ii $\lambda$4686 disappears at
B1 (and is just visible at B0.7). HD 77581 shows a weak He ii $\lambda$4542
absorption, and the ratio of the Si iii to Si iv lines is intermediate between
those at B0 and B0.5. On the other hand, the ratio of the Si iii to He i lines
and the H$\beta$ to H$\delta$ profiles indicates a high luminosity class of Ia
(but not Ia+ because H$\beta$ does not show a P-Cygni profile). Therefore, we
revise the GOSSS spectral classification from B0.5 Ia to B0.2 Ia.

###  5.5 New determination of stellar parameters

Figure 18: CHORIZOS mode SED. Blue error bars are used for the input
photometry that was used for the fit (Gaia DR2 $G+G_{\rm BP}+G_{\rm RP}$,
Johnson $U+B+V$, and 2MASS $J+H+K_{\rm s}$) and red error bars for the
discarded input photometry (WISE $W1+W2+W3+W4$). In each case, the vertical
error bar is the photometric uncertainty, and the horizontal error bar is a
representation of the wavelength coverage of each filter. Green stars are used
for the model photometry.

Based on newer data and our distance estimate (see Section 5.1), we used the
latest version of CHORIZOS (Maíz Apellániz, 2004) to again derive some of the
parameters of HD 77581. CHORIZOS is a Bayesian photometric code that computes
the likelihood that a set of photometric data is compatible with a family of
SEDs.

Photometric data. We collected the magnitudes from Gaia DR2 $G+G_{\rm
BP}+G_{\rm RP}$ (Evans et al. 2018 with the corrections and calibration from
Maíz Apellániz & Weiler 2018), Johnson $U+B+V$ (Mermilliod et al. 1997 with
the calibration from Maíz Apellániz 2006, 2007), 2MASS $J+H+K_{\rm s}$
(Skrutskie et al. 2006 with the calibration from Maíz Apellániz & Pantaleoni
González 2018), and WISE $W1+W2+W3+W4$ (Cutri & et al., 2013). As the object
has a significant IR excess, we eliminated the four photometric points with
the longest effective wavelengths from the runs (the four WISE bands), but
included them in the analysis to evaluate the excesses.

Models. We used the $T_{\rm eff}$-luminosity class grid with solar metallicity
from Maíz Apellániz (2013b). The temperature of interest lies in the region
for which the grid uses the TLUSTY models (Lanz & Hubeny, 2007) in the optical
and the Munari models (Munari et al., 2005) in the NIR. The reason for using
these combined models is that the TLUSTY SEDs for B stars yield incorrect
photometric values by up to several hundredths of a magnitude in the NIR. The
Munari models do not have this problem.

Parameters. For our run we fixed the value of the logarithmic distance $\log
d$ to 3.299, and we left four free parameters: $T_{\rm eff}$, luminosity
class, $E(4405-5495)$ (amount of extinction), and $R_{5495}$ (type of
extinction). We used the family of extinction laws of Maíz Apellániz et al.
(2014), see Maíz Apellániz (2004, 2013a) and Maíz Apellániz & Barbá (2018) for
an explanation of why monochromatic quantities need to be used instead of
band-integrated quantities to characterize extinction. Because we fit nine
photometric points, we had five degrees of freedom.

Results. The results of the CHORIZOS run are given in Table 6, and the mode
SED is shown in Fig. 18. The reduced $\chi^{2}$ of the mode SED is 0.75,
indicating an excellent fit of the model to the data. This is shown in Fig.
18, where all green stars up to $K_{\rm s}$ are within the blue error bars or
very close to them. The first five lines of Table 6 give the results for the
model parameters (mean, uncertainty, and mode), and the last nine lines list
the results for the derived parameters (mean and uncertainty). We discuss the
model parameters first.

Table 6: Detailed results of the CHORIZOS analysis. Quantity | Units | Mean | Unc. | Mode
---|---|---|---|---
Model parameters
$T_{\rm eff}$ | kK | 33. | 7 | 5. | 2 | 30. | 9
Lum. class | — | 0. | 65 | 0. | 31 | 0. | 84
$R_{5495}$ | — | 3. | 505 | 0. | 071 | 3. | 519
$E(4405-5495)$ | mag | 0. | 689 | 0. | 018 | 0. | 681
$\log d$ | log pc | 3. | 299 | fixed | 3. | 299
Derived parameters from evolutionary tracks
Luminosity |  $10^{6}$ solar | 1. | 05 | 0. | 46 |
$\log g$ | log cgs | 3. | 28 | 0. | 19 |
ZAMS mass | solar | 71. | 2 | 22. | 6 |
age | Myr | 2. | 75 | 0. | 78 |
$A_{5495}$ | mag | 2. | 367 | 0. | 064 |
$A_{V_{J}}$ | mag | 2. | 445 | 0. | 063 |
$E(B-V)$ | mag | 0. | 691 | 0. | 018 |
$V_{J,0}$ | mag | 4. | 478 | 0. | 057 |
$M_{V_{J,0}}$ | mag |  $-$7. | 017 | 0. | 057 |

$T_{\rm eff}$ is poorly constrained because the main discerning criterion is
the Balmer jump, which is determined almost exclusively by the Johnson $U-B$
(Maíz Apellániz et al., 2014). The value is above those obtained from
spectroscopy, but within 2$\sigma$. The luminosity class is the photometric
equivalent to the spectral luminosity class; it varies from 0.0 (hypergiants)
to 5.5 (ZAMS). The value found by CHORIZOS (which is dependent on the fixed
distance) indicates it is a bright supergiant, which is consistent with the Ia
spectral luminosity classification. The extinction parameters indicate that
the value of $R_{5495}$ is above the canonical value of 3.1, which is
relatively frequently the case for sightlines with $E(4405-5495)$ (Maíz
Apellániz & Barbá, 2018), and that $E(4405-5495)$ [or $E(B-V)$] is slightly
lower than previous values. We note, however, that the higher value of
$R_{5495}$ leads to only a small change in $A_{5495}$ (or in $A_{V_{J}}$)
compared to a measurement of $E(4405-5495)$ of 0.75 assuming a value of
$R_{5495}$ of 3.1, for example.

Figure 19: Comparison of different derived wind velocity profiles from the
literature (Tab. 7). See text for more details of the models. These are mostly
the solutions for the wind _without_ accounting for the ionization of the
neutron star, except in the case of two of the models by Sander et al. (2018)
and, at the bottom, the velocity profile calculated by Watanabe et al. (2006)
along the line connecting donor and neutron star. The range of possible
distances covered by the neutron star on its eccentric orbit (Bildsten et al.,
1997) and the orbital velocity range, mainly driven by the uncertainty in
inclination, are indicated as well. In the upper left corner we indicate
approximately how wind (shown as horizontal lines) and orbital velocity
(vertical) would combine at the mean neutron star distance from HD 77581 for
three of the reported velocity profiles. These differences have a strong
effect on the expected accretion flow close to the neutron star (see Sect.
4.1).

The derived parameters in Table 6 were obtained with the help of the Geneva
(Lejeune & Schaerer, 2001) evolutionary tracks used by the Maíz Apellániz
(2013b) grid and assuming the value of $\log d$ given there. The luminosity we
found is close to one million in solar units with a relatively large
uncertainty caused by the relatively large uncertainty in $T_{\rm eff}$ (and,
hence, on the bolometric correction). The values for $\log g$, ZAMS mass, and
age depend in a complex way on distance, and if different distance were to be
used, they would change in a nontrivial way. In any case, taken at face value
and assuming the validity of the used evolutionary tracks, the CHORIZOS
results indicate that HD 77581 is a very massive star that was born less than
three million years ago and that initially would have been a very early type O
star or even a WNh star. Based on the existence of the neutron star companion
and on the mass transfer that must have taken place in the past, the previous
initial state and evolution of the star are likely to have been different.
Moreover, the uncertainty in the absolute magnitude of the star in the Johnson
$V$ band is very low. This is so because (a) the extinction parameters are
very well determined, (b) for this quantity we are not affected by the
bolometric correction (as we are for the luminosity), and (c) we considered
the distance as a fixed quantity. For a different distance, we would need to
correct for the different value of $5\log d$.

Finally, we analyzed the behavior of the SED in the WISE bands in the MIR.
There are clear excesses there that we quantify as 0.135$\pm$0.066 mag ($W1$),
0.346$\pm$0.035 mag ($W2$), 0.346$\pm$0.027 mag ($W3$), and 2.007$\pm$0.035
mag ($W4$). We note that the Munari models do not include stellar winds.

###  5.6 Stellar wind parameters

Table 7: Parameters of the stellar wind in the Vela X-1 system as derived by different authors, see Sects. 3.2.3 and 3.2.4. When a value was taken from a previous publication, we mark this by a shorthand, e.g., D80 for Dupree et al. (1980). Numbers in italics in the table have been set as fixed in the respective publications. Krtička et al. (2012) and Sander et al. (2018) did not assume a basic $\beta$-law for wind acceleration, and the value given in this table corresponds to an approximate $\beta$-law fit to their results for comparison. For other parameters, see Table 4 and Table 5. Reference |  $v_{\mathrm{esc}}$ [km s-1] |  $v_{\infty}$ [km s-1] | $\beta$ |  $\dot{M}$ [$10^{-6}M_{\sun}\mathrm{yr}^{-1}$] | Comments
---|---|---|---|---|---
Lamers et al. (1976) (L76) | – | – | – | 0.5 | X-ray fluxes & wind accretion
Hutchings (1976) (H76) | – | – | – | 7 | spectroscopy, mass-loss model
Conti (1978) | 519 | 1600 | – | (H76) | Similar study to L76
Dupree et al. (1980) (D80) | 489 |  $\sim$1700 | 0.5 |  $\sim$1 |
McCray et al. (1984) (M84) | – | (D80) | 1 | 4 |
Sadakane et al. (1985) | – |  $\sim$1500 | – | (D80) | multiple components with different velocities
Sato et al. (1986) | – | 1600(a) | (D80) | $2.4\pm 1.2$ |
Prinja et al. (1990) | (590) | 1105 | – | – |  $v_{\mathrm{esc}}$ is mean for spectral type
Sako et al. (1999) | – | (D80) |  0.79$\pm$0.23 | $0.265^{+0.065}_{-0.050}\,{}^{(b)}$ | X-ray DEM modeling
van Loon et al. (2001) | – | $600\pm 100$ | – | – | Modeling of UV lines
Watanabe et al. (2006) | – | 1100(c) | 0.8 | 1.5 – 2.0 | X-ray line MC modeling
Krtička et al. (2012) |  | 750 |  ($\sim$0.8) | 1.5 | see reference for details
Falanga et al. (2015) | – | $640\pm 190$ | – |  $<$(1 – 5.3) | see reference for details
Manousakis & Walter (2015b) | – | (D80) | 0.55 | (M84) | Simulating hard X-ray light curve
Giménez-García et al. (2016) | $436\pm 65$ | $700^{+200}_{-100}$ | 1.0 | 0.4 – 1 |
Sander et al. (2018), no X-rays | – | 532 |  ($\sim$2.2) | 0.65 | see reference for details
– ditto –, moderate X-ray illum. | – | 584 |  ($\sim$2.2) | 0.66 |
– ditto –, strong X-ray illumination | – | 378 |  ($\sim$2.2) | 0.85 |
(c) Referring to Dupree et al. (1980).
${}^{(b})$ Argued to be the mass fraction in hot ionized gas, while most mass
would be contained in cool clumps.
(c) Referring to Prinja et al. (1990).

Over the years, many different approaches have been made to determine
essential parameters of HD 77581 related to its strong stellar wind. They are
collected in Tables 4, 5, and 7. As summarized in section 3.2.3, one common
approach based on the study of massive stars in general is through
(quantitative) spectroscopy. This has ranged from the study of specific lines
that were sometimes compared with samples of example stars (e.g.,
Wickramasinghe et al., 1974; Ammann & Mauder, 1978; Sadakane et al., 1985;
Prinja et al., 1990; Zuiderwijk, 1995; Howarth et al., 1997), to a detailed
SED analysis against grids of stellar model spectra created with increasingly
refined wind model codes (Vanbeveren et al., 1993; Fraser et al., 2010;
Giménez-García et al., 2016; Sander et al., 2018); see also section 4.1.

Other authors explicitly included the X-ray source in the derivation, for
example, by comparing the observed X-ray luminosity with that expected from a
model of stellar wind accretion (see section 2.2 (Lamers et al., 1976; Conti,
1978) ) or full hydrodynamic simulations (Manousakis & Walter, 2015a, b).
Others modeled the Hatchett-McCray effect on line features that vary with
orbital phase (Dupree et al., 1980; van Loon et al., 2001). Yet another
approach is based on the modeling of high-resolution X-ray spectra (section
3.1.4), assuming an ionized stellar wind of a certain structure, with some
parameters taken as given from previous work and others then derived from the
spectral modeling (Sako et al., 1999; Watanabe et al., 2006).

We also note that most of the efforts start from a smooth, nonclumped wind,
even though other physical effects that are taken into account then may break
the original full symmetry.

When we compare all these results, it is clear that the basic parameters of
the mass donor and its stellar wind are still only known within a significant
range of uncertainty. Care must be taken when astrophysical conclusions are
derived from a comparison of model results with any specific value.

For the rotational velocity, the derived values vary by a factor $>2$,
although in all cases, HD 77581 rotates significantly more slowly than the
orbital movement: by a factor of roughly one-third to two-thirds, based on the
values in Table 5. Giménez-García et al. (2016), who adopted the low value of
56 km s-1 from Fraser et al. (2010). We note that the high estimates of
rotational velocity are not consistent with optical lines that are observed to
be unblended. Additional broadening of lines by macroturbulence (e.g., Ryans
et al., 2002) may play a role in these discrepancies.

The mass-loss rate is known to within a factor of a few when more recent
estimates are compared. The spread in values is within an order of magnitude
when a wider range of estimates and their uncertainties from the literature
are considered.

The terminal velocity $v_{\infty}$ shows an apparent trend with rather high
values preferred in earlier studies and lower values in more current ones,
even when we consider that Giménez-García et al. (2016) and Sander et al.
(2018) were based on the same underlying model code. The $\beta$ parameter has
usually not been derived, but was set to a specific value based on assumptions
about the physics of line-driven winds.

Figure 19 compares different descriptions of the stellar wind velocity around
HD 77581. The profiles based on Dupree et al. (1980), McCray et al. (1984),
and Watanabe et al. (2006) are $\beta$-laws (Eq. 1). Krtička et al. (2012)
used a more complex polynomial expansion in their study. In the case of
Giménez-García et al. (2016), a $\beta$-law is connected to an atmosphere
solution, while in the three solutions given by Sander et al. (2018), the
velocity profile is determined as a solution of a hydrodynamically consistent
atmosphere model describing the wind stratification.

###  5.7 Wake structures

Charles et al. (1976) reported observations of Vela X-1 with the MSSL X-ray
instrument. They found indications in their data for increased absorption in a
trailing accretion wake and possibly in a leading bow shock, noting that these
structures varied from cycle to cycle and probably during each cycle.
Referring to these results, Conti (1978) noted that a wake like this would
mean that the wind velocity at the distance of the source would be on the
order of the orbital velocity, that is, at much lower values than assumed in
some later publications (see Fig. 19). The existence of this wake is also the
commonly accepted main contribution to the high column densities observed at
late orbital phases (Sect. 3.1.3 and Fig 4).

Indications for a photoionization wake (see Sect. 4.1.2) have been found by
Kaper et al. (1994) from high-resolution optical spectroscopy of the hydrogen
Balmer lines and He lines and by van Loon et al. (2001) from orbital
modulation through additional absorption of UV lines (see Sect. 3.2.4). Barziv
et al. (2001) also explained part of their spectroscopic data with the
presence of this wake (see Sect. 3.2.2). From their simulations of stellar
wind disruption, Blondin et al. (1990) noted that for a photoionization wake
to contribute significantly to the integrated column density, the wind would
have to be ionized most of the way to the companion star.

In summary, the trailing wakes caused by the movement of the neutron star
through the dense stellar wind and the added effect of the ionizing X-ray
radiation play an important role in the modulation of the observed emission at
many wavelengths. When models are compared to data, it is important to keep in
mind that these structures can be expected to significantly vary from one
orbital cycle to the next, and that the photoionization wake is also strongly
affected by the significant variations in the X-ray flux (Sect. 3.1.5).

###  5.8 Strength of the magnetic field of the neutron star

When we identify the 25 keV line that is visible in the X-ray spectrum of Vela
X-1 (Section 3.1.7) as the fundamental cyclotron line, we can use Eq. 7 to
infer a magnetic field strength in the line-forming region of
$2\times(1+z)10^{12}$ G. As has been shown by Fürst et al. (2010) and La
Parola et al. (2016), the first harmonic (at about 55 keV) shows a correlation
with luminosity, indicating that the line-forming region changes in altitude
above the surface of the neutron star. Using this correlation and theoretical
calculations by Becker et al. (2012), Fürst et al. (2010) inferred a surface
magnetic field of $B_{\star}=2.6\times 10^{12}$ G. However, this value should
be taken only as an estimate because it depends on the exact physical
conditions within the accretion column, in particular, on the dominating
process of decelerating the in-falling material, as well as on the mass and
radius of the neutron star.

In principle, we could also infer constraints on the magnetic field strength
from the observed torques, as has been done for other accreting X-ray sources.
This type of estimate depends on many model assumptions (see Section 4.2.2),
however, especially for a source with a potentially variable accretion
geometry, and thus does not add information for Vela X-1.

Table 8: Selected estimates of the maximum spin period or frequency changes reported in the literature, as given in the reference. For comparison, they are converted into $\dot{P}$ in ms per day. Reference | Values given in reference | Units used in |  Expressed as $\dot{P}$ in ms / day
---|---|---|---
| Spin-up | Spin-down | reference | Spin-up | Spin-down
Nagase et al. (1984b) | $-4\times 10^{-3}$ | $+3\times 10^{-3}$ |  $\dot{P}/P$ yr-1 | $-3.1$ | $+2.3$
Boynton et al. (1984) |  up to $\pm 5.8(1.4)\times 10^{-3}$ |  $\dot{\Omega}/\Omega$ yr-1 |  up to $\pm 4.5(1)$
Deeter et al. (1989) | +2.3(0.5) |  $-$2.4(0.6) |  $\dot{\Omega}$ in pico-rad s-1 |  $-$2.5(0.6) | +2.6(0.7)
İnam & Baykal (2000) | $+7.1(2.2)\times 10^{-9}$ | $-6.7(2.2)\times 10^{-9}$ |  $\dot{\nu}$ in Hz day-1 |  $-$0.6(0.2) | +0.5(0.2)
Liao et al. (2020) | $+2.3(1.3)\times 10^{-13}$ | $-1.5(8)\times 10^{-13}$ |  $\dot{\nu}$ in s-2 |  $-$1.6(0.9) | +1.0(0.6)

###  5.9 Observed torques on the neutron star

Changes in the observed spin period, or equivalently, the spin frequency
(Section 3.1.8), can in principle be used to derive information on the torques
acting on the neutron star through interactions in the magnetosphere (Section
4.2.2). It is important to note, however, that the reported periods or
frequencies are often average values over timescales of hours, days, and
sometimes many weeks, which averages out short-term fluctuations. The reported
changes are thus at least on the same timescales, or between data points even
farther apart. In contrast, model predictions are instantaneous. A comparison
between observed and predicted changes therefore implies that the
characteristic timescale of a change in spin-period differences is on the same
timescale as the observations. Otherwise, even for a correct model, the
absolute values of the observed differential change will always be lower than
the differential change inferred from the model. On the other hand, attempts
to derive pulse periods on shorter time intervals within long observations
will be affected by the pulse-to-pulse variations in the observed pulse
profile (Section 3.1.9).

Table 8 collects a number of reported spin-period or frequency changes in the
literature, without aiming at completeness. It demonstrates rather consistent
values with the largest observed changes in pulse period on the order of a few
milliseconds per day.

##  6 Outlook to future observations or studies

In this section we consider possible future observations or studies that would
help understand the Vela X-1 system better and thus allow deeper insights into
the physics of this archetype, which may shed light on many other systems as
well. The following ideas do not claim to be an exhaustive list, but reflect
topics raised while compiling this overview.

###  6.1 Evolutionary scenario and origin

The precise evolutionary stage of the donor star is still unclear. Is the blue
supergiant expanding toward becoming a red supergiant? Or is the donor already
evolving blueward and might become a Wolf-Rayet star? The precise
determination of the stellar properties will be key to answering these
questions and constrain evolutionary modeling efforts. These will need to be
provided by analyzing the donor with hydrodynamically consistent models that
sufficiently reproduce the spectrum, but also take the X-ray irradiation of
different parts of the wind into account. The numerical costs currently
prevent full 3D scenarios from being developed, but sufficient approximations
using multiple cones could become available in the near future and extend the
first efforts by Krtička et al. (2012) and Sander et al. (2018).

In order to test the association of Vela X-1 with the Vela OB1 association
(Sect. 5.1.2), it might be attempted to retrace the system trajectory using a
Galactic potential. The difficulty lies in an accurate determination of the
center-of-mass velocity of the binary, which requires low uncertainties for
both the tangential proper motion and the absolute RV. While the first can be
obtained with high precision from Gaia data, the wide absorption lines
affected by line infilling from the stellar wind and pulsations make the
second a challenge that requires multiepoch high-resolution spectroscopy.

###  6.2 Multiwavelength observations

It is essential to coordinate systematic multiwavelength observations to
follow variations on timescales of days or faster. Excellent data exist at all
wavelengths, but for the most part, they were not coordinated in the past and
are thus not able to connect the short-term variations across the spectral
bands.

Another angle that could be further explored are possible correlations between
the longer-term X-ray variations that have been captured for decades now by
X-ray monitors, and the known significant variations in optical brightness
around the mean curve determined by the orbit. In the area of robotic survey
telescopes, obtaining regular short snapshots of a star as bright as HD 77581
should be rather easy. We also note that Vela X 1 has occasionally been in the
field of view of the TESS111111https://heasarc.gsfc.nasa.gov/docs/tess/
mission with its excellent photometry.

The advances of optical and near-IR interferometry, especially recently with
VLTI/Gravity (Gravity Collaboration et al., 2017), allow resolving HMXBs in
the Galaxy at submilliarcsecond resolution and thus at the spatial scale where
accretion takes place (Waisberg et al., 2017). Spectral differential
interferometry can provide spatial information on scales as small as $\sim$1
to 10 $\mu$as, but the only published interferometric results for Vela X-1
(Choquet et al., 2014, see Sect 3.2.5) were unable to detect differential
visibility signatures beyond the noise level. At the distance of Vela X-1, a
resolution of 10 $\mu$as would translate into approximately 10% of the orbital
separation, resolving the scale of the accretion radius (see Fig. 2) and the
level of the uncertainty of the radius of HD 77581 (see Tab. 4). Reaching the
scale of the magnetosphere would require further improvements by at least
another order of magnitude, which appears beyond the scope of improvements in
interferometry foreseen in the nearer future (Eisenhauer, 2019).

In X-rays, no interferometry mission is in direct planning. However, mission
concepts exist, and high-mass X-ray binaries feature prominently among the
possible science cases (Uttley et al., 2019). In particular, even a single
spacecraft interferometer could reach a resolution of $\sim$100 $\mu$as and
thus resolve the components of a bright high-mass X-ray binary such as Vela
X-1 in an observation of a few hours (Uttley et al., 2020). Higher resolution
would likely require formation flight of multiple spacecrafts, but would allow
us to determine the location of different emission components.

In order to study the winds on timescales of seconds, data from the fast and
sensitive optical spectrographs available today will be required from
observations that ideally would cover different orbital phases and also longer
timescales.

At the X-ray energy range, X-ray calorimeters such as that on board the X-ray
Imaging and Spectroscopy Mission (XRISM) (XRISM Science Team, 2020) or the
X-ray Integral Field Unit (Barret et al., 2018) of the Athena observatory
(Barcons et al., 2017) will definitely usher in a new era on the X-ray line
diagnostics of the source in combination with more precise atomic data from
laboratory studies. Compared to the current grating spectrometers, they will
allow studying lines in detail on shorter timescales and thus will enable us
to better follow the dynamics.

The advent of X-ray polarimetry with the upcoming Imaging X-ray Polarimetry
Explorer (IXPE)121212https://ixpe.msfc.nasa.gov/ mission (Sgrò & IXPE Team,
2019) and later in the decade, the planned enhanced X-ray Timing and
Polarimetry (eXTP) mission (in’t Zand et al., 2019), will offer a
complementary opportunity of obtaining information about the geometry of the
intrinsic X-ray emission, as well as about the system geometry and structures
in the stellar wind when compared with detailed models (e.g., Kallman et al.,
2015; Caiazzo & Heyl, 2021, see also below). Additional information can be
obtained from hard X-ray polarimetry, such as has been proposed for the
balloon experiment XL-Calibur (Abarr et al., 2021), which will also cover the
crucial CRSF energy range.

The unexpected detection of Vela X-1 as a radio source has opened a
tantalizing new window to explore the system properties. At the time when this
article was completed, further radio observations have taken place that were
coordinated with X-ray observatories, but the results remain to be examined.

###  6.3 Tracks for future models

Future modeling efforts are still required to account for the features
observed in Vela X-1 and distinguish those that are found in other wind-fed
HMXBs from those that are specific to this system. This will be achieved both
by neatly tailoring the spatial scales together, from the stellar photosphere
all the way down the accretion columns at the neutron magnetic poles, and by
including additional physical ingredients at each level.

In Vela X-1, the orbital eccentricity has long been considered as negligible
by modelers. However, the increasingly improving performances reached by time-
resolved spectroscopy invite us to investigate the effect it could have on the
accretion flow (e.g., presence of a wind-captured disk, shape of the accretion
wake, and mass transfer rate). Eccentricity induces a periodic modulation at
the orbital period that might produce systematic differences as a function of
the orbital phase. This is worth looking for in observations and simulations.

The gravitational pull of the neutron star is an essential factor in shaping
the wind flow. Thus, improved simulations of the flow compared to
observational diagnostics, expanding on the approach of Manousakis et al.
(2012) (see Sect. 5.3), may be able to provide alternative, albeit model-
dependent constraints on the masses of the system components.

Owing to the Roche potential and to its own spin, the donor star is not
spherical, which means that the temperature of its photosphere is not uniform.
For stars with a radiative envelope such as HD 77581, gravity darkening is
well described by the von Zeipel law, which connects local effective gravity
and temperature (von Zeipel, 1924). Given the high stellar filling factor, if
the star were in synchronous rotation, the local effective temperature of the
photosphere near the inner Lagrangian point could be up to $\sim$50% lower
than at the poles (Espinosa Lara & Rieutord, 2012). Nonsynchronous rotation
mitigates this effect, but given how dependent line-driven winds are on the
radiative field, the local mass-loss rate probably varies throughout the
stellar photosphere (Friend & Castor, 1982; Hadrava & Čechura, 2012; Čechura &
Hadrava, 2015). An accurate description of the 3D wind-launching mechanism
requires taking the nonuniform conditions in the photosphere into account,
along with clump formation. This research track echoes the questions brought
up by the multiple evidence in favor of tidally induced nonradial oscillations
of the stellar photosphere.

Another aspects this relates to is the pressing need for multidimensional
models of radiative transfer for UV (e.g., for resonant line-absorption of
stellar photons by the outflowing wind) and for X-rays (e.g., for
photoionization and for the altered conditions in the X-ray irradiated stellar
photosphere). The computational cost of radiative transfer in 3D is very high,
however, therefore it does not seem feasible, even in the longer term, to
fully couple radiation and matter in time-dependent radiative-hydrodynamics
simulations without resorting to simplifying assumptions. Laboratory X-ray
spectroscopic data will be an invaluable compass to determine which hypotheses
are the least harmful (Hell et al., 2016). Notwithstanding the radiative
feedback on the flow dynamics, radiative computations could be performed on
simulation snapshots obtained by using an elementary representation of
radiation-matter coupling (e.g., for wind launching and for the inhibited
line-acceleration beyond a certain critical ionization parameter). It will be
of tremendous importance to interpret the absorbing column density observed by
the upcoming generation of X-ray instruments. The process could even be
iterated to obtain reasonably self-consistent radiative-hydrodynamics
simulations.

In the light of the upcoming X-ray polarimetry instruments, another promising
perspective is the possibility of gaining insights into the geometry and
properties of the system based on detailed modeling of polarized radiation
(see, e.g., Kallman et al., 2015; Caiazzo & Heyl, 2021, further described in
Sect. 3.1.6). Future studies may build on these approaches, considering also
the effect of the emitted X-rays on the ambient medium in multiscale X-ray
radiative transfer computations in order to provide an accurate view of the
polarization of the emitted radiation and the system parameters driving its
observed variations. Polarization signatures like this can be probed with
upcoming missions such as IXPE and eXTP (see Sect. 6.2).

In order to obtain an accurate description of the accretion process itself, we
eventually need better magnetohydrodynamics models, guided by the observed
pulse profiles and the changing neutron star spin period. The latter
encapsulates the information on the coupling between the accreted plasma and
the magnetosphere, while the former tells us about the X-ray emission and
scattering processes in the immediate vicinity of the neutron star surface.
Regarding the accretion-induced torques, significant progress has been made
for disk-magnetosphere coupling around millisecond pulsars (Parfrey &
Tchekhovskoy, 2017) and around T Tauri stars and magnetized white dwarfs in
cataclysmic variables (Romanova et al., 2003; Zanni & Ferreira, 2009). They
pave the way for a future model adapted to a system in which the accretor is a
fast-spinning highly magnetized neutron star and in which wind-captured disks
are transient at best. Alternatively, the quasi-spherical subsonic settling
models introduced by Shakura et al. (2012) might help to connect the evolution
of the X-ray luminosity and of the spin period of the neutron star.

In order to connect X-ray luminosity and mass accretion rate, the accretion
column model needs to be compared to the pulse profiles and their variations
in Vela X-1, however. The complex pulse profile structure and the wealth of
observational data at various luminosity levels promise significant diagnostic
power, first to distinguish plausible emission geometries, and in a further
step, to connect observed short-term variations to changes in the emission
regions and/or in the intervening material.

###  6.4 Key parameters for future studies

Table 9: Key parameters of the Vela X-1 system for future studies. Because the scatter in the values found in the literature is often significant, we refrain from giving specific values here, but refer to the discussion in the corresponding section(s) and provide a rough assessment of the robustness of the knowledge of the parameters. We also note sections describing system elements, observational diagnostics or modeling, where these parameters are relevant. Parameter | Current knowledge | See also sections | Robustness, uncertainty
---|---|---|---
System
Distance |  Sect. 5.1 |  3.1.1, 3.2.3 |  Quite robust, known to $\sim$5%
Orbital phase parameters ($T_{\textrm{ecl}}$ or $T_{\pi/2}$, $P_{\textrm{orb}}$) |  Sect. 5.2, Tab. 3 |  3.1.3, 3.1.5, 3.2.1, 3.2.2, 3.2.4 |  Very robust, uncertainty $<$1%
Eccentricity $e$, longitude of periastron $\omega$ |  Sect. 5.2, Tab. 3 |  3.2.1, 4.2.1, 4.5.1, 4.5.2 |  Robust, known to $\sim$1%,
Inclination $i$ |  Sect. 5.3, Tab. 4 |  2.2, 4.2.1, Fig. 17 |  Within $\sim$20∘ wide interval
HD 77581
Stellar mass $M_{\mathrm{opt}}$ |  Sect. 5.3, Tab. 4 |  3.2.2, 4.6 | Systematics often 10–-20%
Stellar radius $R_{\mathrm{opt}}$ |  Sect. 5.3, Tab. 4 |  2.1, 4.2.1 |  Known to $<$10%
Eff. temperature $T_{\star}$ and surface gravity $g_{\star}$ |  Sect. 5.5, Tab. 5 |  2.1, 3.2.3 | Dep. on spectral model codes
Mass-loss rate $\dot{M}$ |  Sect. 5.6, Tab. 7 |  2.1, 3.3, 4.5.1, 4.6 | Dep. on distance & spectral models
Terminal wind velocity $\varv_{\infty}$ |  Sect. 5.6, Tab. 7 |  2.1, 3.2.3, 3.3 |  Known to $\sim$5%, phase-dependent
Neutron star
Mass $M_{\mathrm{NS}}$ |  Sect. 5.3, Tab. 4, Fig. 17 |  2.2, 2.4, 3.2.2, 4.2.1 | Only known to a few 10%
Magnetic field $B$ |  Sect. 5.8 |  2.3, 2.4, 3.1.7, 4.2.2, 4.4 |  $\sim$10% in CRSF formation zone

Finally, we give an overview of key parameters of the Vela X-1 system that it
would be useful to study and constrain further, in our opinion. It is
important to note that these parameters are typically not independent of each
other. Refining existing determinations or having independent constraints for
any of these parameters will affect the ranges of several other parameters and
thus lead to a better-defined baseline for the physics of the system. The
parameters are summarized in Tab. 9.

The distance to the Vela X-1 system is now better determined based on Gaia
EDR3 data and consistent with other methods (Sect. 5.1). Still, the range
determined in this work leads to a $\sim$30% uncertainty in absolute
luminosities just from the distance factor.

More recent orbital ephemerides agree very well in general, but there is a
slight difference between the last detailed values in the literature
(Kreykenbohm et al., 2008; Falanga et al., 2015) that is significant at the
quoted uncertainties (Sect. 5.2, Tab. 3). An updated ephemeris with a recent
zero-point might resolve this difference and would avoid the systematic
uncertainty from extrapolating over many years. The orbital eccentricity,
which is frequently ignored in modeling efforts, plays an important role for
the varying shape of the mass donor (Sect. 3.2.1) and the mass transfer in
this system, where the donor star is close to filling its Roche lobe (Sect.
4.2.1). For an exact description of these effects, the longitude of periastron
$\omega$ also needs to be known well.

While it is well-constrained for an X-ray binary, the inclination $i$ of the
system remains a strong source of uncertainty for the determination of
absolute values for the masses of the two partners and the radius of HD 77581
(Sect. 5.3, Tab. 4). Conversely, if the supergiant size could be tightly
constrained, the inclination and thus other system parameters would be better
constrained.

The stellar parameters may also be improved in two ways. First, through
detailed SED modeling of high-resolution spectroscopic data, which would yield
more accurate and precise values of the effective temperature and gravity and
would investigate possible composition anomalies. Second, using well-
calibrated photometry and/or spectrophotometry to independently estimate the
effective temperature from the Balmer jump.

With its multiple constraints from the orbit on the one hand and the spectra
and photometric data on the other hand, HD 77581 acts as a testbed for our
understanding of the winds and mass loss of hot, massive stars. To launch a
stellar wind, the ratio of luminosity to mass as well as the radius or the
effective temperature of HD 77581 are vital ingredients. With the considerable
change in the presumed luminosity of the donor star due to the revised
spectral classification, the absolute mass-loss rate needs to be redetermined,
while the terminal wind velocity must be known much better due to its direct
spectral imprint. New quantitative studies with hydrodynamically consistent
atmosphere models building on these updated values are required to determine
whether a coherent solution for the observed appearance and the current
concept of radiation-driven winds for OB stars can be found. The solution for
HD 77581 will also mark an important constraint for the still highly uncertain
winds of binary evolution products.

Because Vela X-1 is frequently quoted as an example of a heavy neutron star,
having a firmer constraint of $M_{\mathrm{NS}}$ would be of high interest
beyond the studies of accreting X-ray pulsars per se. As Tab. 4 and Fig. 17 in
Sect. 5.3 demonstrate, the effective uncertainty of comparing different
determinations remains high, however, and the values depend on other
assumptions about the system, especially the inclination $i$.

The magnetic field of the neutron star is another key parameter that drives
the accretion physics (Sect. 2.3) and the physics of the emission region
(Sect. 2.4), among other elements. While a good estimate of the magnetic field
strength in the main X-ray emission region exists from cyclotron resonance
scattering features (Sect. 3.1.7) constrained to about 10%, it is not clear if
this value can be trivially extrapolated to apply everywhere based on a simple
dipole geometry, however.

##  7 Summary

In this extensive review of a well-known X-ray binary system that shares many
properties with other often less well studied systems, we have first described
the essential elements and basic physics of this archetypical system. Second,
we summarized the observational diagnostics and results from the radio to the
hard X-rays, also including caveats that need to be considered. Third, we
discussed the options and challenges of modeling this X-ray binary in detail
and reviewed the model efforts untertaken so far. Fourth, we summarized the
known properties of the system, which often have a larger spread in their
reported values than is commonly assumed, including essential parameters such
as the neutron star mass or the velocity of the stellar wind. Fifth, we
demonstrated that the system may be at least close to RLOF at periastron and
that the eccentricity is high enough to imply a significant variation of mass
transfer with orbital phase. Sixth, we derived an updated distance estimate to
the Vela X-1 system based on the latest Gaia Data Release131313Gaia EDR3,
https://www.cosmos.esa.int/web/gaia/release. It has tighter margins than
before. This distance estimate is closer to previous estimates and corrects
for the systematic shift that affects the distance given in Bailer-Jones et
al. (2018). Seventh, we revisited the proposed connection to the Vela OB1
association and derived an updated Gaia distance estimate to this association.
Eight, we newly determined the stellar parameters of HD 77581, finding B0.2 Ia
as the most probable spectral type and luminosity class. Ninth, we provided an
overview of opportunities for future observations and modeling efforts to
improve our knowledge of the physics of this specific system and by
extrapolation of related, less well studied systems, together with an overview
of especially relevant key parameters. We encourage similar efforts to group
and discuss the actual knowledge about system parameters also for other well-
known systems that are used as prototypes in the literature in order to create
a reliable framework for future studies.

###### Acknowledgements.

We thank the referee for positive and succinct, but very pertinent comments
that identified some remaining gaps in the presented material and led us to
several additions to the text of this publication. IEM has received funding
from the Research Foundation Flanders (FWO), from the European Union’s Horizon
2020 research and innovation program under the Marie Skłodowska-Curie grant
agreement No 665501 and from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No 863412). SMN acknowledges funding by the Spanish Ministry MCIU
under project RTI2018-096686-B-C21 (MCIU/AEI/FEDER, UE), co-funded by FEDER
funds and by the Unidad de Excelencia María de Maeztu, ref. MDM-2017-0765. VG
is supported through the Margarete von Wrangell fellowship by the ESF and the
Ministry of Science, Research and the Arts Baden-Württemberg. IEM and PK
acknowledge support from the ESA/ESAC Faculty Visiting Scientist Programme.
JMA acknowledges support from the Spanish Government Ministerio de Ciencia
through grant PGC2018-$095\,049$-B-C22. FJE acknowledges support from ESCAPE -
The European Science Cluster of Astronomy & Particle Physics ESFRI Research
Infrastructures, which received funding from the European Union’s Horizon 2020
research and innovation programme under Grant Agreement no. 824064. JvdE is
supported by a Lee Hysan Junior Research Fellowship awarded by St. Hilda’s
College.
The simulations were conducted on the Tier-1 VSC (Flemish Supercomputer Center
funded by Hercules foundation and Flemish government). GraphClick (© Arizona
Software) and WebPlotDigitizer (© Ankit Rohatgi) have been used to digitise
data from figures in older publications.
We thank Eva Laplace for very helpful comments on stellar evolution in binary
systems. PK thanks Mercedes Hainzl for help in digitizing older data.

## References

  * Abarr et al. (2021) Abarr, Q., Awaki, H., Baring, M. G., et al. 2021, Astroparticle Physics, 126, 102529
  * Abarr et al. (2020) Abarr, Q., Baring, M., Beheshtipour, B., et al. 2020, ApJ, 891, 70
  * Alsing et al. (2018) Alsing, J., Silva, H. O., & Berti, E. 2018, MNRAS, 478, 1377
  * Amato et al. (2021) Amato, R., Grinberg, V., Hell, N., et al. 2021, A&A subm.
  * Ammann & Mauder (1978) Ammann, M. & Mauder, H. 1978, Mitteilungen der Astronomischen Gesellschaft Hamburg, 43, 219
  * Araya & Harding (1999) Araya, R. A. & Harding, A. K. 1999, ApJ, 517, 334
  * Arons et al. (1987) Arons, J., Klein, R. I., & Lea, S. M. 1987, ApJ, 312, 666
  * Arons & Lea (1976) Arons, J. & Lea, S. M. 1976, ApJ, 207, 914
  * Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58
  * Bar-Shalom et al. (2001) Bar-Shalom, A., Klapisch, M., & Oreg, J. 2001, J. Quant. Spec. Radiat. Transf., 71, 169
  * Barcons et al. (2017) Barcons, X., Barret, D., Decourchelle, A., et al. 2017, Astronomische Nachrichten, 338, 153
  * Barret et al. (2018) Barret, D., Lam Trong, T., den Herder, J.-W., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10699, Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, S. Nikzad, & K. Nakazawa, 106991G
  * Barziv et al. (2001) Barziv, O., Kaper, L., Van Kerkwijk, M. H., Telting, J. H., & Van Paradijs, J. 2001, A&A, 377, 925
  * Basko & Sunyaev (1975) Basko, M. M. & Sunyaev, R. A. 1975, A&A, 42, 311
  * Bautz et al. (1983) Bautz, M., Howe, S., Gorecki, A., et al. 1983, ApJ, 266, 794
  * Becker et al. (2012) Becker, P. A., Klochkov, D., Schönherr, G., et al. 2012, A&A, 544, A123
  * Becker & Wolff (2005a) Becker, P. A. & Wolff, M. T. 2005a, ApJ, 621, L45
  * Becker & Wolff (2005b) Becker, P. A. & Wolff, M. T. 2005b, ApJ, 630, 465
  * Becker & Wolff (2007) Becker, P. A. & Wolff, M. T. 2007, ApJ, 654, 435
  * Becker et al. (1978) Becker, R. H., Rothschild, R. E., Boldt, E. A., et al. 1978, ApJ, 221, 912
  * Belczynski et al. (2012) Belczynski, K., Bulik, T., & Fryer, C. L. 2012, arXiv e-prints, arXiv:1208.2422
  * Bessell et al. (1975) Bessell, M. S., Vidal, N. V., & Wickramasinghe, D. T. 1975, ApJ, 195, L117
  * Bildsten et al. (1997) Bildsten, L., Chakrabarty, D., Chiu, J., et al. 1997, ApJS, 113, 367
  * Blaauw (1961) Blaauw, A. 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Blaha & Humphreys (1989) Blaha, C. & Humphreys, R. M. 1989, AJ, 98, 1598
  * Blondin et al. (1990) Blondin, J. M., Kallman, T. R., Fryxell, B. A., & Taam, R. E. 1990, ApJ, 356, 591
  * Blondin & Raymer (2012) Blondin, J. M. & Raymer, E. 2012, ApJ, 752, 30
  * Blondin et al. (1991) Blondin, J. M., Stevens, I. R., & Kallman, T. R. 1991, ApJ, 371, 684
  * Bondi & Hoyle (1944) Bondi, H. & Hoyle, F. 1944, MNRAS, 104, 273
  * Börner et al. (1987) Börner, G., Hayakawa, S., Nagase, F., & Anzer, U. 1987, A&A, 182, 63
  * Boynton et al. (1986) Boynton, P. E., Deeter, J. E., Lamb, F. K., & Zylstra, G. 1986, ApJ, 307, 545
  * Boynton et al. (1984) Boynton, P. E., Deeter, J. E., Lamb, F. K., et al. 1984, ApJ, 283, L53
  * Bozzo et al. (2020) Bozzo, E., Ducci, L., & Falanga, M. 2020, arXiv e-prints, arXiv:2012.00392
  * Bozzo et al. (2008) Bozzo, E., Falanga, M., & Stella, L. 2008, ApJ, 683, 1031
  * Bozzo et al. (2016) Bozzo, E., Oskinova, L., Feldmeier, A., & Falanga, M. 2016, A&A, 589, A102
  * Bozzo et al. (2009) Bozzo, E., Stella, L., Vietri, M., & Ghosh, P. 2009, A&A, 493, 809
  * Braun & Langer (1995) Braun, H. & Langer, N. 1995, A&A, 297, 483
  * Bulik et al. (1995) Bulik, T., Riffert, H., Mészáros, P., et al. 1995, ApJ, 444, 405
  * Burnard et al. (1983) Burnard, D. J., Arons, J., & Lea, S. M. 1983, ApJ, 266, 175
  * Bussard et al. (1985) Bussard, R. W., Meszaros, P., & Alexander, S. 1985, ApJ, 297, L21
  * Caballero et al. (2011) Caballero, I., Kraus, U., Santangelo, A., Sasaki, M., & Kretschmar, P. 2011, A&A, 526, A131
  * Caiazzo & Heyl (2021) Caiazzo, I. & Heyl, J. 2021, MNRAS, 501, 109
  * Cappallo et al. (2017) Cappallo, R., Laycock, S. G. T., & Christodoulou, D. M. 2017, PASP, 129, 124201
  * Castor et al. (1975) Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157
  * Castor & Lamers (1979) Castor, J. I. & Lamers, H. J. G. L. M. 1979, ApJS, 39, 481
  * Charles et al. (1976) Charles, P. A., Mason, K. O., Culhane, J. L., Sanford, P. W., & White, N. E. 1976, The X-ray variability of Vela X-1 (3U 0900-40)., Vol. 389, 629–642
  * Charles et al. (1978) Charles, P. A., Mason, K. O., White, N. E., et al. 1978, MNRAS, 183, 813
  * Cherepashchuk (1982) Cherepashchuk, A. M. 1982, Soviet Astronomy Letters, 8, 82
  * Choquet et al. (2014) Choquet, É., Kervella, P., Le Bouquin, J. B., et al. 2014, A&A, 561, A46
  * Coleiro & Chaty (2013) Coleiro, A. & Chaty, S. 2013, ApJ, 764, 185
  * Conti (1978) Conti, P. S. 1978, A&A, 63, 225
  * Corbet & Krimm (2013) Corbet, R. H. D. & Krimm, H. A. 2013, ApJ, 778, 45
  * Cutri & et al. (2013) Cutri, R. M. & et al. 2013, VizieR Online Data Catalog, 2328
  * Daugherty & Harding (1986) Daugherty, J. K. & Harding, A. K. 1986, ApJ, 309, 362
  * Davidson (1973) Davidson, K. 1973, Nature Physical Science, 246, 1
  * Davidson & Ostriker (1973) Davidson, K. & Ostriker, J. P. 1973, ApJ, 179, 585
  * Davies & Pringle (1981) Davies, R. E. & Pringle, J. E. 1981, MNRAS, 196, 209
  * de Kool & Anzer (1993) de Kool, M. & Anzer, U. 1993, MNRAS, 262, 726
  * de Lima et al. (2020) de Lima, R. C. R., Coelho, J. G., Pereira, J. P., Rodrigues, C. V., & Rueda, J. A. 2020, ApJ, 889, 165
  * de Mink et al. (2008) de Mink, S. E., Pols, O. R., & Yoon, S. C. 2008, in American Institute of Physics Conference Series, Vol. 990, First Stars III, ed. B. W. O’Shea & A. Heger, 230–232
  * Deeter et al. (1987) Deeter, J. E., Boynton, P. E., Lamb, F. K., & Zylstra, G. 1987, ApJ, 314, 634
  * Deeter et al. (1989) Deeter, J. E., Boynton, P. E., Lamb, F. K., & Zylstra, G. 1989, ApJ, 336, 376
  * Deeter et al. (1987) Deeter, J. E., Boynton, P. E., Shibazaki, N., et al. 1987, AJ, 93, 877
  * Doroshenko (2017) Doroshenko, R. 2017, PhD thesis, Eberhard Karls Universität Tübingen
  * Doroshenko et al. (2013) Doroshenko, V., Santangelo, A., Nakahira, S., et al. 2013, A&A, 554, A37
  * Doroshenko et al. (2011) Doroshenko, V., Santangelo, A., & Suleimanov, V. 2011, A&A, 529, A52
  * Driessen et al. (2019) Driessen, F. A., Sundqvist, J. O., & Kee, N. D. 2019, A&A, 631, A172
  * Ducci et al. (2009) Ducci, L., Sidoli, L., Mereghetti, S., Paizis, A., & Romano, P. 2009, MNRAS, 398, 2152
  * Dupree et al. (1980) Dupree, A. K., Gursky, H., Black, J. H., et al. 1980, ApJ, 238, 969
  * Eadie et al. (1975) Eadie, G., Peacock, A., Pounds, K. A., & Watson, M. 1975, MNRAS, 172, 35p
  * Edgar (2004) Edgar, R. 2004, New Astronomy Reviews, 48, 843
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368
  * Eisenhauer (2019) Eisenhauer, F. 2019, in The Very Large Telescope in 2030, 30
  * El Mellah et al. (2020a) El Mellah, I., Bolte, J., Decin, L., Homan, W., & Keppens, R. 2020a, A&A, 637, A91
  * El Mellah & Casse (2015) El Mellah, I. & Casse, F. 2015, MNRAS, 454, 2657
  * El Mellah & Casse (2017) El Mellah, I. & Casse, F. 2017, MNRAS, 467, 2585
  * El Mellah et al. (2020b) El Mellah, I., Grinberg, V., Sundqvist, J. O., Driessen, F. A., & Leutenegger, M. A. 2020b, A&A, 643, A9
  * El Mellah et al. (2019a) El Mellah, I., Sander, A. A. C., Sundqvist, J. O., & Keppens, R. 2019a, A&A, 622, A189
  * El Mellah et al. (2018) El Mellah, I., Sundqvist, J. O., & Keppens, R. 2018, MNRAS, 475, 3240
  * El Mellah et al. (2019b) El Mellah, I., Sundqvist, J. O., & Keppens, R. 2019b, A&A, 622, L3
  * Espinosa Lara & Rieutord (2012) Espinosa Lara, F. & Rieutord, M. 2012, A&A, 547, A32
  * Evans et al. (2018) Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4
  * Fabricius et al. (2020) Fabricius, C., Luri, X., Arenou, F., et al. 2020, arXiv e-prints, arXiv:2012.06242
  * Falanga et al. (2015) Falanga, M., Bozzo, E., Lutovinov, A., et al. 2015, A&A, 577, A130
  * Falkner (2018) Falkner, S. 2018, PhD thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg
  * Farinelli et al. (2012) Farinelli, R., Ceccobello, C., Romano, P., & Titarchuk, L. 2012, A&A, 538, A67
  * Farinelli et al. (2016) Farinelli, R., Ferrigno, C., Bozzo, E., & Becker, P. A. 2016, A&A, 591, A29
  * Fender et al. (2004) Fender, R. P., Belloni, T. M., & Gallo, E. 2004, MNRAS, 355, 1105
  * Fender & Hendry (2000) Fender, R. P. & Hendry, M. A. 2000, MNRAS, 317, 1
  * Foglizzo & Ruffert (1997) Foglizzo, T. & Ruffert, M. 1997, A&A, 320, 342
  * Forman et al. (1973) Forman, W., Jones, C., Tananbaum, H., et al. 1973, ApJ, 182, L103
  * Francischelli et al. (2002) Francischelli, G. J., Wijers, R. A. M. J., & Brown, G. E. 2002, ApJ, 565, 471
  * Fransson & Fabian (1980) Fransson, C. & Fabian, A. C. 1980, A&A, 87, 102
  * Fraser et al. (2010) Fraser, M., Dufton, P. L., Hunter, I., & Ryans, R. S. I. 2010, MNRAS, 404, 1306
  * Friedemann et al. (1996) Friedemann, C., Guertler, J., & Loewe, M. 1996, A&AS, 117, 205
  * Friend & Castor (1982) Friend, D. B. & Castor, J. I. 1982, ApJ, 261, 293
  * Fryer et al. (2013) Fryer, C. L., Belczynski, K., Berger, E., et al. 2013, ApJ, 764, 181
  * Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273
  * Fryxell & Taam (1988) Fryxell, B. A. & Taam, R. E. 1988, ApJ, 335, 862
  * Fryxell et al. (1987) Fryxell, B. A., Taam, R. E., & McMillan, S. L. W. 1987, ApJ, 315, 536
  * Fürst et al. (2010) Fürst, F., Kreykenbohm, I., Pottschmidt, K., et al. 2010, A&A, 519, A37+
  * Fürst et al. (2014) Fürst, F., Pottschmidt, K., Wilms, J., et al. 2014, ApJ, 780, 133
  * Gaia Collaboration (2018) Gaia Collaboration. 2018, VizieR Online Data Catalog, I/345
  * Gaia Collaboration et al. (2020) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2020, arXiv e-prints, arXiv:2012.01533