11institutetext: European Space Agency (ESA), European Space Astronomy Centre
(ESAC), Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid,
Spain 11email<EMAIL_ADDRESS>22institutetext: Centre for
Mathematical Plasma Astrophysics, KU Leuven, Celestijnenlaan 200B, 3001
Leuven, Belgium 33institutetext: Institute of Astronomy, KU Leuven, 3001
Leuven, Belgium 44institutetext: Institut de Planétologie et d’Astrophysique
de Grenoble, 414 Rue de la Piscine, 38400 Saint-Martin-d’Hères, France
55institutetext: Instituto de Física de Cantabria (CSIC-Universidad de
Cantabria), E-39005, Santander, Spain 66institutetext: Quasar Science
Resources S.L for European Space Agency (ESA), European Space Astronomy Centre
(ESAC), Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid,
Spain 77institutetext: Institut für Astronomie und Astrophysik (IAAT),
Universität Tübingen, Sand 1, 72076 Tübingen, Germany. 88institutetext: Armagh
Observatory and Planetarium, College Hill, Armagh BT61 9DG, Northern Ireland,
UK 99institutetext: Anton Pannekoek Institute for Astronomy, University of
Amsterdam, Science Park 904, NL-1098 XH Amsterdam, the Netherlands
1010institutetext: Department of Physics, Astrophysics, University of Oxford,
Oxford, UK 1111institutetext: Centro de Astrobiología, CSIC-INTA, Campus ESAC,
Camino bajo del castillo s/n, E-$28\,692$, Villanueva de la Cañada, Spain
1212institutetext: Vitrociset Belgium for European Space Agency (ESA),
European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, 28692
Villanueva de la Cañada, Madrid, Spain 1313institutetext: Aurora Technology BV
for European Space Agency (ESA), European Space Astronomy Centre (ESAC),
Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain

# Revisiting the archetypical wind accretor Vela X-1 in depth

Case study of a well-known X-ray binary and the limits of our knowledge

P. Kretschmar 11 I. El Mellah 223344 S. Martínez-Núñez 55 F. Fürst 66 V.
Grinberg 77 A. A. C. Sander 88 J. van den Eijnden 991010 N. Degenaar 99 J.
Maíz Apellániz 1111 F. Jiménez Esteban 1111 M. Ramos-Lerate 1212 E. Utrilla
1313

(Received 31 December 2020; accepted 22 April 2021)

###### Abstract

Context. The Vela X-1 system is one of the best-studied X-ray binaries because
it was detected early, has persistent X-ray emission, and a rich phenomenology
at many wavelengths. The system is frequently quoted as the archetype of wind-
accreting high-mass X-ray binaries, and its parameters are referred to as
typical examples. Specific values for these parameters have frequently been
used in subsequent studies, however, without full consideration of
alternatives in the literature, even more so when results from one field of
astronomy (e.g., stellar wind parameters) are used in another (e.g., X-ray
astronomy). The issues and considerations discussed here for this specific,
very well-known example will apply to various other X-ray binaries and to the
study of their physics.

Aims. We provide a robust compilation and synthesis of the accumulated
knowledge about Vela X-1 as a solid baseline for future studies, adding new
information where available. Because this overview is targeted at a broader
readership, we include more background information on the physics of the
system and on methods than is usually done. We also attempt to identify
specific avenues of future research that could help to clarify open questions
or determine certain parameters better than is currently possible.

Methods. We explore the vast literature for Vela X-1 and on modeling efforts
based on this system or close analogs. We describe the evolution of our
knowledge of the system over the decades and provide overview information on
the essential parameters. We also add information derived from public data or
catalogs to the data taken from the literature, especially data from the Gaia
EDR3 release.

Results. We derive an updated distance to Vela X-1 and update the spectral
classification for HD 77518. At least around periastron, the supergiant star
may be very close to filling its Roche lobe. Constraints on the clumpiness of
the stellar wind from the supergiant star have improved, but discrepancies
persist. The orbit is in general very well determined, but a slight difference
exists between the latest ephemerides. The orbital inclination remains the
least certain factor and contributes significantly to the uncertainty in the
neutron star mass. Estimates for the stellar wind terminal velocity and
acceleration law have evolved strongly toward lower velocities over the years.
Recent results with wind velocities at the orbital distance in the range of or
lower than the orbital velocity of the neutron star support the idea of
transient wind-captured disks around the neutron star magnetosphere, for which
observational and theoretical indications have emerged. Hydrodynamic models
and observations are consistent with an accretion wake trailing the neutron
star.

Conclusions. With its extremely rich multiwavelength observational data and
wealth of related theoretical studies, Vela X-1 is an excellent laboratory for
exploring the physics of accreting X-ray binaries, especially in high-mass
systems. Nevertheless, much room remains to improve the accumulated knowledge.
On the observational side, well-coordinated multiwavelength observations and
observing campaigns addressing the intrinsic variability are required. New
opportunities will arise through new instrumentation, from optical and near-
infrared interferometry to the upcoming X-ray calorimeters and X-ray
polarimeters. Improved models of the stellar wind and flow of matter should
account for the non-negligible effect of the orbital eccentricity and the
nonspherical shape of HD 77581. There is a need for realistic multidimensional
models of radiative transfer in the UV and X-rays in order to better
understand the wind acceleration and effect of ionization, but these models
remain very challenging. Improved magnetohydrodynamic models covering a wide
range of scales are required to improve our understanding of the plasma-
magnetosphere coupling, and they are thus a key factor for understanding the
variability of the X-ray flux and the torques applied to the neutron star. A
full characterization of the X-ray emission from the accretion column remains
another so far unsolved challenge.

###### Key Words.:

X-rays: individuals: Vela X-1 – X-rays: binaries – Stars: winds, outflows –
Accretion, accretion disks

##  1 Introduction

Vela X-1 (4U 0900-40) is an eclipsing high-mass X-ray binary (HMXB) at a
distance of about 2 kpc (Section 5.1). The binary system consists of an
accreting neutron star that emits X-ray pulses with a period of $\sim$283 s
(Sect. 3.1.8) in a close, mildly eccentric orbit with a period of $\sim$8.964
days (Section 5.2, Table 3) around the supergiant HD 77581, also known as GP
Vel, which is most frequently described as B0.5 Ia supergiant. Because the
orbital separation of $\sim 1.7R_{\star}$ (Section 5.2) is small, the
accreting neutron star is deeply embedded in the dense stellar wind of the
supergiant star, with a mass loss of about $10^{-6}M_{\sun}\,\mathrm{yr}^{-1}$
(Sect. 5.6). Figure 1 illustrates the system geometry. In this figure and in
the remainder of this paper, the orbital phases quoted are expressed with
respect to the mid-eclipse time $T_{\mathrm{ecl}}$ defined in Kreykenbohm et
al. (2008) as zero-point; the different ephemerides used in the literature are
discussed in Section 5.2. The orbital eccentricity is high enough for the
orbital phase 0.5 to be reached significantly after inferior conjunction of
the neutron star. This detail is sometimes ignored when the system with its
“almost circular” orbit is discussed.

Figure 1: (top panel) Sketch of the Vela X-1 system as seen from Earth based
on the ephemeris of Kreykenbohm et al. (2008), using an intermediate assumed
inclination within the known constraints (see Section 5.3 and Table 4). Points
mark the position of the neutron star at given orbital phases. (bottom panel)
Top view of the system.

The intrinsic X-ray luminosity of Vela X-1 is only moderate ($\sim 4\times
10^{36}\mathrm{erg\leavevmode\nobreak\ s}^{-1}$ Nagase et al., 1986), but
because of its proximity, it is still one of the brightest persistent point
sources in the X-ray sky. Together with a wide variety of interesting
phenomenology, this has made Vela X-1 a well-studied target ever since its
detection. It has frequently been called an “archetype” of wind-accreting
HMXB.

Because of its eclipsing nature, the extent of the donor star with respect to
the orbital separation is well constrained, and so is the mass transfer
mechanism. Because the orbit of the slowly spinning neutron star around its
blue supergiant companion is almost circular ($e\sim 0.0898$), the stellar
line-driven wind provides a continuous source of material to the compact
object. The sources of stochastic variability of the intrinsic X-ray emission
induced by, for instance, overdense small-scale regions in the wind (hereafter
called clumps), can be isolated from the orbital modulation. The cooperative
behavior of Vela X-1 from an observational and theoretical point of view has
made it a privileged target of study for modelers, as further detailed in
Sect. 4.

As an example of a high-mass neutron star, with a mass that has been
determined to be clearly higher than the Chandrasekar limit (Sect. 5.3), Vela
X-1 is also frequently used as an example for models of compact stars and
equation-of-state calculations or studies of neutron star mass distributions
(e.g., Gangopadhyay et al., 2013; Özel & Freire, 2016; Mafa Takisa et al.,
2017; Gedela et al., 2019). Over the years, however, the knowledge about many
system parameters has evolved, and different authors have used different
conventions to report their results. Even more importantly, quite different
assumptions about essential parameters of the system have been used by
different authors in the interpretation of their observations or in their
modeling, for instance, the velocity profile of the stellar wind, the nature
and shape of larger structures, clumps and their parameters, or the role of
ionization.

In this article we revisit a large sample of observational findings. We treat
the information in a consistent manner as far as possible, and we compare the
observations with information obtained from recent detailed modeling,
especially of the wind structure (Sander et al., 2018) and of the flow of
matter around the neutron star (El Mellah et al., 2018). The remainder of this
paper is structured as follows: Section 2 summarizes the different elements of
the Vela X-1 system in order to describe the framework in which observational
diagnostics and modeling studies take place. Section 3 gives an overview of
the observational diagnostics at different wavelengths and describes which
parts of the system they explore. Section 4 describes semianalytical and
numerical models of Vela X-1 or close analogs, and each time we pay particular
attention to how different physical ingredients are handled. Section 5
summarizes the knowledge about basic system properties, as found in the
literature, with an emphasis on clearly describing the differences between the
results obtained in different studies. In Section 6 we give our view on how
further observations and modeling efforts may help us to answer the open
questions and add to our knowledge of this prototype system. Finally, Section
7 summarizes the essential points we raised throughout this overview.

##  2 Elements of the Vela X-1 system

The rich observational phenomenology of an HMXB system such as Vela X-1 is the
result of a complex interplay of its various elements, where the relevant
physics cover more than six orders of magnitude in size, as indicated in
Figure 2. In the following we describe these elements from larger to smaller
size scales. For the sake of clarity, we introduce the main physical concepts
relevant to this and similar systems when we discuss observations and modeling
efforts in later sections. For a more complete discussion of stellar wind and
HMXB physics see, for example, Martínez-Núñez et al. (2017).

###  2.1 Supergiant and its stellar wind

HD 77581 looses significant mass through a line-driven wind whose launching
mechanism, first identified by Lucy & Solomon (1970) and Castor et al. (1975),
relies on the resonant line-absorption of UV photons by partly ionized metal
ions. The amount of mass lost and the intrinsic acceleration by the radiative
field of the supergiant evidently strongly depend on parameters such as
temperature, luminosity, radius, mass, and effective surface gravity, which
are in themselves not perfectly well known; see Sections 5.3 and 5.6. One
major factor of uncertainty is that the effective temperature of the
supergiant star (about 25,000 K, as found in most studies; see Table 5), lies
in a sharp transition called the bistability jump between slow winds on the
cool side and fast winds on the hot side (Vink et al., 1999).

In a spherically symmetric configuration, theoretical arguments and
observational insights support a velocity profile for line-driven winds that
follows a so-called $\beta$-law (Puls et al., 2008),

$v(r)=v_{\infty}\left(1-\frac{R_{\star}}{r}\right)^{\beta},$ (1)

with $R_{\star}$ the stellar photosphere radius, $v_{\infty}$ the terminal
wind speed, and $\beta$ a positive exponent: its value determines how quickly
the wind reaches $v_{\infty}$, with a more gradual acceleration for a higher
value of $\beta$. In the case of Vela X-1, quite different values for
$v_{\infty}$ ($<$400 to $\sim$1700 km s-1) and $\beta$ (0.8 to $\sim 2.2$)
have been reported by different authors, and the usefulness of a $\beta$-law
as a description has been questioned (Section 5.6).

In line-driven winds, internal shocks are prone to develop because of the
line-deshadowing instability (Lucy & White, 1980; Owocki & Rybicki, 1984). The
local density departs from the smooth wind value, with an overall density
ratio of $\sim$100 between the most and least dense regions at a given
distance from the donor star. In the case of Vela X-1, current work (see
Section 4.1 for details and references) suggests that at one stellar radius
above the photosphere, most of the wind mass is contained in clumps. This wind
clumping has often not been taken into account for discussions of the overall
wind structure in the system. The internal shocks in line-driven winds that
are responsible for the clump formation also produce X-rays, although the
X-ray luminosity due to this process is low ($\sim$1033erg s-1). The sizes of
clumps in stellar winds and in HMXBs have been estimated with widely varying
ranges in the past (see the corresponding discussion in Martínez-Núñez et al.,
2017). Recent simulations for hot and massive stars (Driessen et al., 2019)
suggest that most of the wind mass is contained in clumps that are expected to
appear as coherent structures of mass 1017-18 g and of a size of about 1% of
the stellar radius. This is about $10^{10}$ cm for HD 77581.

The presence of the orbiting neutron star strongly affects the outflowing
stellar wind through its gravity, the X-ray emission, and the orbital
movement. Bessell et al. (1975) first highlighted the anisotropic distribution
of material in the orbital plane induced by tidal forces, in agreement with
the observational indications for a trailing wake that was reported in the
early papers on Vela X-1 (see references in Conti, 1978). As further described
in Section 4.1, simulations of the Vela X-1 system (or close proxies) find a
wind focused toward the neutron star, forming an unsteady bow shock that leads
to an accretion wake. This wake trails the neutron star throughout its orbit
(e.g, Blondin et al., 1991; Manousakis et al., 2012).

The line-driven acceleration of the stellar wind is affected by the intense
X-ray emissions that are emitted from the immediate vicinity of the neutron
star. This modifies the ionization structure of the wind in which it is
embedded and thus the ability of the stellar radiative field to accelerate the
wind through resonant line-absorption (Hatchett & McCray, 1977). An additional
ionizing contribution, which is expected to be very low for Vela X-1, however,
may come from the X-rays that are produced in shocks in the wind. The local
degree of ionization in an optically thin wind in thermal balance can be
evaluated using the ionization parameter (Tarter et al., 1969),

$\xi=\frac{L_{\text{X}}}{nr_{\text{X}}^{2}},$ (2)

where $L_{\text{X}}$ is the X-ray luminosity, $n$ the local atomic number
density of the gas, and $r_{\text{X}}$ the distance to the X-ray point source.
Above a certain critical ionization parameter $\xi_{\text{crit}}$, most of the
elements responsible for the wind acceleration (e.g., C, N, O, and Fe) are
expected to be fully ionized and the wind acceleration is essentially
quenched. Depending on the shape of the irradiating X-ray spectrum and on the
details of the resonant line-absorption mechanism, the value of
$\xi_{\text{crit}}$ is expected to range between 102 and
104erg$\cdot$cm$\cdot$s-1. Shocks between the radiatively driven wind and
stagnant gas around the compact object can in principle lead to a trailing
accretion wake (Fransson & Fabian, 1980), but for Vela X-1 and similar
systems, this is expected to be a lesser contribution to the gas density
(Blondin et al., 1990). It is important to note, however, that the “on/off
switch” for wind acceleration set by the critical ionization parameter is a
simplifying assumption, as is further explained in the models that we present
in section 4.1.2. Moreover, wind clumping weakens the effect of X-ray
ionization from the accreting compact object through the increased
recombination inside the clumps (Oskinova et al., 2012).

Another large-scale feature that is common around massive hot stars are
corotating interaction regions (Mullan, 1984). None has been reported in Vela
X-1, however.

###  2.2 Mass transfer to the neutron star

A fraction of the mass lost by HD 77581 is captured and accreted by the
neutron star and fuels the intense X-ray emission. While the system is
generally labeled a “wind accretor”, the exact mechanisms of mass transfer in
the Vela X-1 system are not fully determined and have again become a topic of
more intense research; see Section 4.2.1. Uncertainties remain on the orbital
inclination, and to a lesser extent, on the mass ratio (Section 5.3),
therefore it can currently not be excluded in principle that HD 77581 would
fill its Roche lobe and thus lose mass through Roche-Lobe overflow (RLOF).
However, the moderate X-ray luminosity of Vela X-1, the long and erratically
varying spin period (Section 3.1.8), the absence of signatures for a permanent
accretion disk, and the stable orbital period (Falanga et al., 2015) make an
RLOF scenario very unlikely (see Section 4.2.1).

It is therefore generally assumed that the main mass transfer in Vela X-1
occurs through gravitational capture from the stellar wind. The length scale
at which the wind is beamed toward the accretor by its gravitational field (on
the order of the Roche-lobe radius of the neutron star) is much larger than
the extent of the neutron star itself or its magnetosphere (see Figure 2). In
the fast-wind limit, the so-called Bondi-Hoyle-Littleton (BHL) formalism
(Hoyle & Lyttleton, 1939; Bondi & Hoyle, 1944; Edgar, 2004) provides a robust
framework. In this description, the fraction of the wind with an impact
parameter lower than the accretion radius $R_{\text{acc}}$ will be captured,

$R_{\text{acc}}=2GM_{\text{NS}}/v_{\text{rel}}^{2},$ (3)

where $M_{\text{NS}}$ is the mass of the neutron star and $v_{\text{rel}}$ is
the relative speed of the wind with respect to the neutron star. In this case,
the mass accretion rate is approximately given by

$\dot{M}_{\text{acc}}=\frac{4\pi\rho\left(GM_{\text{NS}}\right)^{2}}{v_{\text{rel}}^{3}},$
(4)

where $\rho$ is the wind density at the orbital separation. Subsequently,
empirical refinements were proposed to correct $v_{\text{rel}}$ for slower
winds (Davidson & Ostriker, 1973). Still, the mass accretion rate is very
sensitive to $v_{\text{rel}}$ , and to a lesser extent, to $\rho,$ and a wide
range of X-ray luminosities can be derived from smaller changes in these
parameters, which are only relatively loosely constrained from observations
and modeling (Sections 4.1 and 5.6). In addition, the geometry of the problem
itself might differ from the planar picture of the BHL framework because
orbital effects are strong enough in Vela X-1 to significantly bend the flow
to a point that it could form a disk-like structure before it reaches the
magnetosphere.

Figure 2: Overview of characteristic length scales in the Vela X-1 system.
Some of them are relatively well determined, on the order of 10%, and can be
considered effectively fixed over the timescale of known observations. Other
scales, shown with different types of hatching or shading, strongly depend on
different parameters, as explained in the text. They may vary on short
timescales of sometimes over an order of a magnitude or more, and a variation
in one may drive changes in another. The ranges shown in this figure for these
parameters are only to be taken as indicative. For comparison, the orbital
separation corresponds to approximately 1/4 AU.

###  2.3 Accretion close to and within the magnetosphere

The mass transfer rate is an upper limit on the mass accretion rate at the
neutron star surface that eventually produces the observed X-ray luminosity. A
fraction of the captured material could either stall in intermediate regions
or might in principle be ejected, although there are no reports on outflows
from close to the X-ray source in Vela X-1 so far (but see Section 3.3). When
the flow reaches the magnetosphere, it is highly ionized by the intense X-ray
flux that is emitted from the immediate vicinity of the neutron star surface.
It behaves as a plasma, and its dynamic is entirely determined by the neutron
star magnetic field, which channels the material all the way down to the
neutron star magnetic poles (Davidson & Ostriker, 1973).

For a qualitative understanding of the different possible accretion regimes
within the magnetosphere, the following length scales play a role. First, the
corotation radius which is the distance at which the neutron star
magnetosphere, which rotates with the neutron star spin period, has the same
angular velocity as the local Keplerian velocity. The mass, size, and spin
period of the neutron star are fairly constant over timescales shorter than
years, therefore the corotation radius can be considered as fixed. Second, the
magnetospheric radius111For the sake of simplicity, we ignore the distinction
between magnetospheric and Alfvén radii. which is the distance to the neutron
star below which the magnetic field dominates the dynamics of the inflowing
plasma. A rigorous evaluation of this radius depends on the geometry of the
inflow, but it is typically estimated by comparing the magnetic pressure to
the ram pressure such that it increases with decreasing mass accretion rates.
Other factors such as mass, radius, and magnetic field strength are again
taken to be essentially fixed. Third, the circularization radius which is the
radius of the Keplerian orbit that has the same specific angular momentum as
the orbit of the accreted flow. Semianalytic and numerical estimates exist to
evaluate the upper limit of the specific angular momentum of the inflow, but
they vary by more than an order of magnitude between each other (Illarionov &
Sunyaev, 1975; Shapiro & Lightman, 1976; Wang, 1981; Livio et al., 1986).

Different estimates exist for the angular momentum accretion rate, but the
upper limit commonly reads

$l=\frac{1}{2}\Omega R_{\textrm{acc}}^{2},$ (5)

where $\Omega=2\pi/P$ is the angular orbital speed and $P$ is the orbital
period. We can then deduce the circularization radius,

$R_{\textrm{circ}}=\frac{1+q}{4}\left(\frac{R_{\text{acc}}}{a}\right)^{4}a,$
(6)

where $q$ is the ratio of the mass of the donor star to the mass of the
neutron star, and $a$ is the orbital separation. With realistic values for the
orbital separation, the mass ratio, and the accretion radius, we obtain
circularization radii in Vela X-1 that range between $\sim 5\times 10^{7}$ cm
and $\sim 2\times 10^{9}$ cm for wind speeds at the orbital separation that
range between 600 and 300 km$/$s.

In the case of pure wind mass transfer, the flow carries a negligible amount
of angular momentum, and the circularization radius is much smaller than the
other two length scales (fast-wind case). A bow shock forms ahead of the
neutron star with an opening angle smaller for higher Mach numbers of the
inflow (El Mellah & Casse, 2015). The flow is deflected by the shock, and
provided its impact parameter is smaller than the accretion radius, it returns
to the neutron star, essentially from the back hemisphere (for simulations of
BHL accretion onto a magnetic dipole, see Toropina et al., 2012; Lee et al.,
2014). If the mass accretion rate is high enough, the flow efficiently cools
down downstream of the shock through bremsstrahlung or Compton processes, and
the material freefalls at supersonic speeds onto the neutron star
magnetosphere (Burnard et al., 1983). This regime is highly unsteady and is
referred to as the direct accretion regime. Otherwise, the shock is adiabatic
and the flow forms a quasi-spherical shell of hot material that subsonically
settles onto the neutron star magnetosphere (Davies & Pringle, 1981; Shakura
et al., 2012). This regime has ramifications for different models depending on
the intensity of the plasma-magnetosphere coupling (referred to as strong,
moderate, and weak coupling in Shakura et al., 2012, 2018). In general, the
plasma penetrates the magnetosphere through the interchange instability, which
is the magnetohydrodynamics counterpart of the Rayleigh-Taylor instability
(Arons & Lea, 1976). For a more comprehensive description of the different
accretion regimes onto a neutron star magnetosphere in wind-fed HMXBs in
general, see Bozzo et al. (2008) and Martínez-Núñez et al. (2017).

For slower winds, the wind-RLOF mass transfer mechanism becomes more
realistic. The bow shock extends up to a significant fraction of the Roche
lobe of the neutron star ($\sim$10%). As detailed in Section 4.2.1, in this
regime, the accreted flow acquires and retains enough angular momentum through
the shock so that the circularization radius is larger than the magnetospheric
radius and a centrifugally maintained structure, called a wind-captured disk,
can form within the shocked region. A disk like this would still be truncated
at its inner edge by the magnetosphere, which is a few hundreds times larger
than the neutron star radius, and thus such a disk would not radiate
sufficiently at higher energies to be detected. Indirect evidence has emerged
in recent years, however, for the existence of transient disks in wind-fed
HMXBs in general (Hu et al., 2017; Taam et al., 2018; Taani et al., 2019) and
for in Vela X-1 in particular (Liao et al., 2020). In this configuration, the
disk and the magnetosphere are coupled by magnetic reconnection, the Kelvin-
Helmholtz instability, and Bohm diffusion (Ghosh & Lamb, 1979a). The ratio of
the corotation radius and the magnetospheric radius then controls accretion
and the transfer of angular momentum between the neutron star and the accreted
matter. In the most simple cases, plasma that rotates faster at the
magnetosphere than at the magnetic field lines that are anchored in the
neutron star can dissipate kinetic energy, accrete, and spin the neutron star
up, while at the other extreme, in the propeller regime, plasma that rotates
more slowly than the magnetic field lines at the magnetosphere is expected to
be expelled by the centrifugal force at the magnetocentrifugal barrier. This
reduces the angular momentum of the neutron star (Illarionov & Sunyaev, 1975).
Intermediate cases are discussed in Ghosh & Lamb (1979a) and many publications
based on this work (see Bozzo et al., 2009, for an overview and application to
different sources). Strong torques for spin-up and spin-down can also be
transmitted in quasi-spherical accretion (Shakura et al., 2012). A comparison
of disk and quasi-spherical model predictions for various accreting pulsars is
presented in Malacaria et al. (2020). The two geometries are sketched in
Figure 3, and the implications for Vela X-1 are further discussed in Section
4.2.2.

Figure 3: Disk (a) and quasi-spherical (b) accretion geometries at the outer
edge of the neutron star magnetosphere (represented in green as a dipole),
adapted from Ghosh & Lamb (1978) and Shakura et al. (2012), respectively. In
each of these ideal geometries, different accretion regimes are expected
according to semianalytic derivations, depending on how the plasma couples to
the magnetic field.

###  2.4 Accretion column and X-ray emission

When the plasma has coupled to the magnetic field, it is funnelled to the
magnetic poles of the neutron star where it forms a polar cap or accretion
columns and emits most of the X-ray emission we observe (Davidson & Ostriker,
1973; Lamb et al., 1973; Shapiro & Salpeter, 1975; Arons et al., 1987). No
hysteresis is expected in the sense that a given mass accretion rate will
correspond to a specific X-ray luminosity (averaged over the spin period): The
energy acquired by the flow in its journey from the stellar photosphere to the
neutron star surface has to be released, except when there are significant
outflows from the immediate vicinity of the neutron star. This has not been
observed at this point. When the spin and magnetic axis of the neutron star
are not aligned, which is frequently the case, regular pulses can be observed
at the rotation frequency of the neutron star (see, e.g., the catalog of Liu
et al., 2006, where 66 of 114 HMXBs in the Galaxy have an identified pulse
period). The intensity and shape of the observable pulse profile, that is, the
distinct pattern of the light curve folded with the pulse period, results from
a complex interplay of multiple factors. The emission itself is expected to be
nonuniform and is often highly beamed (Basko & Sunyaev, 1975; Mushtukov et
al., 2018). The emitted radiation is therefore strongly affected by
relativistic light that in the curved space time bends close to the neutron
star (e.g., Meszaros & Riffert, 1988; Riffert & Meszaros, 1988; Kraus et al.,
1995; Odaka et al., 2014; Sotani & Miyamoto, 2018; Falkner, 2018, and
references therein). Close to the neutron star surface, nondipolar components
of the magnetic field might also have a significant effect on the dynamics of
the plasma (Pétri, 2015, 2017; de Lima et al., 2020). Even when the magnetic
field is purely dipolar, there could be many regions on the neutron star
surface to which the plasma is funnelled, and they do not necessarily
correspond to the magnetic poles (Romanova et al., 2004).

The emitted X-ray spectrum of the column is dominated by Compton scattering of
thermal seed photons. Cyclotron emission and scattering processes play an
important role as well (Schwarm et al., 2017b, a, and references therein). A
substantial fraction of accreting X-ray pulsars show relatively broad
absorption line features in their X-ray spectra (see Staubert et al., 2019,
for a recent review). These features result from resonant scattering of X-ray
photons on electrons, whose energies perpendicular to the magnetic field are
quantized into so-called Landau levels. This causes the the plasma to become
optically thick for X-ray photons at these energies and triggers the formation
of cyclotron resonant scattering features (CRSFs) or often simply cyclotron
lines. The spacing of these lines is determined by the energy difference
between adjacent Landau levels, leading to centroid line energies of

$E_{\mathrm{CRSF},n}=\frac{n}{1+z}\frac{\hbar
eB}{m_{e}c}\approx\frac{n}{1+z}11.6\times B_{12}\,\mathrm{keV},$ (7)

where $n$ is number of Landau levels involved in the scattering, $z$ is the
gravitational redshift due to the neutron star mass, $e$ and $m_{e}$ are the
electron charge and mass, $B$ is the magnetic field in the scattering region,
$c$ is the speed of light, and $B_{12}$ is the magnetic field strength in
units of $10^{12}$ Gauss. Determining cyclotron line energies thus gives a
relatively direct measure of the magnetic field strength of the neutron star.
CRSF studies for Vela X-1 are discussed in Section 3.1.7, and the implications
for the neutron star magnetic field are summarized in Section 5.8.

Despite significant efforts and some progress, as further detailed in Sections
4.3 and 4.4, there is currently no self-consistent, physics-based description
that can correctly describe the spectra of accreting X-ray pulsars under all
circumstances.

##  3 Observational diagnostics

In this section, we review the different methods which were used to constrain
the parameters in Vela X-1. We treat the different wavelength domains
separately.

###  3.1 X-ray diagnostics

####  3.1.1 Continuum spectrum

The X-ray luminosity of Vela X-1 is dominated by the X-ray continuum produced
in the accretion column. For a distance of 2 kpc, the median luminosity is
between 4–5$\times 10^{36}$ erg s-1 in the 20–60 keV energy band (Nagase et
al., 1986; Fürst et al., 2010). This approximately translates into a
bolometric luminosity of about $1\times 10^{37}$ erg s-1. However, the Vela
X-1 luminosity varies strongly (e.g., Kreykenbohm et al., 2008) and follows a
log-normal distribution (Fürst et al., 2010). Figure 6 (top) shows the
distribution of the observed hard flux, which is a good tracer of the
intrinsic luminosity. Even on the day-long timescales sampled by Swift/BAT,
large flares of 5–11 times the average luminosity are regularly observed.

The X-ray emission of Vela X-1 outside the eclipse can be phenomenologically
described, as in many other accreting X-ray pulsars, by a power-law shape with
a rollover at higher energies, beyond 15–30 keV (e.g., Nagase et al., 1986;
Kreykenbohm et al., 1999; Orlandini et al., 1998; Odaka et al., 2013).
However, the exact shape of this rollover is unknown. In the simplest
approach, it can be modeled with an exponential cutoff alone, but observations
with a high signal-to-noise ratio at energies beyond 30 keV often require a
more complex shape of the rollover, and different phenomenological models have
been used in the literature (e.g., Orlandini et al., 1998; Odaka et al.,
2013). The continuum is modified at lower energies by variable amounts of
absorption (see Section 3.1.5), a significant iron fluorescence line with an
occasional iron edge (Nagase et al., 1986) and a frequent thermal soft excess
below $\sim$3 keV (Sato et al., 1986; Pan et al., 1994; Haberl, 1994), as well
as further fluorescence lines from elements such as silicon, magnesium, and
neon (e.g., Sako et al., 1999; Goldstein et al., 2004). At higher energies,
the continuum is modified by CRSFs between 25 keV and 30 keV and at $\sim$55
keV, as further detailed in Sect. 3.1.7.

The overall continuum slope and the CRSFs vary strongly with the phase of the
X-ray pulse (e.g., Kreykenbohm et al., 2002; La Barbera et al., 2003), but
within the scope of this paper, we do not discuss this further.

During the eclipse, the residual X-ray emission shows no pulsations and is
dominated by line emission in addition to a flat continuum Nagase et al.
(1994); Schulz et al. (2002).

In general, the X-ray flux observed from Vela X-1 varies significantly in a
largely random pattern around the mean value when the clearly visible X-ray
eclipse is excluded. These variations are driven by variations in the
intrinsic X-ray emission and by strong variations in the absorbing material
between observers and the neutron star, especially for instruments in the soft
X-ray band.

Figure 4: Variations in derived absorbing column density ($N_{\mathrm{H}}$ /
cm2) throughout the binary orbit of Vela X-1 as determined by various
observations and different satellites. See the text for caveats when absolute
values in $N_{\mathrm{H}}$ are compared. The shaded gray areas mark the
eclipse (dark gray) and the ingress and egress (light gray) as determined by
Sato et al. (1986). Binary phase values have been corrected for the
differences in phase-zero definition (using $T_{\mathrm{ecl}}$) and orbital
period relative to Kreykenbohm et al. (2008). The uncertainty in this
conversion is 1–4$\times 10^{-3}$ in phase, i.e., usually within the symbol
size. Care has been taken to separate measurements taken during different
binary orbits by the same instrument and team. Where possible, this is marked
by different symbols. The data have been taken from Sato et al. (1986, Fig. 5
(Tenma)), Haberl & White (1990, Fig. 2 (EXOSAT)), Lewis et al. (1992, Fig. 2
(Ginga)), Pan et al. (1994, Fig. 3 (TTM/HEXE)), Watanabe et al. (2006, Tab. 2
(Chandra 2001)), Martínez-Núñez et al. (2014, Fig. 7 (XMM-Newton)) and Liao et
al. (2020, Tab. 1 (Chandra 2017)). Panel a) gives an overview of the overall
distribution, demonstrating the large scatter of results at intermediate
orbital phases and that the larger structures driving $N_{\mathrm{H}}$ are not
stable from one orbit to the next (see also Malacaria et al., 2016). Panels b)
to e) show subsets of the data in more detail to better visualize short-term
variations. In all these panels, the X-axis covers the same range, and the
Y-axis covers the same factor between minimum and maximum.

####  3.1.2 Eclipse timing

Despite the sparse observations of the time, it was soon realized in the early
years after the detection that Vela X-1 was a variable source (Giacconi et
al., 1972, and references therein). Based on long observations with the OSO-7
UCSD X-ray telescope, Ulmer et al. (1972) found evidence for periodic flux
variations with phases of zero flux at a period of $8.7\pm 0.2$ d and proposed
that the X-ray source might be a member of an eclipsing binary system. This
result was confirmed, and the periodicity was refined to $8.95\pm 0.02$ d by
Forman et al. (1973), based on Uhuru data. In subsequent years, the eclipse
timing and derived orbital parameters were refined by various authors, for
example, by Ögelman et al. (1977), using 20 days of COS-B data, or by Watson &
Griffiths (1977), using Ariel V data covering 17 binary cycles, see Table 3.
Most recently, Falanga et al. (2015) combined eclipse times derived from
INTEGRAL observations with the earliest data to investigate a possible
evolution of the orbital parameters (Sect. 5.2).

####  3.1.3 Variations in the X-ray absorption

Early observations of Vela X-1 (Eadie et al., 1975; Watson & Griffiths, 1977;
Becker et al., 1978) already found indications for significant and variable
absorption in the system, together with brightness variations. Because of
instrumental limitations, the brightness variations could not always be
clearly ascribed to variations in absorbing material or in the intrinsic X-ray
flux. Later observations with improved X-ray continuum spectra allowed better
measurements of the actual absorption and found strong variability on many
timescales. In Tenma data, Ohashi et al. (1984) found large variations of the
absorbing column density $N_{\mathrm{H}}$ that ranged from 2 to $50\times
10^{22}$ cm-2 on a timescale of about an hour.

Over the years, very many observations have among other properties followed
the temporal evolution of the absorbing column density in the Vela X-1 system;
see Figure 4 and Table LABEL:tab:observations_xray. When they are presented in
the same figure as a function of the orbital phase, different sets of
$N_{\mathrm{H}}$ measures can be misleading because they often do not belong
to the same orbit and because the column density at a given orbital phase is
not stable from one orbit to the next. The overall behavior can be summarized
by the following statements. First, from orbital phase
$\phi_{\mathrm{orb}}\sim 0.1$ (eclipse egress) to $\sim 0.25$,
$N_{\mathrm{H}}$ generally decreases strongly, often by more than one order of
magnitude. This is usually explained by suggesting that the X-ray source
shines through the extended atmosphere of its companion. Occasionally, shorter
deviations from this trend have been found (Nagase et al., 1986; Martínez-
Núñez et al., 2014). These are most probably caused by absorbing material
closer to the X-ray source. Second, between $\phi_{\mathrm{orb}}\sim 0.25$ and
$\sim 0.6,$ strong variations of $N_{\mathrm{H}}$ on timescales of hours to
days are observed. At any given phase, high and low absorbing column densities
have been found in different binary orbits, sometimes using the same
instrument and analysis method (e.g., Nagase et al., 1986). This also applies
to $\phi_{\mathrm{orb}}\sim 0.5$, where Haberl & White (1990) and Watanabe et
al. (2006) found high $N_{\mathrm{H}}$ values, which the latter modeled as
broad cold cloud, but where Nagase et al. (1986) found an intermediate
minimum, as shown in Figure 4. In a systematic study of MAXI data, Malacaria
et al. (2016) found that $\sim 15$% of all orbital light curves showed a dip
around $\phi_{\mathrm{orb}}\sim 0.5$, which was best explained by Thomson
scattering in an extended and ionized accretion wake. Third, for late orbital
phases ($\phi_{\mathrm{orb}}>0.6$) up to the eclipse ingress, the observed
absorbing column densities always remain high, but are still significantly
variable by a factor of a few.

The following caveats apply when results from different publications on the
absolute values of $N_{\mathrm{H}}$ are applied. First, spectral modeling
results are usually somewhat degenerate for the parameters defining the
spectral slope (e.g., a power-law index) versus the amount of absorbing
material. The impact of this depends on the energy range that is covered, on
the spectral response of each detector, and on the assumptions about permitted
models in the fit procedure. For historical missions, it is frequently
impossible to repeat the data reduction and analysis, even if time and
knowledge were available for lack of suitable data in archives. Second and
similarly, the inclusion or exclusion of an assumed soft excess can have a
strong effect on the derived $N_{\mathrm{H}}$. Third, while many publications
have used a single, unique absorbing column to obtain their results, in some
particularly detailed studies, multiple absorption components were required to
fit all of the available spectra (Haberl & White, 1990; Martínez-Núñez et al.,
2014). In this case, one absorbing component may have significantly higher
$N_{\mathrm{H}}$ values, but this only affects a fraction of the total
observed radiation. Fourth, as detailed in Martínez-Núñez et al. (2017, Sect.
4.4.1,) various different absorption models have been published in the
literature and used over the years by different authors. These models rely on
different assumptions for the interstellar medium cross-sections and
abundances of the different elements, which can lead to differences of more
than 20% in the obtained $N_{\mathrm{H}}$ .

Figure 5: Effect of variable absorption on the overall shape of the spectrum
of Vela X-1 in the 0.5–10 keV range. Following Martínez-Núñez et al. (2014),
we assume a spectral model consisting of three power-law components with the
same slope, but different normalizations, each absorbed by a different
absorption component. In accordance with trends seen by Martínez-Núñez et al.
(2014), we then vary the total amount of absorption that the second of these
power-law components experiences. The dot-dashed gray line shows the first
power-law component. The dashed lines show the second power-law component, and
different colors represent different absorption strengths, as indicated in the
plot. The dot-dot-dashed line shows the third power-law component. The sum of
the three components, i.e., the total spectrum, is shown as solid lines, and
colors again represent different absorption levels for the second component.

We show the effect of variable absorption on the spectrum of Vela X-1 in Fig.
5. In particular, we demonstrate that even a combination of relatively simple
continuum models can lead to complex variability in spectral shape. To do so,
we used the model that was introduced to describe the 0.5–10 keV spectrum of
Vela X-1 by (Martínez-Núñez et al., 2014), namely a sum of three power-law
components that have different normalizations, but the same spectral index.
Each of the power laws is absorbed by a distinct absorption component. A
possible physical interpretation of this model is a complex partial coverer,
possibly with an additional soft excess contribution (e.g., Grinberg et al.,
2020). Martínez-Núñez et al. (2014) found that the main driver for the
variable absorbed spectral shape of Vela X-1 during a flaring episode was the
absorption of the second, most prominent power-law component. Fig. 5 is
calculated for typical values of the power-law index and relative power-law
component contributions. It clearly shows how changes in absorption in this
component can lead to strong nonlinear variability in the overall spectral
shape.

Grinberg et al. (2017) have explored a different temporal and spatial scale.
In order to test the effect of clumps in the wind on short-term variations in
the absorbing column density, they compared data from a Chandra observation of
Vela X-1 at orbital phases $\sim$0.21 to $\sim$0.25 with a toy model that
consisted of a constant X-ray source and a clumpy wind, based on the
simulations by Sundqvist et al. (2018). They found that undisturbed
inhomogeneities in the model wind could not account for the observed enhanced
absorption (see also El Mellah et al., 2020b, for an extended model),
indicating that more complex factors such as shock-clump interaction, the
possibility of transient disk structures, or the effect of the compact object
and its X-ray emission on the wind itself must be included.

####  3.1.4 X-ray fluorescence lines and radiative recombination continua

Narrow spectral features such as lines, edges, and radiative recombination
continua (RRCs) are imprinted onto the X-ray continuum by the material
surrounding the neutron star. They can in turn be used to constrain the
properties of the surrounding (wind) plasma, such as ionization and density
structure.

Especially the iron lines in the spectral region of 6–7 keV, the so-called Fe
complex, constitute a fundamental tool for probing the physical conditions of
the material in the close vicinity of the X-ray source (George & Fabian,
1991). In the spectra of accreting X-ray binaries, these lines are prominent
because of their intrinsic X-ray brightness and ubiquitous stellar material.
In Vela X-1, iron lines were first reported by Becker et al. (1978) and were
studied by several teams since then, including studies with Chandra and XMM-
Newton (see, e.g., Goldstein et al., 2004; Watanabe et al., 2006; Torrejón et
al., 2010; Giménez-García et al., 2016; Grinberg et al., 2017). Fe K$\alpha$
and K$\beta$ have been detected with multiple instruments and are narrow and
usually unresolved, even with the highest currently available resolution using
the Chandra transmission gratings (e.g., Watanabe et al., 2006; Goldstein et
al., 2004; Grinberg et al., 2017; Amato et al., 2021). The line energy is
consistent with either neutral or lowly ionized ions, below $\sim$Fe xii.
Three reprocessing sites have been suggested based on observational trends,
such as deviations from curve of growth for spherical geometry. All three
likely contribute to some extent: the stellar wind as a whole, the companion
atmosphere, and cold material in the vicinity of the neutron star, for
example, in the accretion wake (e.g., Sato et al., 1986; Watanabe et al.,
2006).

A plethora of other narrow features, in particular, multiple fluorescent lines
from different instruments, have already been detected with ASCA (Sako et al.,
1999) and were resolved in detail using grating instruments on board Chandra
and XMM-Newton. Detections include prominent lines in both emission and
absorption, edges, and RRCs of sulfur, silicon, magnesium, neon, and oxygen at
different ionization levels (Schulz et al., 2002; Goldstein et al., 2004;
Watanabe et al., 2006; Grinberg et al., 2017; Lomaeva et al., 2020; Amato et
al., 2021). The strength of different features can change with orbital phase
(Goldstein et al., 2004; Amato et al., 2021), with the observed amount of
absorption (Grinberg et al., 2017), and with changing irradiation through the
neutron star, for example, during a flare (Lomaeva et al., 2020). While the
total normalization of the spectrum changes dramatically between eclipse and
$\phi_{\mathrm{orb}}\approx 0.5$, striking similarities between the relative
strength of fluorescence lines have been seen in two Chandra/HETG observations
at these orbital phases covering the same orbit (Sako et al., 2003; Watanabe
et al., 2006). This behavior allows us to rule out not only a smooth wind
model, but also several simple clump distributions, such as a population of
identical clumps (Sako et al., 2003).

Multiple modeling approaches have been employed to describe the observed
spectra. Sako et al. (1999) determined the differential emission measure
distribution at $\phi_{\mathrm{orb}}\approx 0$ from modeling contributions of
individual ions and from comparison to theoretical emission measure
distribution for a neutron star immersed in the stellar wind that are
parameterized by wind parameters (mass loss and velocity) alone. They can
obtain a good description of the observed spectrum and constrain the mass
loss. However, their model predicts a lower ionization at
$\phi_{\mathrm{orb}}\approx 0.5$ that has subsequently not been observed; the
discrepancy is due to the assumption of spherical symmetry, which is clearly
not present in Vela X-1 (Sako et al., 2003). Lomaeva et al. (2020) and Amato
et al. (2021) attempted to describe XMM-Newton/pn and Chandra/HETG high-
resolution spectra, respectively, using photoionization models and comparing
different photoionization models, namely XSTAR, CLOUDY, and SPEX/pion. While
the fits are satisfactory, the residuals clearly point toward the existence of
more than one component, in agreement with the assumption of a multiphase,
clumpy medium. Their results clarify that a more sophisticated treatment of
the system is necessary. Watanabe et al. (2006) performed a quantitative
modeling and spectral analysis using Chandra HETG during three orbital phase
ranges within a single binary orbit in 2001: $\phi_{\mathrm{orb}}\approx
0.25$, $\phi_{\mathrm{orb}}\approx 0.5$, and $\phi_{\mathrm{orb}}\approx 0$
(during eclipse). Their aim was to constrain the ionization parameters and the
geometrical distribution of the material in the X-ray irradiated stellar wind.
They computed the ionization structure of a wind whose density and temperature
distribution is given by a simple empirical model. With a Monte Carlo
simulation of X-ray photons propagating through the wind, they computed
synthetic X-ray spectra and established a self-consistent radiative-
hydrodynamics picture of the system, but assuming a smooth and not a clumpy
wind. For a more thorough discussion of simulations of the plasma properties
in the Vela X-1 system of these and other authors, we refer to Sec. 4.1.2.

Figure 6: Brightness distribution of Vela X-1 throughout its orbit as measured
with Swift/BAT (15–50 keV, top panel) and MAXI (2–20 keV, bottom panel). In
comparison, 1 Crab corresponds to values $\sim$0.22 for Swift/BAT and $\sim$4
for MAXI. We used the full history for each instrument, binned to a resolution
of 1 d. The mean orbital flux profile is superimposed in black. The
probability for a measurement to fall into the respective histogram bin is
color-coded according to the color scales on the right. The profile is
repeated once for clarity. For details, see text.

####  3.1.5 Overall flux and continuum variations

In addition to the effect of absorption variations, any longer observation of
Vela X-1 has found significant variations of the source flux on many different
timescales, hours or tens of minutes (see, e.g., Haberl & White, 1990;
Kreykenbohm et al., 2008), down to changes from one X-ray pulse to the next
(e.g., Inoue et al., 1984; Börner et al., 1987). There are many examples of
flares, with timescales ranging from a single pulse to several hours (e.g.,
Watson & Griffiths, 1977; Börner et al., 1987; Odaka et al., 2013; Martínez-
Núñez et al., 2014).

Off-states outside eclipse with no or very little detectable X-ray emission
have been reported by various authors (Inoue et al., 1984; Lapshov et al.,
1992; Kreykenbohm et al., 1999, 2008; Doroshenko et al., 2011; Sidoli et al.,
2015). These states have been explained by different combinations of possible
low-density zones or accretion that is choked or gated by magnetospheric
effects (see Section 4.5). One caveat for a comparison of different
publications is that the effective sensitivities of the instruments used and
the definitions of “off” vary significantly.

Fürst et al. (2010) used data taken in the hard X-ray band (20–60 keV) with
INTEGRAL for a systematic study of the brightness distribution in an energy
range that is not affected by photoabsorption below $\sim 3\times
10^{23}\,\mathrm{cm}^{-2}$. They found a log-normal distribution of the X-ray
brightness and that the power spectral density of the light curve follows a
red-noise power law. They suggested that a mixture of a clumpy wind, shocks,
and turbulence can explain the measured distribution.

The strong variability of Vela X-1 can be seen in the orbital-phase-resolved
histograms shown in Fig. 6. While the average orbital profile is stable and
well defined, large variations around that mean value are also present. This
is indicative of large flares and low states. Notably, the hard Swift/BAT
profile is much more symmetric between eclipses than the softer one as
measured by MAXI, which indicates a strong rise in absorption (on average) at
later orbital phases. The difference between the mean value and the median
that is apparent through the color-coding is expected for a log-normal flux
distribution.

Table 1: Overview of off-states or low-flux observations of Vela X-1 reported in the literature. References | Date | Instrument | Energy range | Off-state luminosity
---|---|---|---|---
Kallman & White (1982), Hayakawa (1984) | 9 May 1979 | Einstein | 2-10 keV | _not quantified_
Inoue et al. (1984), Hayakawa (1984) | 12 Mar 1983 | Tenma | 3-9 keV | _$ <$10% of preceding_
Lapshov et al. (1992) | 9 Jan 1991 |  Granat/WATCH | 8-15 keV | _not quantified_
Kreykenbohm et al. (1999) | 23 Feb 1996 |  RXTE/PCA | 3-30 keV | _$ <$15% of normal_
Kretschmar et al. (2000) | 22 Jan 1998 |  RXTE/PCA | 2-60 keV | _not quantified_
Kreykenbohm et al. (2008) | 8 Dec 2003 |  INTEGRAL/ISGRI | 20-40 keV | _not quantified_
Doroshenko et al. (2011) | 17–18 Jun 2008 |  Suzaku/XIS | 0.4–70 keV |  $\sim 2.4\times 10^{35}$ erg s-1
Sidoli et al. (2015) | 2002–2012 |  INTEGRALISGRI | 22-50 keV |  $\lesssim 3\times 10^{35}$ erg s-1

Hand in hand with the flux variations, a hardening of the spectrum has been
observed in Vela X-1, that is, there is an anticorrelation between the photon
index $\Gamma$ and the observed flux (Fürst et al., 2014). This behavior
agrees with the expected subcritical accretion regime of Vela X-1 and is also
seen in sources of similar luminosity (Klochkov et al., 2011).

####  3.1.6 X-ray polarization

For lack of sensitive X-ray polarimeters in space, no observations of X-ray
polarization parameters in Vela X-1 or similar systems from X-ray satellites
have so far been made. A notable recent exception was the balloon-borne hard
X-ray calorimeter, X-Calibur, which observed Vela X-1. Because the flight
duration was much shorter than anticipated, however, no constraining data
could be obtained (Abarr et al., 2020).

Clear polarization signatures are expected from the physics of the radiation
transfer in the emission region as well as in the system as a whole, however.
Photon-scattering opacities within the accretion column will depend strongly
on energy, on the direction of propagation, and and on polarization, as
discussed by Meszaros et al. (1980) and in various subsequent publications.
Thus, a significant variation of polarization parameters is expected as
function of pulse phase and energy for the X-ray emission originating at the
neutron star.

Recently, Caiazzo & Heyl (2021) have presented a detailed calculation of the
intrinsic polarization of the X-ray point source for accreting pulsars in
X-ray binaries. They took the structure and dynamics of the accretion region
into account and included relativistic beaming, gravitational lensing, and
quantum electrodynamics in their treatment of the propagation of the radiation
toward the observer.

Complementary to this, Kallman et al. (2015) computed scattering of the
emitted X-rays within the dense winds of HMXBs. They showed that it may lead
to a variation in polarization with orbital phase, energy, and system
geometry, especially if large-scale structures are present in the wind. They
described how they obtained the three Stokes parameters of linearly polarized
light from this scattering. However, their calculations assumed an
unpolarized, isotropic X-ray emitter, which is likely not the case for
accreting X-ray pulsars. They found fractional polarization values of up to
$\sim$10% for a spherically symmetrical wind and $>$20% in their more
realistic hydrodynamical model with a focused wind.

####  3.1.7 Cyclotron resonance scattering features

Kendziorra et al. (1992) first reported CRSFs (Section 2.4) in Vela X-1 at 54
keV and possibly 27 keV, based on pulse-phase-resolved spectra obtained from
the High-Energy X-ray Experiment (HEXE) on board the Mir space station. This
made Vela X-1 one of the early identified sources with such line features
(Staubert et al., 2019). Evidence for these features was also given by
Makishima & Mihara (1992). While the detection of line features was not
challenged, it was debated in the literature for at least a decade whether the
feature between 25 keV and 30 keV was the fundamental cyclotron line (e.g.,
Kretschmar et al., 1997; Kreykenbohm et al., 1999) or rather the harmonic
(Orlandini et al., 1998) because the lower-energy feature was not always
apparent, especially in phase-averaged spectra. Kreykenbohm et al. (2002)
found in a deep pulse-phase-resolved analysis of spectra taken with RXTE a
fundamental line feature at $23.3^{+1.3}_{-0.6}$ keV during the pulse maxima
that was much less significant during the pulse minima. All line parameters
showed clear pulse-phase dependence.

Observations of Vela X-1 with NuSTAR (Fürst et al., 2014) clearly detected the
fundamental line at $\sim$25 keV together with the more prominent harmonic at
$\sim$55 keV. Fig. 7 shows one of the spectra presented by Fürst et al.
(2014), which is a typical out-of-eclipse Vela X-1 spectrum. The data
presented here were extracted with NUSTARDAS v2.0.0 and CALDB v20201217. The
model is based on the best-fit model presented by Fürst et al. (2014), which
uses an absorbed Fermi-Dirac cutoff continuum spectrum, modified by two
multiplicative lines with Gaussian optical depth profiles to describe the
CRSFs. In addition, two Gaussian emission lines are used to model the Fe
K$\alpha$ and K$\beta$ lines, and a Gaussian absorption line was added to
describe the 10 keV feature (La Barbera et al., 2003). In time-resolved
spectroscopy, Fürst et al. (2014) found that the strengths of the two CRSFs
appeared to be anticorrelated, which could be explained by photon spawning
(Schönherr et al., 2007).

Earlier studies found no significant variation of the CRSF energy with flux.
Fürst et al. (2014) found a positive correlation of the harmonic (higher-
energy) line with flux, as expected for a source in the subcritical accretion
regime as defined by Becker et al. (2012). La Parola et al. (2016) undertook a
long-term study of Swift/BAT data of Vela X-1. While the fundamental line was
not evident in their 15–150 keV spectra, the first harmonic of the CRSF was
clearly detected between $\sim$53 and $\sim$58 keV, and with an apparent
positive correlation of the line energy with luminosity, flattening above
$L_{\textrm{1--150\leavevmode\nobreak\ keV}}\approx 5\times 10^{36}$ erg s-1.
Recently, Ji et al. (2019) reported a secular decrease of the energy of the
harmonic line, again based on Swift/BAT data . If confirmed, this would
indicate a slow change in the magnetic field configuration close to the polar
caps.

One important caveat when results from different studies are compared is that
different choices of mostly phenomenological spectral models for the continuum
and for the cyclotron line parameters themselves lead to fit parameters that
are intrinsically not directly comparable. Even the line centroids determining
$E_{\text{CRSF},n}$ (Eq. 7) can be systematically different between different
line models by a few keV for the same fitted spectrum, as explained in detail
in Staubert et al. (2019, Sect. 3).

Figure 7: a) X-ray spectrum of Vela X-1 as measured with NuSTAR (ObsID
30002007003). Data from the two instruments of NuSTAR, FPMA and FPMB, are
shown in red and blue, respectively. The best-fit model as presented by Fürst
et al. (2014) is shown in black. The dashed green lines indicates the
continuum model without the CRSF components. b) Residuals in terms of
$\chi^{2}$ between the data and the best-fit model. c) Residuals between the
data and the model without the CRSF components. The harmonic line around 55
keV is much more prominent than the fundamental line around 25 keV.

####  3.1.8 Pulse period variations

Figure 8: Evolution of the pulse period of Vela X-1 as observed by diverse
instruments (gray squares, see Table LABEL:tab:pulse_periods for details), by
CGRO/BATSE (blue circles) or by Fermi/GBM (red circles). BATSE and GBM period
determinations are derived over two and three intervals of the orbit,
respectively. Note the scale. The longest reported pulse period differs by
$\sim$0.3% from the shortest known pulse.

The long-term pulse period evolution of Vela X-1 is usually described as a
random walk. In an in-depth study of HEAO-1 data, Boynton et al. (1984) found
that the power density spectrum covering timescales from 0.25 to 2600 days can
be described by assuming white noise in angular acceleration, or equivalently,
a random walk in pulse frequency. Short-term angular accelerations were found
to be as large as $\dot{\Omega}/\Omega=(5.8\pm 1.4)\times
10^{-3}\textrm{yr}^{-1}$. According to their analysis, apparent changes in
secular trends and frequency variations on shorter timescales of at least days
are all consistent with the same random noise process in angular acceleration.
This result was reinforced and refined in subsequent papers (Boynton et al.,
1986; Deeter et al., 1989) and also by a similar study based on the Hakucho
and Tenma satellites (Deeter et al., 1987). de Kool & Anzer (1993) extended
these studies by combining data from Nagase (1989) and Raubenheimer & Ögelman
(1990), again confirming that the pulse period behavior is very well fit by a
random walk. Boynton et al. (1984) also noted that the short-term angular
accelerations reported in their study were much larger than expected for
accretion from a uniform wind, assuming the wind speeds given in Dupree et al.
(1980), and that the flow may form a ring or disk around the neutron star.

####  3.1.9 Pulse profiles

The first observations discovering the periodic pulsations of Vela X-1 with
the SAS-3 X-ray observatory have found a complex pulse profile structure with
five visible peaks at energies below 6 keV and two broad peaks at the higher
energies above $\sim$10 keV, which are both somewhat asymmetric, with a
sharper trailing than leading edge (Rappaport & McClintock, 1975; McClintock
et al., 1976). The overall structure of the pulse profile has been confirmed
in many subsequent studies of the system (e.g., Nagase et al., 1983;
Raubenheimer, 1990; Kreykenbohm et al., 1999; La Barbera et al., 2003), across
significant variations in X-ray brightness. It can be considered a fingerprint
of the X-ray pulsar, like in many other systems. In other words, while some
variations in the observed profile are visible when different observations are
compared, especially in the peak fluxes of different pulse components, the
overall shape and thus the underlying emission geometry is evidently very
stable. This is also demonstrated in Figure 9. The main visible differences at
the lower energies can be caused by the sometimes very strong absorption and
scattering (Section 3.1.3), which can smear out especially the low-energy
profile (Nagase et al., 1983). For detailed pulse profiles across a wide
energy band, see especially Raubenheimer (1990, Fig. 1) and La Barbera et al.
(2003, Figs. 1–4).

Figure 9: Comparison of the complex low-energy pulse profile in the roughly 2
to 6 keV energy range from observations taken decades apart with three
different X-ray missions, SAS3, EXOSAT, and BeppoSAX, respectively. The data
points have been taken from McClintock et al. (1976); Raubenheimer (1990); and
La Barbera et al. (2003), but shifted in phase and scaled arbitrarily for a
visual match.

When different observations are compared, we recall that like in many other
accreting X-ray pulsars, there are clear variations between individual pulses,
as demonstrated in the first detailed study of the hard X-ray pulsations
(Staubert et al., 1980) and later studies (Kendziorra et al., 1989; Kretschmar
et al., 2014). These variations may affect the mean pulse profile of
observations that cover only a limited number of pulse period intervals.

In earlier publications, pulsations were found to cease during off-states
(Inoue et al., 1984; Kreykenbohm et al., 1999, 2008). This behavior might be
interpreted as the onset of the propeller regime (Section 2.3). In contrast,
Doroshenko et al. (2011) analyzing three off-states observed with Suzaku, did
observe remaining pulsations, but with significant changes in the pulse
profile. They interpreted this behavior as gated accretion, with some
remaining accretion.

A very peculiar observation was reported by Kretschmar et al. (2000): For a
duration of several hours, the source flux diminished, but was still
significantly above the background level, while pulsations ceased practically
completely. After this low-state phase, there was a high state with strong,
flaring pulsations, but a very similar spectrum to the preceding nonpulsating
state. At the time, an explanation by a very thick clump in the wind was put
forward, but in the light of newer results on wind properties and realistic
clump sizes (see Sect. 2), this seems very unlikely.

Liao et al. (2020) reported another case of a low-flux state with ceasing
pulsations. Based on spectral evidence, this was interpreted as formation of a
(temporary) accretion disk in this case. We also note the explanation put
forward by King & Lasota (2020) for suppressed pulsations in ultraluminous
X-ray sources by scattering in a thick beaming funnel. This role would be
played in Vela X-1 by the transient disk-like structure.

###  3.2 Optical, UV, and IR diagnostics

####  3.2.1 Photometric light curve

The study of photometric data, and in particular, the shape and amplitude of
the optical light curve, is a tool that beyond determining the photometric
period can be used through comparison with models to infer information on
system parameters such as the size of the components relative to the orbit or
the mass ratio. For a review, see Wilson (1994a); a visualization of a model
light curve for Vela X-1 is included in Wilson (1994b). Several studies of
this type have been performed in the 1970s and 1980s. All of them reported
variability on the photometric light curve from one cycle to another. At later
times, system parameters were instead obtained from X-ray studies, which
allowed a more accurate determination of the system parameters.

The first photometric study of HD 77581 was carried out by Vidal et al. (1973)
in the V band over a time span of 16 days. They observed four minima in the
light curve that were consistent with the reported period of $8.95\pm 0.02$
days by Forman et al. (1973). They reported a small amplitude of the light
curve of $\sim 0.14$ mag and found some evidence of variability with the
period, particularly at the primary minimum. This behavior was attributed by
Vidal et al. (1973) to the important role played by tidal effects.

The same periodicity was confirmed by Jones & Liller (1973) using UBV
photometric observations of the companion over 27 consecutive nights,
including two maxima and two minima. They reported an unexpected
nonrepeatability of the optical light curve from cycle to cycle with color-
free erratic changes. They explained this as due to complex changes in the
upper layers of the star that are comparable to the gravity-darkening effect.

Petro & Hiltner (1974) combined the observations of Vidal et al. (1973) and
Jones & Liller (1973) with three new campaigns to obtain a more precise
determination of the photometric period of 8.972(1) days. They noted
systematic deviations between the data and the best-fitting circular-orbit
light curve, which might be explained by an elliptical orbit.

Vidal (1974) presented a more extensive set of observations in V band with 126
observations in total, including those reported in Vidal et al. (1973), which
were obtained during the period of 22 January to 19 June 1973. They reported
that the most interesting feature of the light curve was the high scatter of
the data throughout the whole period and claimed that this variability was
intrinsic to the source. They observed that the rising and descending branches
of the light curve during primary minimum were different and interpreted the
shape of the light curve as due to different aspects of a tidally distorted
companion.

Zuiderwijk et al. (1977b) presented a four-color photometric study of HD 77581
in the Strömgren uvby system. They observed the source on 46 nights in three
campaigns in 1975. The magnitudes are listed in tabular form in Zuiderwijk et
al. (1977a). These observations also showed changes in the light curve from
one orbit to another. They found a regular wave-like variation in the color
index c1 and possibly in b-y in phase with the light curves and twice noted
the disappearance of a maximum during their observations, as also reported by
Jones & Liller (1973). The observed variability in the light curve and the
colors were interpreted with a model of a tidally distorted rotating primary.

van Genderen (1981) carried out a VBLUW photometric study in 1976 and 1977.
They obtained a photometric period of $8.9615\pm 0.0025$ days, which is closer
to the contemporary measured period of the X-ray light curve. They reported
the peculiarities of the mean light curve, which had a total amplitude of
$\sim$ 0.10 mag, with a difference in height between the maximum and the
minimum of $\sim$ 0.04 mag. The minimum seemed to be delayed by $\sim$0.1 in
phase with respect to the center of the X-ray eclipse. The high peaked maximum
varied in height and possibly in phase with a timescale of about two years.
They concluded that the observed intrinsic variability of the light curve was
correlated with the appearance and disappearance of hotter or cooler areas in
the stellar atmosphere.

Khruzina & Cherepashchuk (1982) performed an analysis of all published
photometry data and revealed a regular long period of 93.3 days of the optical
light curve (B band). This period was later confirmed by Cherepashchuk (1982)
using independent new UBV photometric data obtained at Siding Spring
Observatory in February-June 1980. Khruzina & Cherepashchuk (1982) and
Cherepashchuk (1982) associated this long-term period with forced precession
of the rotation axis of the optical companion. This induces changes in the
eclipsing gas streams and in the accretion structure. While super-orbital
periods have been detected for other HMXBs (Corbet & Krimm, 2013), no such
periodicity has been reported for Vela X-1.

Figure 10: V-band magnitudes of HD 77581 as compiled by Tjemkes et al. (1986,
Fig. 8), demonstrating the intrinsic scatter, which is much larger than
measurement uncertainties. For visual comparison, we plot the mean orbital
profile in X-rays as observed by the MAXI instrument from Figure 6 in the
background.

.

Tjemkes et al. (1986) developed a simple geometric model to analyze the
optical light curves of several X-ray binaries and compared their predictions
with real data. In the case of HD 77581, they used all published data in
addition to data obtained in 1979 and 1980 to compute an average visual light
curve. Ellipsoidal light variations were observed in the light curve, which is
indicative of a donor star close to filling its Roche lobe. The light curve
had two brightness maxima and two minima, and a photometric period of 8.965
$\pm$ 0.001 was reported. However, the light curve presented some
particularities: 1) the two maxima are different in shape and brightness
level, 2) the minimum is somewhat asymmetric and is shifted with respect to
the mid-time of the X-ray eclipse by 0.05, and 3) the scatter around the
average light curve is $\sim$ 0.02, which is much larger than the accuracy of
the photometric data. These variations occurred on timescales from hours to
many days. They carried out a search of the long-term period reported by
Khruzina & Cherepashchuk (1982) and Cherepashchuk (1982) and detected some
indications in favor of a minimum variance near 93.2 d, but no significant
evidence to support the existence of this super-orbital period. Overall,
Tjemkes et al. (1986) found serious discrepancies between their model and the
data for the minima and concluded that the model did not describe the region
near the inner Lagrangian point well. They suggested that basic model
assumptions, in particular the instantaneous adjustment of the primary shape
to the equipotential surface, are invalid for an eccentric orbit.

Wilson & Terrell (1994) fit optical light curves, optical velocities, and
X-ray pulse arrival times simultaneously to reproduce the observed X-ray
eclipse duration. They assumed a binary star model on equipotentials with an
eccentric orbit and adjusted the varying tidal effect. The companion star
rotated at a constant angular rate. They concluded that the optical star
rotates subsynchronously and presents phase shifts that the authors
interpreted as tidal lags. Both of these properties have major implications
for understanding the scenario of the system evolution. According to Wilson &
Terrell (1994), the system could be on the verge of a common evolution
scenario.

####  3.2.2 Radial velocity curve

In an eclipsing X-ray binary system, the procedure for determining the neutron
star mass is based on the radial velocity (RV, hereafter) curve of the
companion star, as determined from shifts in the positions of various observed
spectral lines, knowing the orbits of the two components, and knowing or
assuming the inclination of the system. This procedure is valid under the
assumption that RV curve variations truly represent the orbital motion of the
star. In the case of HD 77581, this assumption is highly questionable for two
reasons: the prominent line profile variability (van Paradijs et al., 1977b;
van Kerkwijk et al., 1995; Barziv et al., 2001; Quaintrell et al., 2003), and
the asymmetries over the orbital cycle of the RV curve (van Paradijs et al.,
1977b; van Kerkwijk et al., 1995; Barziv et al., 2001).

The first RV study of HD 77581 was performed by Hiltner et al. (1972) using a
Coudé spectrograph with a total of 85 spectrograms. They found variability on
the source, and the velocity curve was not sufficiently defined to warrant a
detailed analysis. They concluded that the system was a spectroscopic binary
with a provisional period of 7d and a minimum mass of the companion near 1.4
$M_{\sun}$.

In 1974, several publications (Wickramasinghe et al., 1974; Mikkelsen &
Wallerstein, 1974; Zuiderwijk et al., 1974; Hutchings, 1974) carried out an
analysis of the RV curve. No consistent orbital elements could be derived,
however, leading to estimates of the mass of the compact object that were too
large to be consistent with those of a neutron star.

van Paradijs et al. (1976) performed the first relatively accurate
determination of masses for both partners using RV measurements. They used 26
Coudé spectrograms to determine the total mass of the system as well as the
masses of each component. In contrast to previous work, they excluded the
lines of hydrogen because they are the most affected by gas motions in the
system. Using He I lines and lines of heavier ions, they derived mean values
of the RV for all the lines and for the He I and heavier-ion lines
independently. They obtained a semiamplitude of the RV curve of
$K_{\textrm{opt}}=20\pm 1\,\textrm{km s${}^{-1}$}$ and derived a total mass of
the system of $\sim 21.6M_{\sun}$, with a mass of the compact object of $\sim
1.6M_{\sun}$, and $\sim 20M_{\sun}$ for the companion (see Table 4).

van Paradijs et al. (1977b) carried out a more detailed optical RV study to
improve the accuracy of previous works and to investigate the possible
presence of the distortion effect predicted in van Paradijs et al. (1977a).
They found short-term autocorrelated nonorbital variations, but because so
many factors were involved, they were unable to estimate the amplitude of the
individual nonorbital contributions. They derived a more accurate
$K_{\textrm{opt}}=21.75\pm 1.15\,\textrm{km s${}^{-1}$}$ than in van Paradijs
et al. (1976) and claimed that the mass of the primary was rather low for its
spectral type and luminosity class (see Table 4), a similar type of
undermassiveness as was found for SK 160 in SMC X-1.

Nearly 20 years later, van Kerkwijk et al. (1995) performed a RV curve
analysis using optical and IUE data. Since the previous mass determination,
the statistical and systematic accuracy of stellar spectroscopy had improved
significantly following the introduction of CCDs detectors. Thus, they
expected to obtain a smaller error on the RV determination. However, this was
hampered by the large deviations from a pure Keplerian RV curve. As did
Hiltner et al. (1972); Zuiderwijk et al. (1974); van Paradijs et al. (1977b),
they reported velocity differences with respect to the orbital fit dominated
by random excursions that they found to be autocorrelated within one night,
but not from night to night. They found substantial and very similar changes
in the shape of the profiles of all the lines. These changes in the lines plus
the velocity variability were interpreted as due to large-scale motion of the
surface induced by tidal forces due to the presence of the neutron star. This
is the same underlying physical mechanism that Tjemkes et al. (1986) proposed
to explain the irregular variations of the optical light curve. van Kerkwijk
et al. (1995) derived a $K_{\textrm{opt}}$ = $18.0-28.2\,\textrm{km
s${}^{-1}$}$ (95$\%$ confidence range) and concluded that the accuracy of the
RV study is limited by three factors: (i) the velocity excursions, (ii)
possible effects induced by the tidal deformation, and (iii) the possible
presence of systematic positive deviations in velocity close to the time of
minimum velocity. To derive more accurate constraints on the orbital
parameters, they suggested that the short- and long-term behavior of the
system should be studied to better understand the tidal interaction.

For a better understanding of the radial-velocity excursions, Barziv et al.
(2001) performed an extensive long-term spectroscopic campaign for about nine
months following one of the recommendations given by van Kerkwijk et al.
(1995). They expected that the velocity excursions would average, thus
allowing a more accurate determination of the RV amplitude than previous
works. They determined the RV from several lines using the same cross-
correlation algorithm as van Kerkwijk et al. (1995). As in previous analyses,
they found strong deviations in the RVs from those expected for a pure
Keplerian orbit, with root-mean-square amplitudes of $\sim 7\,\textrm{km
s${}^{-1}$}$ for strong lines of Si IV and N III near 4100 Å, and up to $\sim
20\,\textrm{km s${}^{-1}$}$ for weaker lines of N II and Al III near 5700 Å.
They found that systematic deviations depended on the orbital phase, with the
largest deviations observed near inferior conjunction of the neutron star and
near the phase of maximum approaching velocity. They also reported a so-called
blue excursion in the H$\delta$ data and interpreted it as a photoionization
wake. They inferred a radial-velocity amplitude $K_{\textrm{opt}}=21.7\pm
1.6\,\textrm{km s${}^{-1}$}$ with an uncertainty not much smaller than was
found in previous analyses, in which the effect of systematic phase-locked
deviations were not taken into account.

To overcome the velocity excursions issue, Quaintrell et al. (2003) carried
out a comprehensive RV study with maximum phase coverage over two consecutive
orbits of the system instead of averaging the velocity excursions over many
orbital periods. They found evidence for tidally induced nonradial
oscillations from the power spectrum of the residuals to the RV curve fit.
Moreover, when they allowed the zero phase of the RV curve to vary instead of
constraining it to the X-ray ephemeris, the fit significantly improved and
they obtained a semiamplitude of the RV curve of $K_{\textrm{opt}}=21.2\pm
0.7\,\textrm{km s${}^{-1}$}$. This is the most accurate value to date (see
Table 4 for a comparison of the determined $K_{\textrm{opt}}$). Quaintrell et
al. (2003) concluded that this apparent shift in the zero-phase may indicate
an additional RV component at the orbital period, which could be another
manifestation of the tidally induced nonradial oscillations and an additional
source of uncertainty in RV studies.

Koenigsberger et al. (2012) have explored the manner in which surface motions
induced by the tidal coupling between the blue supergiant and the neutron star
affect the RV curve (for the theory of tides in general, see Zahn, 1989). They
developed a 2D code that provides the time-dependent shape of the stellar
surface and its surface velocity field for the general case of an elliptic
orbit and asynchronous rotation without considering the effects of a
nonuniform temperature distribution over the stellar surface. They concluded
that the tidal effects on the RV curve produce orbital phase-dependent
variations in the profiles that lead to asymmetries, blue or red wings.
Moreover, the line-profile variability induces a significant variation from a
Keplerian RV curve that artificially enhances the semiamplitude. One point of
note is that a prominent feature appears in their synthetic RV curves: A blue
dip occurs shortly after the maximum. This is caused by the asymmetrical shape
of the line profiles. This feature coincides with the blue excursion reported
by Barziv et al. (2001), and it could be explained by a higher mass outflow
after periastron passage, where the tidal effects are stronger.

####  3.2.3 Quantitative spectroscopy

The stellar and wind parameters of hot massive stars such as HD 77581 are
frequently determined through quantitative spectroscopy, that is, by fitting
synthetic spectral energy distributions and normalized spectra to
observations, mainly optical and UV spectra. The best-fitting spectra and thus
parameters are frequently found by eye in a visual comparison of models and
data or by minimization on a grid of spectra.

The main parameters derived in these studies are the measured color excess
($E_{B-V}$), the (effective) stellar temperature ($T_{\star}$) from the
ionization equilibrium determined from line ratios, the surface gravity
($g_{\star}$), and the projected rotational velocity ($v_{\mathrm{rot}}\sin
i$). Stellar radii and thus luminosities are then derived using the
temperature, absolute magnitude (depending on the assumed or derived
distance), and the bolometric correction. An overview of methods and
diagnostics used for these studies and caveats to consider are given in
Martínez-Núñez et al. (2017).

Early studies in the 1980s were undertaken before codes that describe hot
stellar atmospheres were developed. These studies relied on an examination of
specific lines and a comparison to similar stars. Dupree et al. (1980) carried
out a simultaneous observation program in the X-ray, ultraviolet (IUE)222The
International Ultraviolet Explorer satellite , and optical bands. They mainly
used the comparison of P Cygni lines with the atlas of Castor & Lamers (1979)
and theoretical profiles provided by Olson (privat communication of Dupree et
al.) to estimate the terminal velocity and mass-loss rate of the line-driven
wind. Sadakane et al. (1985) analyzed IUE spectra of HD 77581 by comparing
individual lines with the corresponding lines from well-studied single B0-B1
supergiants. Other examples in which specific line features were used to
estimate some stellar wind or line properties include Prinja et al. (1990) for
terminal velocities of massive star winds or Zuiderwijk (1995) for the
rotation period and rotation velocity of HD 77581.

From the the 1990s onward, nonlocal thermodynamic equilibrium (NLTE)
hydrostatic codes began to be used to genereate synthetic spectra. In a
seminal paper, Vanbeveren et al. (1993) used a plane-parallel code by D. Kunze
for the comparison with optical spectra with high signal-to-noise ratio in the
range 4175–4525 Å to derive stellar parameters. Using the orbital parameters
of Rappaport & Joss (1983), they determined a range of stellar parameters
satisfying the orbital, X-ray eclipse, and stellar atmosphere analysis.

The NLTE codes SYNSPEC333http://tlusty.oca.eu/Synspec49/synspec.html (Hubeny
et al., 1985) and TLUSTY444http://tlusty.oca.eu/ (Hubeny & Lanz, 1995) were
used by Fraser et al. (2010) to calculate model atmosphere grids. These were
then compared with high-resolution ($R\sim 48000$) spectra obtained with
FEROS555The Fiber-fed Extended Range Optical Spectrograph
https://www.eso.org/sci/facilities/lasilla/instruments/feros.html for very
many massive stars, including HD 77581, in order to determine the atmospheric
parameters (effective temperature, surface gravity, and microturbulent
velocity), the surface nitrogen abundances, and the rotational and
macroturbulent velocities. These codes do not take the inhomogeneities of the
wind into account, however, and are optimized for hot stars with no
significant wind.

A stellar wind can considerably alter the spectral appearance and for example
spoil the derived stellar parameters such as $\log g$ if it is not taken into
account. So-called unified model atmospheres are necessary to consistently
describe the outermost layers of the star and their winds. Model atmospheres
like this are inherently in NLTE and have considerably improved over the past
decades through higher computational power and the better understanding of
physical processes in stellar atmospheres (Sander, 2017), for example, by
properly accounting for the complex effects of the millions of iron lines.
While modern codes provide a more accurate determination of the stellar and
wind parameters, a proper comparison of observed and model spectra can only be
automatized to a certain degree (see also Martínez-Núñez et al., 2017). In any
case, spectral analysis studies should be performed using as many lines as
possible because a particular single line can be affected by a variety of
parameters, and it requires detailed knowledge about which features are
affected by which parameters.

Giménez-García et al. (2016) carried out the first detailed study of HD 77581
using models with an expanding stellar atmosphere for line-driven winds
generated by the Potsdam Wolf-Rayet (PoWR)666http://www.astro.physik.uni-
potsdam.de/~wrh/PoWR/ model atmosphere code (Gräfener et al., 2002; Hamann &
Gräfener, 2003). The PoWR code takes line-driven winds, wind clumping, and a
plethora of ions into account. These were compared to the previously mentioned
IUE and FEROS spectra to derive stellar and wind parameters. The rotational
velocity, the $\beta$ exponent for the wind velocity profile (Eq. 1), and the
macroturbulence were set based on insights from previous works and on spectral
line shapes and depths. They also found indications of chemical evolution in
the star, with a moderate overabundance of He and N, together with an
underabundance of C and O.

In recent years, the PoWR code has been extended with the option to
consistently solve the hydrodynamic and statistical equations in 1D, together
with the radiative transfer, in order to obtain a hydrodynamically consistent
atmosphere stratification (Sander et al., 2017). This was applied to the Vela
X-1 system in Sander et al. (2018), who in a simplified manner also included
the effect of the X-ray radiation (see also Section 4.1.2).

Tables 5 and 7 in section 5.6 give an overview of the parameters that were
derived mainly from quantitative spectroscopy. In some cases, these results
were also used to infer information about other properties of the system (see
Sections 5.1 or 5.3 and table 4).

####  3.2.4 UV resonance lines

Hatchett & McCray (1977) predicted that an X-ray source in the stellar wind of
an early star produces changes in the photospheric ultraviolet resonance lines
of ions such as Si iv, C iv, Al iii, or N v. They also predicted that the
lines vary along the orbit and that these variations can be used to determine
the size of the stellar wind region that is affected by the X-ray ionization
induced by the compact object. However, the amplitude of the observed orbital
modulation also depends on the density and velocity of the stellar wind.

Dupree et al. (1980) presented the first UV resonance line study of the system
using IUE observations. They detected very prominent P Cygni profiles in the
resonance lines of Si iv (1393.75$\AA$ and 1402.77$\AA$) and C iv
(1548.19$\AA$ and 1550.76$\AA$). They observed variability in the profiles
over the orbital phase with more dramatic changes in the Si iv than in the C
iv profiles. They reported that when the compact object is in front of the
mass donor, around orbital phase 0.5, the P Cygni absorption is ionized to a
higher state than at eclipse, and the emission component is increased. They
concluded that the observed variability and the variation in the edge velocity
could be qualitatively explained by the Hatchett-McCray effect (hereafter, the
HM effect).

Sadakane et al. (1985) detected P Cygni profiles of the C iv and Si iv
resonance lines in the same way as Dupree et al. (1980), and in addition, they
reported the Al iii resonance doublet at 1854.72$\AA$ and 1862.79$\AA$. They
found similar profiles and variations on the C iv and Si iv as Dupree et al.
(1980). In the case of Al iii, they found a remarkable difference in the
profiles before and after mid-eclipse. They explained this behavior as due to
the complex structure of the stellar wind with two different components, the
high- and low-velocity components. The low-velocity component seemed to be
present only from mid-eclipse and during the second half of the orbital cycle.
They claimed that it was formed near the secondary and that its origin might
be a trailing wake behind the compact object. On the other hand, the high-
velocity component was observed throuhgout the orbit and was interpreted as
originating from a cooler expanding region.

Kaper et al. (1993) carried out a study to investigate the dynamical structure
of the stellar wind in the system using IUE data. They also observed orbital
variations in the blueshifted absorption part in the Si iv and C iv resonance
lines. When the compact object was in front of the optical companion, the Si
iv absorption was weaker not only at high velocities ($\sim$ -700 to -1300 km
s-1) as reported by Dupree et al. (1980), but also at low velocities ($\sim$
-300 to -500 km s-1). In the case of C iv, they found variations at high
velocities, whereas at low velocity, the variations were not significant. Like
Sadakane et al. (1985), they also found variations in the Al iii resonance,
but in this case, also at low and intermediate velocities. They interpreted
the larger Al iii variations as due to ionization effects. In addition, they
reported N v resonance lines at 1238 and 1242 $\AA$ for the first time, with
variations in the absorption part that are consistent with the changes found
in the Si iv and C iv profiles and consistent with the variability predicted
by the HM effect. Regarding the emission of the P Cygni lines, Kaper et al.
(1993) found stronger emission in the Si iv, C iv, and Al iii profiles around
the mid-eclipse phase when the X-ray source is in the line of sight. They
estimated that these orbital variations of the redshifted emission extended to
about 1000 km s-1 for the Si iv and C iv lines. The N v resonance doublet
showed different variability at the short- and in long-wavelength sides. In
the former, the strongest emission was observed during eclipse, and in the
latter, the red emission peak is remarkably stable throughout the orbit. In
conclusion, Kaper et al. (1993) reported an orbital modulation of the P Cygni
profiles not only in the high-velocity part of the profiles, but also at
intermediate and low velocities. This shows the occurrence of the HM effect
for the full range of velocities. Their observations could not be explained
with a monotonic stellar wind combined with partial ionization. They concluded
that the resonance lines can be better understood if in some regions, the wind
decelerates, in agreement with the predictions from hydrodynamical calculation
of radiation driven winds performed by Owocki et al. (1988).

van Loon et al. (2001) presented the first quantitative analysis of the UV
spectral variability in five HMXBs including HD 77581, for which they used
high-resolution spectroscopy obtained with the IUE satellite. Their synthetic
profiles were obtained using a modified version of the Sobolev exact
integration method of Lamers et al. (1987), including a nonmonotonic wind
structure, turbulence, and an extended X-ray photoionized zone. The line
profiles of Si iv, C iv, and Al iii and their orbital modulation were rather
well reproduced by their model, with an additional absorption component due to
a photoionization wake that was previously indicated by Kaper et al. (1994)
using strong optical lines such as the hydrogen Balmer lines and He lines. In
the case of the N v resonance line, their model was unable to properly
reproduce the profile, and it seemed that this line was more affected by the
photoionization wake than the other lines. From these data, they were able to
derive the extent of the region in which line-acceleration is quenched because
the energy levels that are most strongly responsible for it are depopulated by
the X-ray photoionizing feedback.

####  3.2.5 Infrared diagnostics

Observations of Vela X-1 beyond the near-infrared region of the spectrum have
been reported only rarely. Smith et al. (1990) listed HD 77581 with
significant IRAS (Infrared Astronomical Satellite) fluxes in their collection
of 81 X-ray binary counterparts, but did not discuss this further in their
study, which focused on infrared emission from accretion disks. Friedemann et
al. (1996) noted that the IRAS fluxes reveal cold dust and that the binary may
be surrounded by a fossil dust shell.

As detailed in Sect. 5.1.2, Vela X-1 is associated with a wind bow shock in
the interstellar medium. This structure was first discovered in a H$\alpha$
image (Kaper et al., 1997), but has been studied in the infrared based on IRAS
(Huthoff & Kaper, 2002), Spitzer (Gvaramadze et al., 2011), and WISE (Maíz
Apellániz et al., 2018) observations.

Choquet et al. (2014) obtained infrared interferometric observations of Vela
X-1 with the VLTI in the K band (2.2$\mu$m) in 2010 and in the H band
(1.6$\mu$m) in 2012. The second observations covered different orbital phases
(0.89, 0.11, 0.33, and 0.55), but with no significant variation in the
interferometric observables. In the K band, a structure of $8\pm 3R_{\star}$
was resolved, while two years later, a centro-symmetric structure of radius
$2.0^{+0.7}_{-1.2}R_{\star}$ was found in the H band. The authors offered
three different possible explanations for the difference in derived diameters:
(1) A strong temperature gradient in the supergiant wind, where hot material
at 1720 K is more compact than material at 1350 K; (2) a diffuse gaseous shell
observed in 2010, which had diffused two years later; or (3) the structure
observed in the H band was the stellar photosphere and not the supergiant
wind.

The chemical evolution noted for HD 77581 could indicate that shock chemistry
may be dominant in the surrounding interstellar medium. For other regions of
the sky, it has been shown that relevant dynamical and physical parameters of
the shocked gas can be determined by analysing its far-infrared (far-IR)
emission lines (Lerate et al., 2008, e.g.). In the same vein, the study of
far-infrared emission lines could provide a better insight into the wind
clumps and structure.

The region of Vela X-1 has been briefly observed by Herschel in 2012, but no
associated publications are found in the Herschel Science Archive.

###  3.3 Radio observations

On 15 May, 2018, van den Eijnden et al. (in prep.) observed Vela X-1 with the
Australia Telescope Compact Array. This four-hour observation yielded the
first radio detection of the source at flux densities of $92\pm 11$ and
$122\pm 10$ $\mu$Jy at $5.5$ and $9$ GHz, respectively. These measured flux
densities imply a spectral index measurement of $\alpha=0.56\pm 0.36$, where
$S_{\nu}\propto\nu^{\alpha}$. The observation was performed completely during
the eclipse of the neutron star, at orbital phases of 0.92–0.94.

With the currently known radio properties of Vela X-1, it is challenging to
determine the dominant radio emission process. While in LMXBs radio emission
has been established as a signature of synchrotron-emitting jets (except for
the handful of LMXBs hosting a strongly magnetized neutron star), the massive
donor and its wind in HMXBs can also contribute to the observed radio
emission. The properties of Vela X-1 fit both scenarios: The spectral index is
consistent with that of a steady, compact jet, as seen typically in hard
states of persistent low-mass X-ray binaries ($\alpha=0$–$0.7$; Fender et al.,
2004). The radio luminosity of Vela X-1, interpolated to be $\sim 2\times
10^{27}$ erg s-1 at 6 GHz, is also similar to the radio emission of other jets
launched by strongly magnetized (e.g. $B>10^{12}$ G) accreting neutron stars
(van den Eijnden et al., 2018a, b, c, 2019). Similarly, the spectral index
fits the predicted spectral index for thermal radio emission from an ionized
stellar wind ($\alpha=0.6$; Wright & Barlow, 1975; Güdel, 2002). The flux
densities also fit within their uncertainties the predicted brightness of such
a wind following the formalism in Wright & Barlow (1975), assuming a mass-loss
rate of $10^{-6}$ $M_{\odot}$ yr-1 and a terminal wind velocity of $\sim 700$
km s-1 (Grinberg et al., 2017).

Both scenarios could have significant implications: First, strongly magnetized
neutron stars were long thought to be incapable of launching jets (e.g.,
Fender & Hendry, 2000; Migliari & Fender, 2006; Massi & Kaufman Bernadó,
2008), until the recent detection of such a jet showed otherwise (van den
Eijnden et al., 2018c). If Vela X-1 also launches a jet, it would add to the
still small but growing sample of strongly magnetized neutron stars launching
jets (van den Eijnden et al., in prep.). Alternatively, if a wind is observed
in radio, it would only be the second HMXB with a radio wind detection (after
GX 301-2; Pestalozzi et al. 2009). Based on its high inclination and clear
orbital variation of the wind absorption, Vela X-1 would provide a unique and
powerful testbed in which to explore the properties of stellar winds in HMXBs
in radio, while at the same time, radio observations would provide a new and
complementary view of Vela X-1 itself.

At the time of writing, however, the first challenge is to probe the origin of
the detected radio emission: jet, wind, or a combination of both? Detailed
radio and X-ray monitoring over several orbital phases will be vital for this
purpose, revealing the variability of the radio properties and allowing a
comparison between accretion, X-ray, wind, and radio properties along the
orbit.

##  4 Models of Vela X-1

Like most wind-fed HMXBs, Vela X-1 is a system that can only be appreciated
with multiphysics and multiscale models (see Figure 2). Many efforts have been
made to capture the predominant mechanisms at each scale. In simplified
geometries assuming an isotropic stellar wind (1D models) or focusing on the
orbital plane (2D models), some semianalytic results were derived. However,
since the late 1980s, numerical simulations have been the privileged tool for
capturing the whole 3D complexity of this archetypical system and provide
suitable models for interpreting the observations.

###  4.1 Stellar wind

####  4.1.1 Hydrodynamic structure

The donor star is not isolated in Vela X-1. The presence of the orbiting
neutron star breaks the spherical symmetry of the problem, and the wind cannot
be fully described in a 1D framework. If the wind speed were much higher than
the orbital speed, the wind would only be affected by the gravitational field
of the neutron star in a region far smaller than the orbital separation (see
the discussion of the accretion radius in Section 4.2). Except in the thin
wake of the accretor, which is axisymmetric with respect to the axis joining
the two bodies, the 1D representation would hold. On this implicit assumption
Watanabe et al. (2006) relied when they empirically introduced a cone-shaped
dense and cold cloud in the wake of the neutron star to explain the excess of
absorption at inferior conjunction. In the same spirit, Fryxell et al. (1987)
performed wind tunnel simulations: They ran 2D axisymmetric hydrodynamic
simulations on a spherical grid of a planar adiabatic flow coming from
infinity with a supersonic relative speed that was deflected by the
gravitational field of an accretor. This configuration corresponds to the
problem that was first addressed by Hoyle & Lyttleton (1939) and Bondi & Hoyle
(1944) in ballistic terms, the BHL configuration. Fryxell et al. (1987) found
a good agreement with the BHL predictions and that an axisymmetric bow shock
forms around the accretor, with an opening angle depending on the Mach number
of the flow. Provided the orbital effects are negligible, that is, when the
wind speed is much higher than the orbital speed, the flow in a wind-fed HMXB
can be divided into two parts: A purely radial wind flowing away from the
donor star, and in the immediate vicinity and in the wake of the compact
object, a flow that is well described by the BHL solution or by its
hydrodynamics counterpart.

Nevertheless, in Vela X-1, a more realistic wind speed at the orbital
separation of about or even lower than the orbital speed, such as the one used
by Watanabe et al. (2006) and Sander et al. (2018) (see Section 4.1.2),
indicates an anisotropic wind that surrounds the donor star. The outflow is
strongly shaped by the orbital motion, as is corroborated by the systematic
asymmetries observed in the absorbing column density profile777An eccentric
orbit in an isotropic wind would also result in an asymmetric $N_{\textrm{H}}$
profile, although of more limited contrast than observed. between orbital
phases 0-0.5 and 0.5-1 (see Secion3.1.3). In this case, a fully 3D model is
required to describe the hydrodynamic structure of the wind and interpret the
evolution of the observations with the orbital phase in a high-inclination
HMXB such as Vela X-1. Bessell et al. (1975) first highlighted the anisotropic
distribution of material in the orbital plane induced by the tidal forces. It
was later confirmed by 2D hydrodynamic simulations in the equatorial plane of
a spherical mesh by Blondin et al. (1990, 1991). They characterized the
gravitational and radiative impacts of the neutron star on the wind: They
showed that the nonsteady and overdense envelope of the bow shock that is
formed upstream trails the neutron star throughout its orbit. This structure
is referred to as the accretion wake. Blondin et al. (1991) further showed
that in simulations of Vela X-1, the wind was highly focused into a steady
tidal stream between the star and the neutron star. This feature is clearly
visible in the 2D simulations by Manousakis et al. (2012). On the other hand,
wind tunnel simulations centered on the neutron star but with a nonplanar
inflow were previously performed on 2D cylindrical $(r,\phi)$ grids (Taam &
Fryxell, 1988; Fryxell & Taam, 1988; Taam & Fryxell, 1989). The aim was to
parameterize the shearing in velocity and/or mass density of the inflow that
is caused by the orbital effects. They identified an oscillation of the
accretion wake around the main axis of the cylinder, which they called the
flip-flop instability (see also Matsuda et al., 1987, 1991, for the seminal
simulations in which this effect was determined). First thought to be
responsible for the formation of transient disks in wind-fed HMXBs, this
feature was later found to be a numerical artifact due to the unrealistic
geometry of the mesh and to the unphysically large size of the sink particle
standing for the accretor (Ishii et al., 1993; Ruffert, 1999; Blondin &
Raymer, 2012). In addition to the tidal stream and the accretion wake, a
spiral-shaped density enhancement is to be expected at the side of the star
opposite to the neutron star due to the tidal forces (El Mellah et al.,
2019b). With the 3D framework they worked with, these authors also solved the
dynamics of the wind off the orbital plane and found a significant flattening
of the wind in the orbital plane when the ratio of the wind speed to the
orbital speed was similar to what is expected for Vela X-1.

In Vela X-1, the microstructure of the wind is still poorly constrained,
although Giménez-García et al. (2016) measured a density contrast in the wind
that is consistent with theoretical predictions. Observational diagnostics
that depend quadratically on density (e.g., H$\alpha$ line and thermal radio
emission) lead to overestimated mass-loss rates when clumping is not accounted
for (Sundqvist et al., 2011). Conversely, the mass-loss rates deduced from
diagnostics that depend linearly on density (e.g., UV resonance lines) can be
underestimated (Sundqvist & Puls, 2018). Therefore the current uncertainties
on the clumping factor prevent us from drawing a definitive conclusion
concerning the stellar mass-loss rate. In Vela X-1, the size of the clumps
derived from simulations and from the two main observational diagnostics at
our disposal, absorption variability and X-ray flares, differs. While the
first is thought to be due to unaccreted clumps passing the line of sight
(see, e.g., the observations and model by Grinberg et al., 2017; El Mellah et
al., 2020b, respectively), the second have been proposed to be due to the
serendipitous capture of a clump by the neutron star (Fürst et al., 2010;
Martínez-Núñez et al., 2014). However, hydrodynamics simulations by El Mellah
et al. (2018) showed that the bow shock around the accretor significantly
lowers the variability compared to a purely ballistic model of the clump
capture. It suggested that although intrinsic X-ray flares might be triggered
by clump capture, the amount of mass involved in a flare was probably much
more important than the amount of mass contained in a single clump. A better
understanding of the whole accretion process, from the orbital scale all the
way down the accretion columns where most of the X-rays we observe originate,
is necessary to determine the connection between the flares and the clumpiness
of the wind.

####  4.1.2 Ionization structure and inhibition of wind acceleration

It is commonly accepted that the X-ray irradiation from the neutron star
significantly affects the wind velocity profile, but it has proven very
challenging to quantify the extent of this effect on the dynamics of the flow.
In a seminal paper, Fransson & Fabian (1980) described how the photoionization
of the wind by the accretor in high-mass X-ray binaries could lead to the
formation of a shock between the accelerating wind and the stalling
photoionized plasma. Blondin et al. (1990) quantified this effect with 2D
numerical simulations. They modeled the effect of the X-ray ionizing feedback
on the hydrodynamic structure of the wind. Assuming that wind acceleration was
unaltered below $\xi_{\text{crit}}=10^{2.5}$erg cm s-1 and fully inhibited
above, they showed that the wind speed in the direction of the neutron star
was lowered by a factor of $\sim$2\. Because of the stalling of the line-
driven wind that is caused by X-ray photoionization, in some simulations, the
wind formed a trailing spiral density enhancement between the neutron star and
the star. Kaper et al. (1994) relied on this photoionization wake to interpret
an absorption component that they were unable to attribute to the accretion
wake or to the tidal stream.

Although wind acceleration is inhibited in the limit case of a high X-ray
illumination, for intermediate fluxes, it is not a step function nor even a
monotonous dependence of the X-ray luminosity. For instance, as already noted
byMacGregor & Vitello (1982), for intermediate X-ray photoionizing fluxes, the
wind would overall be more efficient at absorbing UV photons, and the line-
driven acceleration would even become higher. In Vela X-1, a few studies
initiated by Sako et al. (1999), went beyond the on/off switch assumption in
different ways. In a preliminary approach, Watanabe et al. (2006) assumed
simplified descriptions of the line acceleration, the radiative transfer, and
the ionization structure to compute 1D velocity profiles. They found a wind
speed at the orbital separation in Vela X-1 of 180km$\cdot$s-1 for
$L_{\text{X}}=$3.5$\times$1036erg$\cdot$s-1 instead of 570km$\cdot$s-1 without
X-ray ionizing feedback. However, they did not iterate the radiative transfer
and computation of the ionization structure and only accounted for a very
limited number of line transitions.

More advanced treatments have emerged in the past decade to determine the
effect X-rays have on the wind acceleration. In these methods, realistic
stellar atmosphere models accounting for the inherent NLTE conditions and the
moving outer layers were irradiated by an external X-ray source in order to
study the imprint of the different elements and ions on the line acceleration.
Krtička et al. (2012) derived wind solutions with the METUJE code (Krtička &
Kubát, 2017) for different axis in the orbital plane and iterated the
radiative transfer and statistical equilibrium equations to obtain a
consistent wind ionization structure. In a follow-up study, Krtička et al.
(2018) accounted for optically thin wind clumping and obtained similar and
very low values of the critical ionization parameter, above which line-
acceleration was significantly inhibited, ranging between 5 and 25
erg$\cdot$cm$\cdot$s-1. In both cases, they computed wind velocities that were
severely lower in the direction of the neutron star in Vela X-1, with a speed
that hardly peaked at 100km$\cdot$s-1. Following a spectral analysis by
Giménez-García et al. (2016) that yielded updated parameters for the donor
star and wind in Vela X-1, Sander et al. (2018) performed an analysis with the
dynamically consistent branch of the PoWR code (Sander et al., 2017) to study
the effect of X-ray irradiation on individual driving contributions and
spectral lines. They showed that with the updated parameters, an X-ray
ionizing source even as high as 1037erg$/$s does not necessarily inhibit wind
acceleration, but might even enhance it upstream the neutron star. While both
METUJE and PoWR perform calculations in 1D, their results demonstrate that the
wind speed could be so low that the outflow significantly departs from a
spherically symmetric wind. However, the computational cost of these methods
precludes any 3D direct computation in the next years.

Finally, it is worth mentioning that a preliminary global model of the
radiation-driven photoionized wind was sketched and applied to Vela X-1 in
Mauche et al. (2007) and Mauche et al. (2008). The authors combined different
numerical tools. First, 2D and 3D hydrodynamical simulations were performed
with FLASH (Fryxell et al., 2000). Second, the wind ionization structure was
computed with the photoionization code XSTAR in order to derive the local
heating and cooling rate, along with the ionization fractions of the different
ions (Kallman & Bautista, 2001). Third, the rich HULLAC database provided the
emission models for the X-ray photoionized plasma (Bar-Shalom et al., 2001).
Finally, a Monte Carlo radiative transport computation enabled them to derive
detailed synthetic X-ray spectra.

###  4.2 Accreted flow

In Vela X-1, the ionization state and thus the whole hydrodynamic structure of
the flow depends on the intensity of the X-ray emission from the neutron star,
fed by the wind it illuminates. It emphasizes the need addressed by various
authors for simultaneously modeling the orbital scale and the accretion
process itself in order to achieve a fully consistent model that can connect
the X-ray luminosity to the mass accretion rate. The accretion of stellar
material by the neutron star proceeds in different steps. A mass transfer
mechanism must first ensure that stellar material is brought into the Roche
lobe of the neutron star if the flow has a negligible amount of specific
kinetic energy, or is within the accretion radius given by equation (3)
otherwise. These two characteristic length scales far exceed the other
relevant scales of accretion, therefore we can first treat the mass transfer
alone (Section 4.2.1). Then, accretion can proceed onto the magnetosphere of
the neutron star and eventually carry the plasma all the way down the X-ray
emitting regions near the magnetic poles (Section 4.2.2).

####  4.2.1 Mass transfer mechanisms

The donor star in Vela X-1 is unlikely to undergo RLOF. In X-ray binaries,
mass transfer through RLOF usually leads to higher X-ray luminosities than in
Vela X-1. Moreover, RLOF leads to the formation of a permanent accretion disk,
which has not been observed in Vela X-1. The torques by such a disk would also
spin the neutron star up, as is the case in Be X-ray binaries, but this
contradicts the measured long neutron star spin period. An additional strong
argument against mass transfer through RLOF comes from the stable orbital
period (Falanga et al., 2015) because no outflow from the vicinity of the
neutron star is observed and conservative mass transfer from a high-mass star
to a low-mass accretor would lead to a quickly shrinking orbit (Quast et al.,
2019; El Mellah et al., 2020a). Instead, accretion of the stellar wind is to
be preferred as the dominant mass transfer mechanism.

Wind accretion as described by the BHL formalism might be a misleading
alternative, however. Equation (4) has been widely used to interpret the
intrinsic X-ray luminosity of Vela X-1 and its variations, but fundamental
caveats must be acknowledged, especially because the wind speed at the orbital
separation is not high compared to the orbital speed. It is plausible that a
hybrid mass transfer mechanism is taking place in Vela X-1 (called wind-RLOF
and introduced by Mohamed & Podsiadlowski, 2007), where the wind is
significantly disrupted and the geometry of the problem no longer obeys that
of the planar BHL (see Section 2.1). For instance, the wind is compressed in
the orbital plane by the motion of the two bodies, so that the density in
equation 4 cannot be deduced from the continuity equation in spherical
geometry. An approach to obtain more accurate mass transfer rates was
introduced by El Mellah & Casse (2017), who assumed a donor star in
synchronous rotation with the orbit and applied their model to Vela X-1:
Because the wind is highly supersonic, ballistic streamlines reveal how beamed
the flow is toward the neutron star, and density-enhanced (depleted) regions
are identified where the streamlines diverge slower (faster) than in the
spherically geometric case. The streamlines that intercept within the Roche
lobe of the neutron star indicate the fraction of the stellar wind that will
take part in the accretion process. They also enable us to compute the net
specific angular momentum of the flow, which significantly departs from zero:
An even more important reason to avoid the use of equation (4) in Vela X-1 is
indeed that the flow might carry enough angular momentum to circularize before
it reaches the neutron star magnetosphere, in contrast to the intrinsic
assumption of cancelling angular momenta in the BHL picture.

The fact that the donor star does not fill its Roche lobe sets a stringent
constraint on the orbital inclination angle $i$. If the eclipse duration
relative to the orbital period is $\Delta t/P$, the ratio of the stellar
radius $R$ to the orbital separation $a$ reads

$\frac{\Delta
t}{P}=\frac{1}{\pi}\arcsin\left(\frac{\sqrt{R^{2}-a^{2}\cos^{2}i}}{a\sin
i}\right),$ (8)

where for now, we neglect the eccentricity of the orbit and the tidal
deformations of the star. The condition of absence of RLOF yields a lower
limit on the inclination angle,

$i>i_{\text{min}}=\arcsin\left[\frac{\sqrt{1-\mathcal{E}^{2}(q)}}{\cos\left(\pi{\Delta
t}/{P}\right)}\right],$ (9)

where $\mathcal{E}$ is the estimate of the ratio of the stellar Roche-lobe
radius to the orbital separation provided by Eggleton (1983) and function of
the mass ratio alone. For a given range of $q$ and $\Delta t/P$, we can thus
deduce a range of minimum orbital inclination angles below which the donor
star would overflow its Roche lobe. In the upper panel in Figure 11, we plot
the minimum orbital inclination angle to avoid RLOF as a function of the mass
ratio for realistic eclipse durations. At low mass ratios and for long eclipse
durations ($q<12$ and $\Delta t/P>0.2$), the donor star would fill its Roche
lobe even if the system were seen edge-on. Currently, the upper limit set on
the eclipse duration is closer to $\Delta t/P=0.19$ (Kreykenbohm et al.,
2008), which would indicate a lower limit for the orbital inclination of
$\sim$74∘. This conclusion must be tempered, however, because the orbit is not
circular. Moreover, the proximity of the donor star to filling its Roche lobe
means that the star should not be assumed to be spherical, and only an
effective stellar radius can be defined.

Figure 11: Minimum inclination angle of the orbit for avoiding RLOF as a
function of the mass ratio and for different durations of the eclipse relative
to the orbital period. In the upper panel, a circular orbit is assumed, and in
the lower panel, we account for an eccentricity of 0.0898 and enforce no RLOF
even at periastron.

Furthermore, the eccentricity of the orbit in Vela X-1 leads to a filling
factor that varies with the orbital phase: Because the orbit is eccentric, the
orbital separation varies, and so does the filling factor. When we enforce the
more severe constraint to avoid RLOF even at periastron, we obtain the minimum
inclination angles displayed in the lower panel in Figure 11, which are
significantly higher. Two conclusions can be drawn from this result. Either
the system presents an orbital inclination and a mass ratio that lie at the
upper edge of the intervals derived from observational diagnostics not based
on the assumption of the absence of RLOF, or there are regular episodes of
intermittent RLOF at periastron. In Figure 12 we show the evolution of the
Roche-lobe filling factor of the donor star as a function of the orbital phase
for three different sets of system parameters derived from observations by
Quaintrell et al. (2003), Rawls et al. (2011), and Falanga et al. (2015). At a
given orbital phase, the uncertainty ranges are upper limits, obtained by
considering the extreme values obtained without accounting for the
correlations between parameters. The striking result is that observationally
derived parameters tend to give filling factors above 1. When we take the
absence of RLOF for granted, it means that we can set stronger constrains on
the system parameters. With an eccentricity of 0.0898, however, the filling
factor necessarily varies by $\sim$15% throughout the orbit. This major
variation implies a mass transfer mechanism that might be modulated throughout
the orbit, with a short phase of mass transfer through RLOF at periastron.

Figure 12: Stellar filling factor, i.e., the ratio of the stellar radius to
the Roche-lobe radius, as a function of the orbital phase (0 is mid-eclipse).
Each colored stripe is obtained for a different set of extreme values
observationally derived by Quaintrell et al. (2003), Rawls et al. (2011), and
Falanga et al. (2015).

####  4.2.2 Magnetospheric accretion and induced torques

Following its formation, a pulsar spins down first due to its autonomous
magnetodipole radiation (Pacini, 1968; Gunn & Ostriker, 1969). After a few
thousand years, it spins down at a higher pace due to the torques associated
with the repelling of infalling stellar material that is ejected during the
propeller phase (Francischelli et al., 2002). After this, the magnetospheric
and corotation radii are expected to be similar because the torques applied by
the accreted flow onto the neutron star magnetosphere tend to spin it up (spin
it down) if the corotation radius is larger (smaller) than the magnetospheric
radius (see Ho et al., 2020, for a study of the long-term evolution toward
spin equilibrium that covers these three phases). Because this equilibrium
phase is reached within only a few 10 kyr, a first attempt to estimate the
neutron star magnetic field from its spin period and the mass accretion rate
can be made assuming that the magnetospheric radius and the corotation radius
match. The corotation radius of the neutron star in Vela X-1 can be accurately
deduced from the mass of the neutron star and from its spin period,

$R_{co}=\left(\frac{GMP^{2}}{4\pi^{2}}\right)^{1/3}\sim 8.1\times
10^{9}\text{cm}\left(\frac{M}{2\,M_{\odot}}\right)^{1/3}\left(\frac{P}{283\,\text{s}}\right)^{2/3}.$
(10)

In addition to being stable over timescales shorter than years, the corotation
radius is known with a 10% precision given its low dependence on the mass and
spin period of the neutron star. On the other hand, the expression that
determines the magnetospheric radius depends on the geometry of the accretion
flow (planar, disk, spherical), but an order-of-magnitude estimate can be
obtained from (Lamb et al., 1973; Wang, 1996)

$\begin{split}R_{mag}\sim 4.9\times
10^{8}\text{cm}&\left(\frac{\dot{M}_{acc}}{2\times
10^{-10}\text{M}_{\odot}\cdot\text{yr}^{-1}}\right)^{-2/7}\left(\frac{M}{2\,M_{\odot}}\right)^{-1/7}...\\\
&...\left(\frac{B}{10^{12}\,\text{G}}\right)^{4/7}\left(\frac{R_{NS}}{12\,\text{km}}\right)^{12/7}\end{split},$
(11)

where the inertia moment of the neutron star is computed assuming a uniform
sphere (i.e., $2MR^{2}/5$), where an accretion rate of $2\times
10^{-10}\text{M}_{\odot}\text{yr}^{-1}$ corresponds to an accretion luminosity
of 4$\times 10^{36}$erg s-1 for a conversion factor of 35%, and where we used
the definition of a magnetic dipole $\mu$ as $BR^{3}/2$.

At equilibrium, the magnetic field should be such that these two radii are
equal. With numerical values coherent with those deduced for Vela X-1, we
obtain an estimate of the magnetic field of $B\sim 2.2\times 10^{14}\text{G}$.
This value is at odds with what was measured relying on the CRSFs diagnostics
(see section 3.1.7). Fürst et al. (2014) measured a cyclotron line energy
corresponding to a magnetic field $B\sim 2.6\times 10^{12}\text{G}$, which is
two orders of magnitude smaller than the value derived assuming spin
equilibrium. The serious mismatch between the magnetic field values
independently deduced from the corotation radius and from CRSFs in Vela X-1
would mean that the neutron star has been spun down to a spin period much
higher than the equilibrium spin period set by the torques predicted by Ghosh
& Lamb (1979b) and corresponding to the estimated mass accretion rate. This is
probably the strongest theoretical point that can be made against the presence
of a permanent disk beyond the neutron star magnetosphere in Vela X-1. The
magnetic field required to account for the observed neutron star spin
variations is also higher than the field measured by the CRSFs diagnostics by
two orders of magnitude (see Staubert et al., 2019, for a detailed
discussion).

Alternatively, this overbraking could be indicative of a phase in the past in
which stellar material was supplied to the neutron star at a much lower rate
than today, and with a high amount of specific angular momentum, which is
suggestive of a slow wind (Wang & Tong, 2020; Ho et al., 2020). The system
would have been in the propeller regime, and the neutron star would have spun
down until the mass accretion rate suddenly increased by so much that the
magnetospheric radius became smaller than the corotation radius, as is
currently observed if the magnetic field measured through CRSFs is accurate.
We must acknowledge, however, that this scenario requires a fine-tuning of the
increase in mass accretion rate, of the initial magnetic field, and of the
magnetic field decay, which is unlikely. In what follows, we rely on the
magnetic field value obtained by direct measures through CRSFs, that is,
$B\sim 2.6\times 10^{12}\text{G}$.

The first conclusion to draw from this value of the magnetic field is that in
Vela X-1, the neutron star is never expected to be in the propeller regime. In
equation (11), the X-ray luminosity that is required for the magnetospheric
radius to be as large as the corotation radius is about 1032erg s-1. This
threshold is too low to be compatible even with the lowest observed X-ray flux
values (Liao et al., 2020). For the propeller effect to be triggered, a drop
in mass accretion rate by 3 to 4 orders of magnitude compared to its median
value would be needed (Bozzo et al., 2008; Doroshenko et al., 2011). For a
realistic wind microstructure, the inter-clump environment is therefore
insufficiently empty to cause such a drop (Sundqvist & Puls, 2018),
notwithstanding the lowered variability due to the mixing of the wind material
at the hydrodynamics shock (El Mellah et al., 2018).

In spite of the absence of the propeller accretion regime in Vela X-1, we do
observe month- to year-long episodes of spin-up and spin-down. Each phase
occurs at a steady spin period rate: $\sim$7$\times$10-4s day-1 to
$\sim\times$10-3s day-1 for spin-up and $\sim$10-4s day-1 for spin-down
(Malacaria et al., 2020). They are separated by sudden torque inversions, and
the net spinning up over the past decades is compatible with zero. In neutron-
star-hosting HMXBs in general and in Vela X-1 in particular, the spinning up