11institutetext: Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001
Leuven, Belgium
11email<EMAIL_ADDRESS>22institutetext: National Solar
Observatory, 22 Ohi’a Ku Street, Makawao, HI 96768, USA 33institutetext: Anton
Pannekoek Institute for Astronomy, University of Amsterdam, Postbus 94249,
1090 GE Amsterdam, The Netherlands 44institutetext: European Southern
Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching bei München,
Germany 55institutetext: European Southern Observatory, Alonso de Cordova
3107, Vitacura, Casilla 19001, Santiago de Chile, Chile 66institutetext:
Institut de Planétologie et d’Astrophysique de Grenoble, 414 Rue de la
Piscine, 38400 Saint-Martin-d’Hères, France

# Binary-object spectral-synthesis in 3D (BOSS-3D)

Modelling Hα emission in the enigmatic multiple system LB-1

L. Hennicker Binary-object spectral-synthesis in 3D (BOSS-3D)Binary-object
spectral-synthesis in 3D (BOSS-3D) N. D. Kee Binary-object spectral-synthesis
in 3D (BOSS-3D)Binary-object spectral-synthesis in 3D (BOSS-3D) T. Shenar
Binary-object spectral-synthesis in 3D (BOSS-3D)Binary-object spectral-
synthesis in 3D (BOSS-3D)Binary-object spectral-synthesis in 3D
(BOSS-3D)Binary-object spectral-synthesis in 3D (BOSS-3D) J. Bodensteiner
Binary-object spectral-synthesis in 3D (BOSS-3D)Binary-object spectral-
synthesis in 3D (BOSS-3D)Binary-object spectral-synthesis in 3D
(BOSS-3D)Binary-object spectral-synthesis in 3D (BOSS-3D) M. Abdul-Masih
Binary-object spectral-synthesis in 3D (BOSS-3D)Binary-object spectral-
synthesis in 3D (BOSS-3D) I. El Mellah Binary-object spectral-synthesis in 3D
(BOSS-3D)Binary-object spectral-synthesis in 3D (BOSS-3D) H. Sana Binary-
object spectral-synthesis in 3D (BOSS-3D)Binary-object spectral-synthesis in
3D (BOSS-3D) J. O. Sundqvist Binary-object spectral-synthesis in 3D
(BOSS-3D)Binary-object spectral-synthesis in 3D (BOSS-3D)

(Received 20 July 2021; Accepted 23 November 2021)

###### Abstract

Context. To quantitatively decode the information stored within an observed
spectrum, detailed modelling of the physical state and accurate radiative
transfer solution schemes are required. The accuracy of the model is then
typically evaluated by comparing the calculated and observed spectra. In the
analysis of stellar spectra, the numerical model often needs to account for
binary companions and 3D structures in the stellar envelopes. The enigmatic
binary (or multiple) system LB-1 constitutes a perfect example of such a
complex multi-D problem. Thus far, the LB-1 system has been indirectly
investigated by 1D stellar-atmosphere codes and by spectral disentangling
techniques, yielding differing conclusions about the nature of the system
(e.g. a B-star and black-hole binary with an accretion disc around the black
hole or a stripped-star and Be-star binary system have been proposed).

Aims. To improve our understanding of the LB-1 system, we directly modelled
the phase-dependent Hα line profiles of this system. To this end, we developed
and present a multi-purpose binary-object spectral-synthesis code in 3D
(BOSS-3D).

Methods. BOSS-3D calculates synthetic line profiles for a given state of the
circumstellar material. The standard $pz$-geometry commonly used for single
stars is extended by defining individual coordinate systems for each involved
object and by accounting for the appropriate coordinate transformations. The
code is then applied to the LB-1 system, considering two main hypotheses, a
binary containing a stripped star and Be star, or a B star and a black hole
with a disc.

Results. Comparing these two scenarios, neither model can reproduce the
detailed phase-dependent shape of the Hα line profiles. A satisfactory match
with the observations, however, is obtained by invoking a disc around the
primary object in addition to the Be-star disc or the black-hole accretion
disc.

Conclusions. The developed code can be used to model synthetic line profiles
for a wide variety of binary systems, ranging from transit spectra of
planetary atmospheres, to post-asymptotic giant branch binaries including
circumstellar and circumbinary discs and massive-star binaries with stellar
winds and disc systems. For the LB-1 system, our modelling provides strong
evidence that each object in the system contains a disc-like structure.

###### Key Words.:

radiative transfer – methods: numerical – stars: emission-line, Be – stars:
black holes – binaries: spectroscopic

##  1 Introduction

The interpretation of observed radiation emerging, for instance, from
galaxies, individual stars, or planetary systems, is key to understanding
fundamental properties of the Universe. Deep insight into the underlying
physics of an observed object can be obtained by analysing the electromagnetic
spectrum of the emitting object. Decoding the information stored within a
spectrum, however, is not a simple task. Deducing the fundamental parameters
of a single star (i.e. $L_{\ast}$, $T_{\rm eff}$, $R_{\ast}$), for example,
requires detailed modelling of the stellar atmosphere, possibly relaxing the
assumption of local thermal equilibrium (LTE) and accounting for wind
outflows. For a given model atmosphere, synthetic spectra then need to be
calculated and compared with the observations. In the OB-star regime (on which
this paper focusses), state-of-the-art atmospheric modelling and spectral
synthesis codes typically assume spherical symmetry (e.g. PHOENIX: Hauschildt
1992; CMFGEN: Hillier & Miller 1998, Hillier 2012; WM-basic: Pauldrach et al.
2001; PoWR: Hamann & Gräfener 2003, Sander et al. 2017; FASTWIND: Sundqvist &
Puls 2018, Puls et al. 2020).

Many stars, however, can deviate from spherical symmetry, with perturbations
typically induced by magnetic fields, surface (and wind) distortions from
rotation, surrounding discs, and binarity or multiplicity effects. Due to the
complexity of including full non-LTE occupation numbers within multi-D
calculations, corresponding spectral synthesis codes are only gradually
developed by applying different solution methods from finite-volume methods
(FVM) to short- and long-characteristics methods (SC and LC) to Monte-Carlo
methods (MC). A non-exhaustive list of 3D radiative transfer codes accounting
for supersonic velocity fields includes Wind3D (Lobel & Blomme 2008, FVM), a
code developed by Hennicker et al. 2020 (SC), Phoenix/3D (Hauschildt & Baron
2006 and other papers in this series, LC), and HDUST (Carciofi & Bjorkman
2006, MC), where the first two examples currently apply a two-level-atom (TLA)
approach, whereas the last two include a multi-level description of an atom
(of typically one particular species). When the atomic level populations are
obtained by such codes, the synthetic spectra can be calculated by solving the
equation of radiative transfer along the direction to the observer, for
instance, by applying the LC method in a cylindrical coordinate system based
on Lamers et al. (1987), Busche & Hillier (2005), Sundqvist et al. (2012) (see
also Hennicker et al. 2018).

Thus far, the above described codes and methods have been tailored to the
spectral synthesis of single stars in 1D or 3D. To our knowledge, there exists
no viable multi-D alternative to calculating synthetic spectra from massive-
star binary systems with circumstellar material around them. For example, the
SPAMMS-code by Abdul-Masih et al. (2020b) relies on patching the emergent
intensities as obtained from several 1D spherically symmetric FASTWIND models.
Thus, the effect of rays propagating through the winds of both individual
stars is neglected. Moreover, such models cannot be used for arbitrary 3D
structures such as circumstellar discs.

In this paper, we thus develop a general purpose spectral synthesis tool for
binary systems (BOSS-3D, i.e. binary-object spectral-synthesis in 3D) by
extending the solution method for single stars from Hennicker et al. (2018)
using simple coordinate transformations. For given atmospheric and
circumstellar structure(s) with (thus far approximate) occupation numbers,
this code accounts for large-scale, possibly structured outflows or discs
within both involved systems, thus circumventing the aforementioned
shortcomings of previous methods. By allowing for different length scales of
the involved objects, the code can also be used to analyse transit spectra of
planetary systems. The basic idea of our method relies on defining two
(independent) coordinate systems for the individual objects, and applying the
standard single-star $pz$-geometry to each of the systems. To correctly
account for overlapping coordinates (e.g. during transits), the individual
systems are merged by triangulation.

As a first application of the newly developed algorithm, we considered the
binary (or multiple) system LB-1. The spectrum of this system clearly shows an
anti-phase behaviour of stellar absorption lines against the line wings of a
broad Hα emission (Liu et al. 2019). Yet the interpretation of these findings
is still under debate. Based on these observations and the corresponding
orbital reconstruction, Liu et al. (2019) proposed a binary system with period
$P\approx 79\,{\rm d}$, consisting of a B star (visible component) and a
$70\,M_{\odot}$ black hole (BH) with associated disc yielding the Hα emission.
Since it is difficult to explain the origin of a $70\,M_{\odot}$ BH in the
Milky Way from stellar evolution theory (though, see Belczynski et al. 2020),
there has been vivid discussion about LB-1 in the literature. While Simón-Díaz
et al. (2020) re-analysed the visible component (B star, primary object) and
found a significantly lower mass for the B star thus also reducing the BH
mass, Abdul-Masih et al. (2020a) and El-Badry & Quataert (2020) showed that
the wobbling Hα wings could also be explained by a superposed broad Hα
absorption component on a static emission profile, thus questioning the BH
hypothesis as a whole. Shenar et al. (2020) disentangled the optical spectra
of LB-1 and proposed a binary system consisting of a stripped He-star (the
primary) in a $79\,{\rm d}$ orbit with a rapidly-rotating Be-star, presumably
the product of a recent mass-transfer event. Alternatively, Rivinius et al.
(2020) proposed that the Be-star is a tertiary object, leaving the secondary
component unidentified. Recently, Lennon et al. (2021) modelled the spectral
energy distribution (SED) of LB-1 in the optical and ultra-violet regime using
1D plane-parallel model-atmospheres, and compared several models (including
the aforementioned stripped-star/Be-star and B-star/BH hypotheses) with
observations from the Hubble Space Telescope. While their models slightly
favour the B-star and BH scenario, the stripped-star and Be-star hypothesis
could not be ruled out entirely. Throughout this paper, these two competing
scenarios will be abbreviated by the B+BH and the strB+Be scenario,
respectively.

In this paper, we aim to provide further constraints on LB-1 by applying the
BOSS-3D code to models that represent the main hypotheses detailed above,
namely the B+BH and the strB+Be binary scenarios. By calculating the
corresponding Hα line profiles, we show that both hypotheses can explain the
observed phase-dependent Hα line profile provided the LB-1 system contains an
additional disc attached to the B-star primary (i.e. either to the B star in
the B+BH scenario, or to the stripped star in the strB+Be scenario). The paper
is structured as follows. In Sect. 2 we review the basic numerical techniques
for spectral synthesis of single objects, in Sect. 3 we extend these methods
to binary systems, and in Sect. 4 we apply the developed code to the LB-1
system. Finally, we summarize our findings in Sect. 5.

##  2 Single-star spectral synthesis

Figure 1: Geometry used for the single-star algorithm (top panels) and for the
binary code (bottom panels). Top left panel: Cylindrical coordinate system
$\Sigma_{\rm cyc}^{(\rm{ray})}$ (blue) with corresponding Cartesian reference
frame $\Sigma_{\rm cac}^{(\rm{ray})}$ (red) for an observer’s direction
$\boldsymbol{n}\parallel\boldsymbol{e}_{\tilde{z}}$ used to obtain the
synthetic spectrum of a given (single) object described in a spherical
coordinate system $\Sigma_{\rm spc}^{(\rm{obj})}$ (magenta) with corresponding
Cartesian reference frame $\Sigma_{\rm cac}^{(\rm{obj})}$ (black). Top middle
panel: The $p-\zeta$ (or $\tilde{x}-\tilde{y}$-plane) of the cylindrical
coordinate system (perpendicular to the observer’s direction). Top right
panel: An arbitrary $p-\tilde{z}$-slice, with $\boldsymbol{e}_{\tilde{z}}$
pointing towards the observer. The grey circles indicate the radial grid of
the spherical coordinate system of a considered (single) object, and the black
arrows correspond to the individual rays for each impact parameter $p$. The
radiative transfer is then performed along each ray with corresponding
discretized $\tilde{z}_{{k}}$ coordinates indicated by the red dots. Bottom
left panel: Two local coordinate systems $\Sigma_{1}$ and $\Sigma_{2}$
describing two individual objects indicated by the cyan and magenta quarter-
circles, respectively. Both coordinate systems are embedded within a global
coordinate system $\Sigma_{0}$ (typically but not necessarily chosen to be the
centre of mass of both objects). The vectors $\boldsymbol{t}_{1}$ and
$\boldsymbol{t}_{2}$ describe the origin of the local coordinate systems in
the global coordinate system. Bottom middle panel: As top middle panel,
however showing the triangulation of the $\tilde{x}-\tilde{y}$-plane for an
arbitrary orbital configuration in a global cylindrical coordinate system
within the binary algorithm. Bottom right panel: As top right panel, but now
for the binary system. The individual $\tilde{z}$-rays propagating through the
system of object 1 (blue) and object 2 (red) potentially need to be merged.

In this section, we briefly summarize the basic solution method for obtaining
synthetic spectra from a given single-star model in 3D. Following, for
instance, Busche & Hillier (2005), the radiation flux at frequency $\nu$
emerging from any unresolved emitting object as seen by an observer at
distance $d$ can be formulated as

$F_{\nu}=\dfrac{1}{d^{2}}\int_{0}^{2\pi}\int_{0}^{R_{\rm
max}}I_{\nu}\left(p,\zeta,\tilde{z}=R_{\rm max}\right)p{\rm d}p{\rm
d}\zeta\,,$ (1)

where $R_{\rm max}$ is the maximum size of the emitting object,
$\left(p,\zeta\right)$ are polar coordinates describing the projected disc of
the emitting object perpendicular to the observer’s direction
$\boldsymbol{n}$, and $I_{\nu}$ is the emergent specific intensity into that
direction evaluated at a distance $\tilde{z}=R_{\rm max}$.

The flux can then be obtained by numerically integrating the emergent specific
intensity over the projected disc, with the intensity to be calculated by
solving the equation of radiative transfer along all rays in the discretized
$(p_{i},\zeta_{j})$ domain of a cylindrical coordinate system $\Sigma_{\rm
cyc}^{({\rm ray})}=\left(p,\zeta,\tilde{z}\right)$ with corresponding
Cartesian reference frame $\Sigma_{\rm cac}^{({\rm
ray})}=\left(\tilde{x},\tilde{y},\tilde{z}\right)$, and with the
$\tilde{z}$-axis being aligned with the direction to the observer (see also
Fig. 1). The corresponding (time-independent) equation of radiative transfer
reads:

$\dfrac{{\rm d}I_{\nu}\left(p_{i},\zeta_{j},\tilde{z}\right)}{{\rm
d}\tilde{z}}=\chi_{\nu}\left(p_{i},\zeta_{j},\tilde{z}\right)\left[S_{\nu}\left(p_{i},\zeta_{j},\tilde{z}\right)-I_{\nu}\left(p_{i},\zeta_{j},\tilde{z}\right)\right]\,,$
(2)

with the opacity $\chi_{\nu}$ and source function $S_{\nu}$ describing the
absorption and emission of all continuum and line processes. For simplicity,
we neglect continuum processes in the following and only consider a single
line transition between lower and upper level $l\leftrightarrow u$ with
transition frequency $\nu_{lu}$. Further, we omit the notation for explicit
spatial dependencies. Then, the opacity and source functions are given by:

$\displaystyle\chi_{\nu}$ $\displaystyle=$ $\displaystyle\dfrac{\pi{\rm
e}^{2}}{m_{e}c}\left(gf\right)\left[\dfrac{n_{l}}{g_{l}}-\dfrac{n_{u}}{g_{u}}\right]\Phi_{lu}\left(\nu\right)$
(3) $\displaystyle S_{\nu}$ $\displaystyle=$
$\displaystyle\dfrac{2h\nu_{lu}^{3}}{c^{2}}\dfrac{1}{\dfrac{n_{l}}{n_{u}}\dfrac{g_{l}}{g_{u}}-1}\,,$
(4)

where $n_{l}$ and $n_{u}$ are the occupation numbers, $g_{l}$ and $g_{u}$ are
the statistical weights, and

$\Phi_{lu}\left(\nu\right)=\dfrac{1}{\sqrt{\pi}\Delta\nu_{\rm
D}}\exp\left[-\left(\dfrac{\nu-\nu_{lu}}{\Delta\nu_{\rm
D}}-\dfrac{\boldsymbol{n}\cdot\boldsymbol{v}}{v_{\rm th}}\right)^{2}\right]$
(5)

is the profile function (approximated here by a Doppler profile) at observers
frame frequency $\nu$ including non-relativistic Doppler shifts from the
observer’s to the comoving frame via the projected velocity
$\boldsymbol{n}\cdot\boldsymbol{v}$. The width of the profile, $\Delta\nu_{\rm
D}=\nu_{lu}\dfrac{v_{\rm th}}{c}$, is given from the thermal velocity of the
considered species with mass number $m_{\rm A}$, accounting for the local
temperature $T$ and a micro-turbulent velocity $v_{\rm micro}$:

$v_{\rm th}=\sqrt{\dfrac{2k_{\rm B}T}{m_{\rm A}}+v_{\rm micro}^{2}}\,.$ (6)

For a given atmospheric model (with known temperatures, velocity vectors,
occupation numbers, and micro-turbulent velocities), the equation of radiative
transfer, Eq. (2), becomes an ordinary differential equation with exact
solution from any point $\tilde{z}_{k-1}$ to $\tilde{z}_{k}$

$I_{\nu}\left(\tilde{z}_{k}\right)=I_{\nu}\left(\tilde{z}_{k-1}\right){e}^{-\Delta\tau_{k}}+{e}^{-\Delta\tau_{k}}\int_{0}^{\Delta\tau_{k}}e^{t}S(t){\rm
d}t\,.$ (7)

We note that
$\Delta\tau_{k}:=\int_{\tilde{z}_{k-1}}^{\tilde{z}_{k}}\chi_{\nu}{\rm
d}\tilde{z}$ is the optical-depth increment between two positions along a
particular ray. Equation (7) is discretized by approximating the source
function by a constant, linear, or quadratic form in optical-depth space (e.g.
Hennicker et al. 2020, Eq. 12) based on Kunasz & Auer (1988, linear and
quadratic interpolations) and Hayek et al. (2010, monotonic Bézier
interpolations).

Since we expect the input model to be given in spherical coordinates
$\Sigma_{\rm spc}^{({\rm obj})}=(r,\Theta,\Phi)$ with corresponding Cartesian
reference frame $\Sigma_{\rm cac}^{({\rm obj})}=(x,y,z)$, and if we assume
that all structures are well resolved by the corresponding discretized grid,
the discretization of $\tilde{z}_{k}$ is easily constructed following the
standard $pz$-geometry approach (e.g. Mihalas 1978, see also Fig. 1):

$\tilde{z}_{k}\in\begin{cases}\left[-\sqrt{r_{q}^{2}-p_{i}^{2}},\sqrt{r_{q}^{2}-p_{i}^{2}}\right]\quad\forall
r_{q}\geq\lvert p_{i}\rvert\quad{\rm and}\quad r_{q}\geq R_{\ast},\\\
\left[\sqrt{R_{\ast}^{2}-p_{i}^{2}},\sqrt{r_{q}^{2}-p_{i}^{2}}\right]\,\,\quad\forall
r_{q}\geq\lvert p_{i}\rvert\quad{\rm and}\quad r_{q}<R_{\ast}\end{cases}\,,$
(8)

where $r_{q}$ refers to the discretized grid of the spherical coordinate
system. Thus, the required coordinates for each ray corresponding to a certain
$(p_{{i}},\zeta_{{j}})$ pair are determined in the cylindrical coordinate
system $\Sigma_{\rm cyc}^{(\rm{ray})}$. All quantities describing the state of
the gas along a ray need to be interpolated from the spherical grid of the
input model. To this end, we need to transform coordinates from $\Sigma_{\rm
cyc}^{(\rm{ray})}\rightarrow\Sigma_{\rm
cac}^{(\rm{ray})}\rightarrow\Sigma_{\rm
cac}^{(\rm{obj})}\rightarrow\Sigma_{\rm spc}^{(\rm{obj})}$:

$\displaystyle\begin{pmatrix}\tilde{x}\\\ \tilde{y}\\\
\tilde{z}\end{pmatrix}=\begin{pmatrix}p\cos\zeta\\\ p\sin\zeta\\\
\tilde{z}\end{pmatrix}$ (9) $\displaystyle\begin{pmatrix}x\\\ y\\\
z\end{pmatrix}=\begin{pmatrix}\boldsymbol{e}_{\tilde{x}}\cdot\boldsymbol{e}_{x}&\boldsymbol{e}_{\tilde{y}}\cdot\boldsymbol{e}_{x}&\boldsymbol{e}_{\tilde{z}}\cdot\boldsymbol{e}_{x}\\\
\boldsymbol{e}_{\tilde{x}}\cdot\boldsymbol{e}_{y}&\boldsymbol{e}_{\tilde{y}}\cdot\boldsymbol{e}_{y}&\boldsymbol{e}_{\tilde{z}}\cdot\boldsymbol{e}_{y}\\\
\boldsymbol{e}_{\tilde{x}}\cdot\boldsymbol{e}_{z}&\boldsymbol{e}_{\tilde{y}}\cdot\boldsymbol{e}_{z}&\boldsymbol{e}_{\tilde{z}}\cdot\boldsymbol{e}_{z}\end{pmatrix}\cdot\begin{pmatrix}\tilde{x}\\\
\tilde{y}\\\ \tilde{z}\end{pmatrix}=\boldsymbol{A}\cdot\tilde{\boldsymbol{r}}$
(10) $\displaystyle\begin{pmatrix}r\\\ \Theta\\\
\Phi\end{pmatrix}=\begin{pmatrix}\sqrt{x^{2}+y^{2}+z^{2}}\\\
\tan^{-1}\left(\sqrt{x^{2}+y^{2}}/z\right)\\\
\tan^{-1}\left(y/x\right)\end{pmatrix}\,,$ (11)

with $\boldsymbol{e}_{x,y,z}$ and
$\boldsymbol{e}_{\tilde{x},\tilde{y},\tilde{z}}$ the unit vectors of the
Cartesian reference frames for the spherical and cylindrical systems,
respectively, and $\boldsymbol{A}$ the transformation matrix from the
Cartesian ray system to the object’s Cartesian reference frame. The unit
vectors $\boldsymbol{e}_{\tilde{x}}$ and $\boldsymbol{e}_{\tilde{y}}$ can be
chosen arbitrarily under the constraint that
$\boldsymbol{e}_{\tilde{z}}\parallel\boldsymbol{n}$ and
$\boldsymbol{e}_{\tilde{x}}\perp\boldsymbol{e}_{\tilde{y}}\perp\boldsymbol{e}_{\tilde{z}}$.
With the obtained coordinates in the spherical system of the object, all
quantities on the discretized ray are calculated by tri-linear interpolations.
To ensure that the line-profile function along a ray is resolved by the
$\tilde{z}$-grid, we refine the $\tilde{z}$-grid if the projected velocity
steps are above $v_{\rm th}/3$.

For a single object, synthetic line profiles are thus obtained for a specific
observer’s direction by performing the following steps:

  1. 1.

Definition of the $(p,\zeta)$ grid, and of the transformation matrix relating
the Cartesian reference frames of the cylindrical (ray) and spherical (object)
coordinate systems.

  2. 2.

Calculation of the $\tilde{z}_{{k}}$-grid for each $(p_{{i}},\zeta_{{j}})$
following Eq. (8).

  3. 3.

Transformation of the coordinates $(p_{{i}},\zeta_{{j}},\tilde{z}_{{k}})$ to
the object’s spherical system by Eqs. (9)-(11), and interpolation of all
required quantities (i.e. opacities, source functions, and velocity
components).

  4. 4.

Solving the discretized form of the radiative transfer equation, Eq. (7),
along each ray to obtain the emerging specific intensity. Boundary conditions
are specified by photospheric line profiles for the specific intensity if a
ray hits the stellar core (implemented via a user-specified subroutine to
couple either line profiles as obtained from Kurucz atmosphere models, e.g.
Castelli & Kurucz 2003, or from FASTWIND)111Alternatively, one could also
explicitly solve the radiative transfer through a pre-specified density and
temperature stratification in the photosphere, of course.. For non-core rays,
the incident intensity at the outer boundary is set to zero.

  5. 5.

Integration of emerging specific intensities over the $(p,\zeta)$-plane by
replacing the integral in Eq. (1) with a discretized sum, in order to obtain
the flux at a specific frequency bin.

  6. 6.

Restarting the procedure at step 2 for the next frequency bin.

##  3 Multi-object spectral synthesis

In this section, we extend the solution scheme for a single object described
above to multiple systems by accounting for the appropriate coordinate
transformations and triangulation techniques. When accounting for a secondary
object, the algorithm above essentially suffers from one key issue, namely the
resolution of both involved objects. For instance, the spherical grid of the
primary object would not resolve the secondary object in the required detail
when simply adding the secondary into the spherical grid. Although a
description in Cartesian coordinates (presumably requiring mesh-refinement
techniques) might help here, we follow a different path within this paper.

Our method is based on the simple observation that each individual object in a
multiple system can be described on an individual spherical grid tailored to
the properties (e.g. length scale) of the object (see Fig. 1). All spherical
grids are then embedded within a global coordinate system with origin chosen
at the centre of mass of the multiple system. To keep track of the coordinate
representations for each object $q$, we define the following transformations
between the global system $\Sigma_{{\rm cac},0}^{(\rm{obj})}$ and the
individual systems $\Sigma_{{\rm cac},q}^{(\rm{obj})}$ (see Appendix A):

$\displaystyle\boldsymbol{r}_{{\rm cac},0}^{(\rm{obj})}$ $\displaystyle=$
$\displaystyle\dfrac{L_{q}}{L_{0}}\boldsymbol{Q}_{q}\cdot\boldsymbol{r}_{{\rm
cac},q}^{(\rm{obj})}+\boldsymbol{t}_{q}$ (12)
$\displaystyle\boldsymbol{r}_{{\rm cac},q}^{(\rm{obj})}$ $\displaystyle=$
$\displaystyle\dfrac{L_{0}}{L_{q}}\boldsymbol{Q}_{q}^{-1}\cdot\left(\boldsymbol{r}_{{\rm
cac},0}^{(\rm{obj})}-\boldsymbol{t}_{q}\right)\,,$ (13)

where $\boldsymbol{r}_{{\rm cac},0}^{(\rm{obj})}$ and $\boldsymbol{r}_{{\rm
cac},q}^{(\rm{obj})}$ describe the coordinates of a given point in the global
and individual systems, respectively, $L_{0}$ and $L_{q}$ are the length
scales of the coordinate systems, $\boldsymbol{t}_{q}$ is the translation
vector to the origin of the individual systems measured in the global system,
and the transformation matrices $\boldsymbol{Q}_{q}$ are given by:

$\boldsymbol{Q}_{q}=\begin{pmatrix}\boldsymbol{e}_{x}^{(q)}\cdot\boldsymbol{e}_{x}^{(0)}&\boldsymbol{e}_{y}^{(q)}\cdot\boldsymbol{e}_{x}^{(0)}&\boldsymbol{e}_{z}^{(q)}\cdot\boldsymbol{e}_{x}^{(0)}\\\
\boldsymbol{e}_{x}^{(q)}\cdot\boldsymbol{e}_{y}^{(0)}&\boldsymbol{e}_{y}^{(q)}\cdot\boldsymbol{e}_{y}^{(0)}&\boldsymbol{e}_{z}^{(q)}\cdot\boldsymbol{e}_{y}^{(0)}\\\
\boldsymbol{e}_{x}^{(q)}\cdot\boldsymbol{e}_{z}^{(0)}&\boldsymbol{e}_{y}^{(q)}\cdot\boldsymbol{e}_{z}^{(0)}&\boldsymbol{e}_{z}^{(q)}\cdot\boldsymbol{e}_{z}^{(0)}\end{pmatrix}\,.$
(14)

The unit vectors $\boldsymbol{e}_{x,y,z}^{(q)}$ of the individual coordinate
systems $\Sigma_{\rm{cac},q}^{(\rm{obj})}$ are described in the global system
$\Sigma_{\rm{cac},0}^{(\rm{obj})}$, and allow for tilts of the individual
objects within the global coordinate system. In regions where individual
coordinate systems are overlapping, the local state of the gas needs to be
described consistently within all grids.

For a given observer’s direction $\boldsymbol{n}_{0}$ measured in the global
frame then, we compute the corresponding representations in the individual
coordinate systems,
$\boldsymbol{n}_{q}=\boldsymbol{Q}_{q}^{-1}\cdot\boldsymbol{n}_{0}$, and
assign cylindrical coordinate systems to each of the objects following the
procedure outlined in Sect. 2 by setting
$\boldsymbol{e}_{\tilde{z}}^{(q)}=\boldsymbol{n}_{q}$. Again, to keep track of
the different coordinate systems, we define the transformation from the local
ray-coordinates to the local object coordinates:

$\boldsymbol{r}_{{\rm
cac},q}^{(\rm{obj})}=\boldsymbol{A}_{q}\cdot\boldsymbol{r}_{{\rm
cac},q}^{(\rm{ray})}\qquad\Longleftrightarrow\qquad\boldsymbol{r}_{{\rm
cac},q}^{(\rm{ray})}=\boldsymbol{A}_{q}^{-1}\cdot\boldsymbol{r}_{{\rm
cac},q}^{(\rm{obj})}\,,$ (15)

with transformation matrix

$\boldsymbol{A}_{q}=\begin{pmatrix}\boldsymbol{e}_{\tilde{x}}^{(\rm{q})}\cdot\boldsymbol{e}_{x}&\boldsymbol{e}_{\tilde{y}}^{(\rm{q})}\cdot\boldsymbol{e}_{x}&\boldsymbol{e}_{\tilde{z}}^{(\rm{q})}\cdot\boldsymbol{e}_{x}\\\
\boldsymbol{e}_{\tilde{x}}^{(\rm{q})}\cdot\boldsymbol{e}_{y}&\boldsymbol{e}_{\tilde{y}}^{(\rm{q})}\cdot\boldsymbol{e}_{y}&\boldsymbol{e}_{\tilde{z}}^{(\rm{q})}\cdot\boldsymbol{e}_{y}\\\
\boldsymbol{e}_{\tilde{x}}^{(\rm{q})}\cdot\boldsymbol{e}_{z}&\boldsymbol{e}_{\tilde{y}}^{(\rm{q})}\cdot\boldsymbol{e}_{z}&\boldsymbol{e}_{\tilde{z}}^{(\rm{q})}\cdot\boldsymbol{e}_{z}\end{pmatrix}\,,$
(16)

and $\boldsymbol{e}_{x}=(1,0,0)$, $\boldsymbol{e}_{y}=(0,1,0)$,
$\boldsymbol{e}_{z}=(0,0,1)$ the representation of the object’s local
coordinate system within the corresponding basis. Accordingly, we also define
a global ‘cylindrical’ coordinate system with
$\boldsymbol{e}_{\tilde{z}}=\boldsymbol{n}_{0}$, since the individual
cylindrical coordinate systems may overlap (e.g. during transits). The flux
integration, Eq. (1), is then performed in the global cylindrical system, and
needs to be adapted to avoid double counting of overlapping areas. To this
end, we consider the global 2D Delaunay-triangulation of all
$(p_{{i}},\zeta_{{j}})_{q}$ points from all individual cylindrical coordinate
systems by using the GEOMPACK2
library222https://people.math.sc.edu/Burkardt/f_src/geompack2/geompack2.html
(Joe 1991). The flux integration, Eq. (1), can then be performed by
numerically integrating the intensities over all triangles.

With Eqs. (12), (13), and (15) representing the formalism to transform
coordinates between the individual and global ray and object systems, the
basic algorithm is then defined as follows:

  1. 1.

Definition of the individual coordinate systems for all objects $q$, that is
length-scales $L_{q}$, translation vectors $\boldsymbol{t}_{q}$, basis vectors
(described in the global coordinate system)
$\boldsymbol{e}_{x,y,z}^{(\rm{q})}$, global velocities
$\boldsymbol{v}^{(\rm{q})}$ (corresponding to the orbital velocities of the
objects), and calculation of the transformation matrices $\boldsymbol{Q}_{q}$
via Eq. (14).

  2. 2.

Definition of the $(p,\zeta)_{q}$ grid for all objects $q$, and calculation of
the transformation matrices between ray and object coordinate systems
$\boldsymbol{A}_{q}$ via Eq. (16).

  3. 3.

Transformation of all $(p,\zeta)_{q}$ coordinates to the global ray-coordinate
system using Eqs. (12) and (15) to obtain the corresponding
$(\tilde{x},\tilde{y})$ coordinates in the global ray coordinate system.

  4. 4.

Calculation of the 2D triangulation in the global ray system.

  5. 5.

Back-transformation of each $(\tilde{x},\tilde{y})$ coordinate to all
individual coordinate systems via Eqs. (13) and (15).

  6. 6.

Step 2 of the single-star algorithm: Setting the $\tilde{z}_{q}$ grid along a
ray within each individual coordinate system following Eq. (8).

  7. 7.

Step 3 of the single-star algorithm: Transformation of all
$(\tilde{x}_{{i}},\tilde{y}_{{j}},\tilde{z}_{{k}})_{q}$ coordinates to the
individual spherical systems and interpolation of all required physical
quantities (i.e. opacities, source functions, velocity fields) onto the ray.

  8. 8.

Transformation of the individual rays to the global ray system via Eqs. (12)
and (15).

  9. 9.

Merging all individual rays to one single ray for a given
$(\tilde{x},\tilde{y})$ point (see Fig. 1).

  10. 10.

Step 4 of the single-star algorithm: Solving the discretized form of the
radiative transfer equation, Eq. (7), along the ray.

  11. 11.

Step 5 of the single-star algorithm: Integration of emerging specific
intensities over the triangulated area in Eq. (1) to obtain the flux at a
specific frequency bin.

  12. 12.

Step 6 of the single-star algorithm: Restarting the procedure at step 6 for
the next frequency bin.

With this algorithm, we are able to calculate synthetic line profiles for any
given orbital configuration and any given state of the surrounding material.
We emphasize that this procedure is independent of the length scales of the
individual systems, and automatically accounts for transits. Thus, one could
also analyse, for example, planetary atmospheres or complex jet configurations
as frequently found in post-asymptotic giant branch binary systems (e.g.
Bollen et al. 2020). Moreover, this algorithm can be extended to multiple
systems, with the computation time typically scaling linearly with the number
of objects, and a quadratic scaling found only for complex situations (e.g.
when all involved objects are transiting).

##  4 LB-1

With the above described algorithm, we can readily tackle binary systems in 3D
(see also Appendix B for some test calculations of the code when applied to
single-star models). In this section we model the Hα line of the binary (or
multiple) system LB-1 as a first application, focussing on the strB+Be
scenario by Shenar et al. (2020) and on the original B+BH scenario (though
with a revised mass for the B star by Simón-Díaz et al. 2020). For both
hypotheses, we adopted the stellar parameters as obtained by these authors.
For the B star in the B+BH scenario, Simón-Díaz et al. (2020) applied the 1D
spherically symmetric NLTE code FASTWIND to deduce the stellar parameters.
Similarly, Shenar et al. (2020) applied the 1D spherically-symmetric NLTE code
PoWR (Hamann & Gräfener 2003, Sander et al. 2017) to the disentangled spectrum
of the Be-star in the strB+Be scenario, while relying on the Grid Search in
Stellar Parameters tool (GSSP, Tkachenko 2015) with synthetic spectra derived
from 1D plane-parallel LTE modelling using the synthv radiative transfer code
(Tsymbal 1996) and a grid of LLmodel atmospheres (Shulyak et al. 2004) when
analysing the stripped star. The updated stellar parameters by Lennon et al.
(2021) using the 1D plane-parallel NLTE code TLUSTY (Hubeny 1988, Hubeny &
Lanz 1995) are typically in a similar range (see Table 1). These studies
relied on investigating photospheric lines that are (at most) very weakly
contaminated by the disc (in contrast to Hα, for instance), so that the
adopted stellar parameters should be good representatives of the underlying
stars. With the given stellar parameters then, we ‘only’ need to define a
suitable disc model for the BH disc and the Be-star disc to calculate
synthetic Hα line profiles.

To this end, we have set up an axisymmetric analytical model for the disc
following Kee et al. (2018) (see also Carciofi & Bjorkman 2006) which is based
on a prescribed density stratification in the equatorial plane by a
parameterised power-law, and assumes hydrostatic equilibrium in the vertical
direction. Then, the density in the disc is:

$\rho\left(r,\Theta\right)=\rho_{0}\left(\dfrac{r\sin\Theta}{R_{\ast}}\right)^{-\beta_{\rm
D}}\exp\left[-\dfrac{GM_{\ast}}{v_{\rm
s}^{2}}\left(\dfrac{1}{r\sin\Theta}-\dfrac{1}{r}\right)\right]\,,$ (17)

with $r,\Theta$ the radial and co-latitudinal coordinates, $M_{\ast}$ the
stellar or BH mass, $R_{\ast}$ the stellar radius or the minimum radius of the
BH disc, and $\rho_{0}$ and $\beta_{\rm D}$ the base density and power-law
index for the density stratification. When $\beta_{\rm D}=15/8$, this
formulation is equivalent to the standard $\alpha$-accretion disc model by
Shakura & Sunyaev (1973) at large distances from the BH, with the
$\alpha$-prescription and accretion rate hidden333Following the formulation by
Frank et al. (2002), the base density can be derived from their Eq. 5.49 at
large distances from the accreting object, yielding $\rho_{0}=3.1\cdot
10^{-8}\alpha^{-7/10}\dot{M}_{16}^{11/20}m_{1}^{5/8}R_{\ast,10}^{-15/8}$, with
$\dot{M}_{16}$ the mass-accretion rate in $10^{16}\,{\rm g/s}$, $m_{1}$ the
mass of the central object in $M_{\odot}$, and $R_{\ast,10}$ the inner disc
radius in units of $10^{10}\,{\rm cm}$. within the parameter $\rho_{o}$. As
such, $\rho_{0}$ is related to the mass flux through the disc, and $\beta_{D}$
can be considered as a combined proxy for the disc viscosity, scale height,
and temperature as a function of radius (within both a BH accretion or a Be-
star decretion disc). For Be-stars, typical values are of the order of
$\beta_{\rm D}\in[1.5,4]$ and $\rho_{0}\in 5\cdot[10^{-10},10^{-13}]\,{\rm
g}\,{\rm cm}^{-3}$ (based on observations and Hα line fitting, e.g. Silaj et
al. 2010, see also the review by Rivinius et al. 2013). With the sound speed,
$v_{\rm s}=\sqrt{k_{\rm B}T/\mu m_{p}}$, evaluated here for a typical
temperature $T=10\,{\rm kK}$ and a mean molecular weight $\mu=0.6$, the disc
density is specified for given stellar parameters and the input parameters
$\rho_{0}$, $\beta_{\rm D}$.

Further, we assumed Keplerian rotation of the disc in the orbital plane of the
binary system, with the azimuthal component of the velocity field given as

$v_{\Phi}\left(r,\Theta\right)=\sqrt{\dfrac{GM_{\ast}}{r\sin\Theta}}\,.$ (18)

Since the temperature stratification crucially depends on the type of the disc
(BH-accretion disc or Be-star disc), the formulation corresponding to each
individual model will be described in Sects. 4.1 and 4.2, respectively. To
calculate Hα opacities and source functions, we applied Eqs. (3) and (4), with
occupation numbers calculated from Saha-Boltzmann statistics assuming full
ionization.

For such disc models, an overarching goal is to qualitatively reproduce the
observed Hα line profiles (see Fig. 2), together with the corresponding
dynamical spectrum and the radial velocity curves with a semi-amplitude for
the primary and secondary of $K_{1}\approx 53\,{\rm km}\,{\rm s}^{-1}$ and
$K_{2}\approx[6,13]\,{\rm km}\,{\rm s}^{-1}$, respectively. For the secondary
object, the interpretation of the radial velocity curve crucially depends on
the applied method, with essentially two choices at hand. Firstly, the semi-
amplitude can be calculated from the barycentric method, where the Hα line
centre at each phase is determined by considering the line wings up to a
certain flux level. Using this method, Liu et al. (2019) estimated
$K_{2}^{\rm(H_{\alpha})}\approx 6\,{\rm km}\,{\rm s}^{-1}$. However, as shown
by Abdul-Masih et al. (2020a), an absorption component of the primary object
could increase the apparent motion of the line wings, thus overestimating the
true orbital velocities. Following this argumentation, an emission component
of the primary object would underestimate the true orbital velocities when
measured with the barycentric method. Alternatively, the orbital solution
could be estimated by applying spectral disentangling techniques, with
corresponding results expected to represent the actual orbital motion of the
objects. Using this method, Shenar et al. (2020) measured
$K_{2}^{\rm(true)}\approx 11\,{\rm km}\,{\rm s}^{-1}$ for LB-1, in good
agreement with a detailed investigation of near-infrared emission lines by Liu
et al. (2020) estimating $K_{2}^{\rm(true)}\in[8,13]\,{\rm km}\,{\rm s}^{-1}$.
The differences between the $K_{2}$ values as obtained from the barycentric
method and the disentangled spectra might then be explained by an emission
component attached to the primary object. Our modelling then should ideally
reproduce the $K_{2}^{\rm(H_{\alpha})}\approx 6\,{\rm km}\,{\rm s}^{-1}$ semi-
amplitude when applying the barycentric method to the synthetic Hα lines,
while also yielding $K_{2}^{\rm(true)}\approx 11\,{\rm km}\,{\rm s}^{-1}$ from
the actual orbit of the system.

Since the eccentricity is small ($e=0.03\pm 0.01$, see Liu et al. 2019), we
assumed a circular orbit. For given masses of the individual objects, $M_{1}$,
$M_{2}$, and an observed orbital period $P\approx 79\,{\rm d}$ then, the orbit
is obtained from the two-body problem:

$\displaystyle a$ $\displaystyle=$
$\displaystyle\left[\dfrac{G\left(M_{1}+M_{2}\right)P^{2}}{4\pi^{2}}\right]^{1/3}$
(19) $\displaystyle a_{1}$ $\displaystyle=$
$\displaystyle\dfrac{M_{2}}{M_{1}+M_{2}}a\qquad
v_{1}=\dfrac{M_{2}}{M_{1}+M_{2}}\dfrac{2\pi a}{P}$ (20) $\displaystyle a_{2}$
$\displaystyle=$ $\displaystyle\dfrac{M_{1}}{M_{1}+M_{2}}a\qquad
v_{2}=\dfrac{M_{1}}{M_{1}+M_{2}}\dfrac{2\pi a}{P}\,,$ (21)

with $a$ the binary separation, $a_{1}$, $a_{2}$ the distance of the objects
to the centre of mass, and $v_{1}$ and $v_{2}$ the absolute orbital
velocities. Thus, for the observed radial velocity curve(s) and period, and
assuming the individual masses of the objects to be given, the inclination $i$
is fixed within our calculations in order to reproduce the semi-amplitude of
the radial velocity curve for the primary object by $v_{1}\cdot\sin(i)=K_{1}$.
For our best models, the orbital parameters are summarized in Table 2.

###  4.1 B-star and BH-disc (B+BH) scenario

Figure 2: Hα line profiles at different phases $\varphi$ (top row) with
corresponding mean-subtracted dynamical spectra (middle row) and radial
velocity curves (bottom row) of the BH or strB companion. The $x_{\rm
obs}$-axis describes the frequency shift from line centre expressed in
velocity space. The first three columns to the left display the solution for
the three different models, PB_SBH (primary: B star, secondary: BH), PST_SBE01
(primary: stripped star, secondary: Be-star, model 1), PST_SBE02 (primary:
stripped star, secondary: Be-star, model 2), with Hα line profiles at two
distinct phases from the actual observations indicated by the black dashed and
dotted lines. The right column additionally displays the corresponding
observed (see Shenar et al. 2020) Hα line profiles, mean-subtracted dynamical
spectrum and radial velocity curves of LB-1, where we distinguish between the
true ones (e.g. from using disentangling methods) indicated in blue and those
obtained from the barycentric method applied to the Hα line wings (red). To
reduce the noise, all observed line profiles have been convolved with a
Gaussian filter of width $10\,{\rm km}\,{\rm s}^{-1}$. The bottom panels
display the radial velocity curves of the secondary when calculated from the
barycentric method applied to the Hα line profile wings (red crosses) and when
using the true orbital parameters (blue crosses). Additionally, we display the
corresponding radial velocity curves as found by Shenar et al. (2020) and Liu
et al. (2019) using disentangling techniques (blue solid line) and the
barycentric method (red solid lines), respectively. At the top, the (phase-
averaged) equivalent width evaluated in velocity space in units of $100\,{\rm
km}\,{\rm s}^{-1}$ and the inclination used within the synthetic-spectra
calculations are indicated.

Figure 3: Disc models in the equatorial plane corresponding to our best-fit
parameters of the LB-1 system (see Table 3). The top and middle panels display
the radial density and temperature stratification, respectively, and the
azimuthal velocity component is shown in the bottom panel. The blue and red
lines indicate the BH accretion disc (model PBH_SB) and the Be-star disc
within the single-disc model (model PST_SBE01), respectively. For model
PST_SBE02 (containing two discs), the Be-star disc is indicated by the green
dotted line (partially overlapping with the red solid line from model
PST_SBE01), and the stripped-star disc is indicated by the green solid line.
Table 1: Stellar parameters for the B+BH and strB+Be scenario used for our
simulations of the LB-1 system, and as found in the literature.

|  | B+BH |  | strB+Be
---|---|---|---|---
| reference | this work | (1) | (2) |  | this work | (3) | (2)
primary | type | B star | B star | B star |  | stripped star | stripped star | stripped star
$T_{\rm eff}^{\rm(1)}\,[{\rm kK}]$ | 14 | $\in\left[13.5,14.5\right]$ | $\in\left[15.0,15.6\right]$ |  | 12.7 | $12.7$ | $\in\left[12.4,12.6\right]$
$R_{\ast}^{(1)}\,[R_{\odot}]$ | 5.3 | $\in\left[4.1,6.5\right]$ | $\in\left[5.4,6.7\right]$ |  | 5.4 | $5.4$ | $\in\left[4.4,5.3\right]$
$M_{\ast}^{(1)}\,[M_{\odot}]$ | 4.2 | $\in\left[1.3,5.5\right]$ | $\in\left[3.8,7.0\right]$ |  | 1.3 | $\in\left[1.1,1.9\right]$ | $\in\left[0.5,1.3\right]$
$\log g^{(1)}$ | 3.5 | $\in\left[3.35,3.65\right]$ | $\in\left[3.4,3.8\right]$ |  | 3.1 | $3.0$ | $\in\left[2.8,3.2\right]$
$v_{\rm rot}^{(1)}\,[{\rm km}\,{\rm s}^{-1}]$ | 21 | $<27^{\rm{(b)}}$ | $21^{\rm{(c)}}$ |  | 11 | $\in[8,14]$ | $11$
$Y_{\rm He}^{(1)}$ | 0.1 | $0.1$ | $0.1$ |  | 0.2 | $0.21$ | $0.2$
secondary | type | BH | - | BH |  | Be star | Be star | Be star
$T_{\rm eff}^{(2)}\,[{\rm kK}]$ | - | - | - |  | 18 | $18$ | $\in\left[18.7,19.1\right]$
$R_{\ast}^{(2)}\,[R_{\odot}]$ | 1 | - | - |  | 3.7 | $3.7$ | $\in\left[2.8,3.4\right]$
$M_{\ast}^{(2)}\,[M_{\odot}]$ | 30 | - | $\in\left[13,30\right]$ |  | 7 | $\in\left[5,9\right]$ | $\in\left[1.6,6.9\right]$
$\log g^{(2)}$ | - | - | - |  | 4.15 | $4.0$ | $\in\left[3.7,4.3\right]$
$v_{\rm rot}^{(2)}\,[{\rm km}\,{\rm s}^{-1}]$ | - | - | - |  | 600 | $\in[400,560]$ | $480$
$Y_{\rm He}^{(2)}$ | - | - | - |  | 0.1 | $0.08$ | $0.1$

444The values for the rotational velocities have been deduced from the
inclination used in this paper (see Table 2, $i=22^{\circ}$ and $i=39^{\circ}$
for the B+BH and strB+Be scenarios, respectively). For consistency, we
calculated $\log g$ from the stellar radii and masses.

(1) Simón-Díaz et al. (2020); (2) Lennon et al. (2021); (3) Shenar et al.
(2020).

Table 2: Orbital parameters for the LB-1 system as calculated from our models
and found in the literature.

| models |  | literature
---|---|---|---
| PB_SBH | PB_SBE01 | PB_SBE02 |  | (1) | (2) | (3)
$i\,[^{\circ}]$ | 22 | 39 | 39 |  | - | - | $\in[35,43]$
$K_{1}\,[{\rm km}\,{\rm s}^{-1}]$ | 53 | 53 | 53 |  | 53 | 53 | 53
$K_{2}^{\rm(true)}\,[{\rm km}\,{\rm s}^{-1}]$ | 7.4 | 9.9 | 9.9 |  | - | $\in[8,13]$ | $\in[10.2,12.2]$
$K_{2}^{\rm(H_{\alpha})}\,[{\rm km}\,{\rm s}^{-1}]$ | 10 | 13 | 8 |  | $\in[5.6,7.2]$ | - | -

(1) Liu et al. (2019); (2) Liu et al. (2020); (3) Shenar et al. (2020).

Table 3: Disc parameters for the three different (best-fit) models as obtained from our simulations of the LB-1 system (see text). | model | PB_SBH | PST_SBE01 | PST_SBE02
---|---|---|---|---
primary | type | B star | stripped star | stripped star
$v_{\rm micro}^{(1)}\,[{\rm km}\,{\rm s}^{-1}]$ | - | - | $30$
$\rho_{0}^{\rm(1)}\,[{\rm g}\,{\rm cm}^{-3}]$ | - | - | $1.3\cdot 10^{-12}$
$\beta_{\rm D}$ | - | - | $7/2$
$T_{\rm disc}\,[{\rm kK}]$ | - | - | Eq. (22)
secondary | type | BH | Be star | Be star
$v_{\rm micro}^{(2)}\,[{\rm km}\,{\rm s}^{-1}]$ | $100$ | $100$ | $100$
$\rho_{0}^{\rm(2)}\,[{\rm g}\,{\rm cm}^{-3}]$ | $3\cdot 10^{-12}$ | $6.9\cdot 10^{-13}$ | $6.4\cdot 10^{-13}$
$\beta_{\rm D}$ | $3/2$ | $3/2$ | $3/2$
$T_{\rm disc}\,[{\rm kK}]$ | $=\rm const.=8$ | Eq. (22) | Eq. (22)

In the following, we describe our model for the B+BH scenario (model PB_SBH,
primary: B star, secondary: BH). For the B star, we used the stellar
parameters as summarized in Table 1. With $K_{1}$ and $K_{2}^{\rm(true)}$
given, the mass of the black hole is constrained via Eqs. (20) and (21)
yielding $M_{2}=M_{1}K_{1}/K_{2}^{\rm(true)}\approx[20,30]\,M_{\odot}$, where
we have chosen the black-hole mass as $M_{2}=30\,M_{\odot}$ and inclination
$i=22^{\circ}$ in order to reproduce the semi-amplitude of the radial velocity
$K_{1}$. With this inclination, we adopted a rotational velocity for the B
star, $v_{\rm rot}^{(1)}=21\,{\rm km}\,{\rm s}^{-1}$ (with the observed value
$v{\thinspace}\sin{\thinspace}i=8\,{\rm km}\,{\rm s}^{-1}$, Lennon et al.
2021). To model the disc, we have set the minimum radius to
$R_{\ast}^{(2)}=1\,R_{\odot}$, since any regions with smaller radii would not
contribute to the emission due to the small emitting area in any case, and
have set the maximum radius of the disc to the distance from the BH to the
first Lagrange point.

For a standard $\alpha$-disc model without irradiation from an external
source, we do not expect any emission in Hα, since the surface of the hot
regions close to the BH is very small, and the B-star companion completely
dominates the total radiation flux. In the B+BH binary system, however, the
disc will also be heated by the B-star companion. In order to avoid complex
multi-D radiation-hydrodynamic simulations, and to keep our model as simple as
possible, we introduced a constant temperature within the disc as a free
parameter. With given stellar parameters and inclination then, the model is
completely specified by four input parameters, namely the disc temperature,
$T_{\rm disc}$, the micro-turbulent velocity in the disc, $v_{\rm
micro}^{(2)}$, the base density of the disc, $\rho_{0}^{(2)}$, and the slope
of the density stratification, $\beta_{\rm D}^{(2)}$.

We then calculated Hα line profiles with the algorithm described in Sect. 3
for different models (see Appendix C.1). As an inner boundary condition for
the specific intensity emerging the B star, we used a photospheric line
profile obtained from FASTWIND. Further, we calculated the radial velocity
curve from the barycentric method considering only the wings of the Hα lines
up to 1/3 of the total height of the Hα line profile as suggested by Liu et
al. (2019). The results for our best fit (by eye) with disc temperature,
$T_{\rm disc}=8\,{\rm kK}$, micro-turbulent velocity $v_{\rm
micro}^{\rm(2)}=100\,{\rm km}\,{\rm s}^{-1}$, base density
$\rho_{0}^{(2)}=3\cdot 10^{-12}\,{\rm g}\,{\rm cm}^{-3}$, and radial density
stratification $\beta_{\rm D}=3/2$ are shown in Fig. 2, with corresponding
density, temperature, and velocity profiles in the equatorial plane shown in
Fig. 3.

The slope of the density stratification is in fairly good agreement with
accretion discs described within the $\alpha$-disc prescription. Indeed, a
shallow slope of the density stratification is also required from a modelling
point of view to obtain a centrally peaked emission profile. For a steeper
decline of the density, for instance, the outer disc regions (where the
velocities are lowest) barely contribute to the low-velocity emission due to
much smaller densities.

On the other hand, the emitting area at high (projected) velocities (i.e. near
to the inner disc radius) is quite small. Thus, a relatively high micro-
turbulent velocity is required to obtain the broad Hα emission wings,
effectively increasing the emitting area (and consequently the emerging flux)
at high velocity shifts. In particular, this also smoothens the typically
double-peaked emission feature of accretion discs even at such inclinations
($i=22^{\circ}$). High turbulent velocities, however, are frequently found
when modelling dynamical systems with simplified stationary models (which by
assumption neglect small-scale turbulent motions in the gas). As just one
example, to mimick the observed velocity dispersion when analysing
magnetically confined winds described with a steady state model, Owocki et al.
(2016) convolved their synthetic Hα line profiles with a $150\,{\rm km}\,{\rm
s}^{-1}$ Gaussian ‘macro-turbulence’ velocity. Since we have neglected
gravitational interactions of the primary object with the disc, we might thus
expect the inferred high turbulent velocities to mimick dynamical effects
within the binary system.

The qualitative behaviour of the observed Hα line profile is well reproduced,
with corresponding phase-averaged equivalent width measured in velocity space
in units of $100\,{\rm km}\,{\rm s}^{-1}$, $\bar{W}_{x}=6.81$, in good
agreement with observations ($\bar{W}_{x}^{\rm(obs)}=6.84$). We emphasize,
however, that the radial velocity curves show some discrepancies when compared
to observations (see also Table 2). While the actual orbital velocities,
$K_{2}^{(\rm true)}=v_{2}\sin\left(i\right)=7.4\,{\rm km}\,{\rm s}^{-1}$, are
slightly below the observed range ($K_{2}^{(\rm true,obs)}\in[8,13]\,{\rm
km}\,{\rm s}^{-1}$), we find a higher semi-amplitude when applying the
barycentric method ($K_{2}^{\rm(H_{\alpha})}=10\,{\rm km}\,{\rm s}^{-1}$).
Indeed, this discrepancy can be explained by the absorption profile of the B
star (see discussion in Abdul-Masih et al. 2020a), and could be solved by an
additional Hα emission component attached, for example, to the B star or to a
tertiary object. Such an additional emission component would reduce the radial
velocity amplitude as measured from the barycentric method, by slightly
modifying the Hα line wings.

More pronounced, we find significant differences between the observed and
synthetic (mean-subtracted) dynamical spectra of the Hα line. While the
observed dynamical spectrum clearly displays an anti-phase behaviour of the
line wings and the line core, the corresponding features in the synthetic
dynamical spectrum are fairly well aligned. As shown below for the strB+Be
scenario, this issue could be solved by introducing also a disc around the B
star though.

A main theme of our model regards the temperature stratification of the BH
disc, where we required $T_{\rm disc}=8\,{\rm kK}$ in order to reproduce the
observations. This temperature region can indeed be considered as a sweet
spot, since lower temperatures ($T_{\rm disc}\lesssim 5\,{\rm kK}$) would
require much higher densities to provide enough emission in Hα. The
corresponding densities would translate to a high mass-accretion rate of the
BH within the standard $\alpha$-disc prescription, thus predicting an X-ray
bright accretion disc that has not been observed. Similarly, too high
temperatures ($T_{\rm disc}\gtrsim 15\,{\rm kK}$) are not allowed either,
because we would need to increase the densities to match the equivalent width
of the observed Hα line profiles (see Appendix C.1), again yielding an X-ray
bright accretion disc. If the sweet-spot temperature as found in our model is
a reasonable substitute for detailed numerical simulations still needs to be
investigated.

###  4.2 Stripped-star and Be-star (strB+Be) scenario

For the stellar parameters of the LB-1 system within the strB+Be scenario, we
adopted the values proposed by Shenar et al. (2020), summarized again in Table
1, and yielding an inclination $i=39^{\circ}$. Motivated by these authors, we
have set the rotational velocity of the stripped star to $v_{\rm
rot}^{\rm(1)}=11$ (with $v{\thinspace}\sin{\thinspace}i=7\,{\rm km}\,{\rm
s}^{-1}$ following Shenar et al. 2020 and Lennon et al. 2021), and assumed the
Be-star to rotate at its critical velocity, $v_{\rm
rot}^{\rm(2)}=\sqrt{GM_{\ast}^{\rm(1)}/R_{\ast}^{\rm(1)}}$. For simplicity, we
neglected effects of gravity darkening.

To set the temperature stratification in the disc, we followed the qualitative
behaviour as found by detailed radiative-transfer calculations in Carciofi &
Bjorkman (2006, Sect. 5.2). In regions close to the stellar surface (where the
vertical electron-scattering optical depth $\tau_{z}(r\sin\Theta)>0.1$), we
adopted an analytical temperature stratification described by a flat blackbody
reprocessing disc (their Eq. (12)):

$T_{\rm
disc}\left(r,\Theta\right)=\dfrac{T_{0}}{\pi^{1/4}}\left[\sin^{-1}\left(\dfrac{R_{\ast}}{r\sin\Theta}\right)-\dfrac{R_{\ast}}{r\sin\Theta}\sqrt{1-\dfrac{R_{\ast}^{2}}{r^{2}\sin^{2}\Theta}}\right]^{1/4}\,,$
(22)

with $T_{0}$ the base temperature of the disc. Since the base densities of our
disc models are typically relatively low, we can neglect back-warming effects
of the photosphere by the disc, and therefore have set the base temperature to
the stellar effective temperature, $T_{0}=T_{\rm eff}$. To mimick the heating
of the disc in (vertically) optically thin regions by the star, we adopted a
constant disc temperature $T_{\rm disc}=0.6\cdot T_{\rm eff}$ where
$\tau_{z}(r\sin\Theta)<0.1$ (again, following the argumentation by Carciofi &
Bjorkman 2006).

We note already here that we additionally allowed for a disc attached also to
the stripped star (see below). With given stellar parameters, a now specified
temperature stratification, and a given inclination to reproduce the semi-
amplitude of the radial velocity curve of the primary, $K_{1}$, the model is
completely specified by six input parameters, namely the micro-turbulent
velocities, $v_{\rm micro}^{\rm(1,2)}$, the base density
$\rho_{0}^{\rm(1,2)}$, and the slope of the density stratification of both
discs, $\beta_{\rm D}^{\rm(1,2)}$. When considering only the Be-star disc, the
number of free parameters is reduced to three.

As for the B+BH scenario above, we calculated synthetic Hα line profiles with
the developed BOSS-3D code for various different models (see Appendix C.2),
with the photospheric line profiles of both stars again obtained from
FASTWIND. The results for our best two models are shown in Fig. 2, with the
corresponding density, temperature, and velocity stratification in the
equatorial plane shown in Fig. 3, and the orbital and disc parameters
summarized in Tables 2 and 3, respectively.

#### Model PST_SBE01.

This model describes a stripped star (primary object) in orbit with a Be-star
(secondary object) hosting a circumstellar disc. The obtained disc parameters,
$\beta_{\rm D}^{\rm(2)}=3/2$ and $\rho_{0}^{\rm(2)}=6.9\cdot 10^{-13}\,{\rm
g}\,{\rm cm}^{-3}$ are both at the lower end of expected values (see above and
Silaj et al. 2010), which can be explained as follows: Firstly, the inner
regions of the disc (i.e. those near to the star with highest rotational
velocities) would become optically thick in the continuum555We have
approximated the continuum by pure electron-scattering opacities and assuming
LTE to determine the source function. for higher base densities, and the broad
line wings would vanish within the continuum flux. Secondly, as described
above for the BH disc, a shallow slope of the density stratification is
required to obtain a centrally peaked emission profile. Again, we also require
a high micro-turbulent velocity to obtain the broad Hα emission wings and to
smooth out the otherwise double-peaked emission profile.

Even for our best fit (by eye) parameters, however, the central emission peak
is not well reproduced with this model. Additionally, the equivalent width
averaged over all phases ($\bar{W}_{x}=6.57$, again measured in $100\,{\rm
km}\,{\rm s}^{-1}$) is somewhat smaller than expected from observations
($\bar{W}_{x}^{\rm(obs)}=6.84$). Moreover, the obtained semi-amplitude of the
radial velocity of the secondary object, $K_{2}^{\rm(H_{\alpha})}=13\,{\rm
km}\,{\rm s}^{-1}$, shows a clear deviation from the actually observed values
($K_{2}^{\rm(H_{\alpha},obs)}\approx 6\,{\rm km}\,{\rm s}^{-1}$). As in the
B+BH scenario, the semi-amplitude is larger than the orbital velocity of the
secondary object ($v_{2}\sin(i)=9.9\,{\rm km}\,{\rm s}^{-1}$), demonstrating
again the impact of the photospheric absorption profile underlying the
stripped star in this case. Finally, the dynamical Hα line profile shows the
same behaviour as for the B+BH scenario, in contrast to observations.

#### Model PST_SBE02.

To improve the model and to increase the low-velocity emission, we
additionally implemented a circumstellar disc attached to the stripped star.
Although speculative, this disc might have formed by re-accretion from the Be-
star (e.g. Shenar et al. 2020) or from earlier mass-transfer phases.

While the best fit (by eye) gives a similar disc for the Be-star as found
above, we predict a slightly higher base density of the stripped-star disc,
$\rho_{0}^{\rm(1)}=1.3\cdot 10^{-12}\,{\rm g}\,{\rm cm}^{-3}$, and a
relatively steep slope of the density stratification, $\beta_{\rm
D}^{\rm(1)}=7/2$. Since the orbital velocities of the stripped-star disc are
quite low in any case (Eq. (18) with a small mass and large radius), this
model gives an increased central emission peak, at least for moderate micro-
turbulent velocities. The qualitative behaviour of the observed Hα line
profile then is very well reproduced, with corresponding equivalent width,
$\bar{W}_{x}^{\rm(obs)}=6.84$, in very good agreement with observations. Since
the Hα line wings are slightly contaminated by the (anti-phase) emission from
the stripped-star disc, the obtained radial velocity curve
($K_{2}^{\rm(H_{\alpha})}=8\,{\rm km}\,{\rm s}^{-1}$) agrees fairly well with
the corresponding observations ($K_{2}^{\rm(H_{\alpha},obs)}=6\,{\rm km}\,{\rm
s}^{-1}$ as derived from the observed Hα line wings). This underestimation of
the true orbital velocities due to the contamination of the line profile by an
emission component from the primary object is tightly connected to the
opposite effect as found by Abdul-Masih et al. (2020a) for an absorption
component.

While we were not aiming at an exact reproduction of the radial velocity
curves due to the simplicity of our disc model (e.g. neglecting disc
inhomogeneities), we show that the disc’s emission component can qualitatively
explain the discrepancies between the true and apparent orbital motion of the
system. We emphasize that any sort of relatively narrow, centrally peaked Hα
emission originating from the stripped star might provide similar results. The
origin of such narrow Hα emission, however, remains unclear. For instance, we
were not able to reproduce the observations with a simple ($\beta$-velocity-
type) stellar wind blown from the stripped star. In contrast, the dynamical
spectrum of our strB+Be model with an additional disc around the stripped star
is in good agreement with observations. Thus, and due to the simplicity of our
disc model (with only three free parameters for each disc), we propose that
LB-1 contains two discs, each attached to the individual objects in a strB+Be
binary system (as shown here), or in a B+BH binary (as qualitatively explained
above). We emphasize that Shenar et al. (2020) found a centrally peaked Hα
emission from the stripped star by using disentangling techniques as well,
suggesting that this residual component might originate from a disc around the
stripped star. Our present results put this interpretation on a much firmer
ground.

##  5 Summary and conclusions

In this paper, we have presented a new method for calculating synthetic
spectra of binary systems, assuming the orbital setup and the atmospheric and
circumstellar structure(s) to be known. This method can be extended to
multiple systems, with the computation time scaling only linearly (or slightly
super-linearly) with the number of involved objects. Moreover, by assigning
individual coordinate systems to each object and accounting for the
appropriate coordinate transformations, the developed method automatically
accounts for the different length-scales of each object. The presented
algorithm and associated code (BOSS-3D) is therefore capable of calculating
synthetic line profiles for a wide variety of astrophysical systems, such as
planetary or stellar transits, jets associated with young stellar objects
(YSO’s) or post-asymptotic giant branch binary systems (see also Bollen et al.
2020), and multiple systems involving (possibly colliding) stellar winds, Be
stars, or black holes (BH’s) with discs.

As a first application of the method, we considered the Hα line formation of
the B-star/BH-disc (B+BH) and stripped-star/Be-star (strB+Be) scenarios for
the enigmatic (multiple) system LB-1 that have previously been proposed (among
other hypotheses, see Sect. 1) by Liu et al. (2019) and Shenar et al. (2020),
respectively. To calculate the density stratification and velocity field, we
applied an analytical, geometrically thin, Keplerian disc model in vertical
hydrostatic equilibrium. Further, we assigned a constant temperature to the BH
disc as a free parameter, and used a prescribed temperature stratification
based on Carciofi & Bjorkman (2006) for the Be-star disc. The level
populations for the Hα line transition have then been computed from Saha-
Boltzmann statistics assuming full ionization.

Under these assumptions none of the models can reproduce the detailed phase-
dependent shape of the observed line profiles, particularly due to differences
in the dynamical spectrum (and missing low-velocity emission in some cases).
Moreover, the radial velocity curves measured from the synthetic Hα line wings
are highly overestimated due to the absorption component of the primary object
(see also Abdul-Masih et al. 2020a). It remains open whether NLTE calculations
can help to explain the detailed phase-dependent shape of Hα line profiles.
Detailed theoretical investigations are needed to test these effects.
Alternatively, we have empirically set an additional disc around the stripped
star within the strB+Be scenario to solve this problem, and found a sound
reproduction of the observed Hα line profiles both from their qualitative
shape and their equivalent widths, as well as of the dynamical spectrum and
the radial velocity curves. This putative disc might be a remnant associated
with earlier mass-transfer phases or could have originated from re-accretion
of material from the Be-star disc.

Similarly, such a disc might be attached to the B star in the B+BH scenario as
well, and the B+BH hypothesis remains a valid possibility to explain the LB-1
system. In any case, our findings provide strong evidence that LB-1 contains a
disc-disc system, with an additional disc either attached to the B star (in
the B+BH scenario) or to the stripped star (in the strB+Be scenario).

###### Acknowledgements.

We thank our referee, Dr. Maria Bergemann, for many helpful comments and
suggestions. LH and JOS gratefully acknowledge support from the Odysseus
program of the Belgian Research Foundation Flanders (FWO) under grant
G0H9218N. JB acknowledges support from the FWO Odysseus program under project
G0F8H6N. The National Solar Observatory (NSO) is operated by the Association
of Universities for Research in Astronomy, Inc. (AURA), under cooperative
agreement with the National Science Foundation. TS acknowledges funding
received from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (grant agreement number 772225:
MULTIPLES) and from the European Union’s Horizon 2020 under the Marie
Skłodowska-Curie grant agreement No 101024605. IEM has received funding from
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (SPAWN ERC, grant agreement No 863412).

## References

  * Abdul-Masih et al. (2020a) Abdul-Masih, M., Banyard, G., Bodensteiner, J., et al. 2020a, Nature, 580, E11
  * Abdul-Masih et al. (2020b) Abdul-Masih, M., Sana, H., Conroy, K. E., et al. 2020b, A&A, 636, A59
  * Belczynski et al. (2020) Belczynski, K., Hirschi, R., Kaiser, E. A., et al. 2020, ApJ, 890, 113
  * Bollen et al. (2020) Bollen, D., Kamath, D., De Marco, O., Van Winckel, H., & Wardle, M. 2020, A&A, 641, A175
  * Busche & Hillier (2005) Busche, J. R. & Hillier, D. J. 2005, AJ, 129, 454
  * Carciofi & Bjorkman (2006) Carciofi, A. C. & Bjorkman, J. E. 2006, ApJ, 639, 1081
  * Castelli & Kurucz (2003) Castelli, F. & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, Vol. 210, A20
  * El-Badry & Quataert (2020) El-Badry, K. & Quataert, E. 2020, MNRAS, 493, L22
  * Frank et al. (2002) Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition
  * Hamann & Gräfener (2003) Hamann, W. R. & Gräfener, G. 2003, A&A, 410, 993
  * Hauschildt (1992) Hauschildt, P. H. 1992, J. Quant. Spec. Radiat. Transf., 47, 433
  * Hauschildt & Baron (2006) Hauschildt, P. H. & Baron, E. 2006, A&A, 451, 273
  * Hayek et al. (2010) Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49
  * Hennicker et al. (2018) Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2018, A&A, 616, A140
  * Hennicker et al. (2020) Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2020, A&A, 633, A16
  * Hillier (2012) Hillier, D. J. 2012, in From Interacting Binaries to Exoplanets: Essential Modeling Tools, ed. M. T. Richards & I. Hubeny, Vol. 282, 229–234
  * Hillier & Miller (1998) Hillier, D. J. & Miller, D. L. 1998, ApJ, 496, 407
  * Hubeny (1988) Hubeny, I. 1988, Computer Physics Communications, 52, 103
  * Hubeny & Lanz (1995) Hubeny, I. & Lanz, T. 1995, ApJ, 439, 875
  * Joe (1991) Joe, B. 1991, Advances in Engineering Software and Workstations, 13, 325
  * Kee et al. (2018) Kee, N. D., Owocki, S., & Kuiper, R. 2018, MNRAS, 474, 847
  * Kunasz & Auer (1988) Kunasz, P. & Auer, L. H. 1988, J. Quant. Spec. Radiat. Transf., 39, 67
  * Lamers et al. (1987) Lamers, H. J. G. L. M., Cerruti-Sola, M., & Perinotto, M. 1987, ApJ, 314, 726
  * Lennon et al. (2021) Lennon, D. J., Maíz Apellániz, J., Irrgang, A., et al. 2021, A&A– accepted, arXiv:2103.14069
  * Liu et al. (2019) Liu, J., Zhang, H., Howard, A. W., et al. 2019, Nature, 575, 618
  * Liu et al. (2020) Liu, J., Zheng, Z., Soria, R., et al. 2020, ApJ, 900, 42
  * Lobel & Blomme (2008) Lobel, A. & Blomme, R. 2008, ApJ, 678, 408
  * Mihalas (1978) Mihalas, D. 1978, Stellar atmospheres (2nd edition) (San Francisco: W. H. Freeman and Co., 1978)
  * Owocki et al. (2016) Owocki, S. P., ud-Doula, A., Sundqvist, J. O., et al. 2016, MNRAS, 462, 3830
  * Pauldrach et al. (2001) Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161
  * Puls et al. (2020) Puls, J., Najarro, F., Sundqvist, J. O., & Sen, K. 2020, A&A, 642, A172
  * Rivinius et al. (2020) Rivinius, T., Baade, D., Hadrava, P., Heida, M., & Klement, R. 2020, A&A, 637, L3
  * Rivinius et al. (2013) Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&A Rev., 21, 69
  * Sander et al. (2017) Sander, A. A. C., Hamann, W. R., Todt, H., Hainich, R., & Shenar, T. 2017, A&A, 603, A86
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 500, 33
  * Shenar et al. (2020) Shenar, T., Bodensteiner, J., Abdul-Masih, M., et al. 2020, A&A, 639, L6
  * Shulyak et al. (2004) Shulyak, D., Tsymbal, V., Ryabchikova, T., Stütz, C., & Weiss, W. W. 2004, A&A, 428, 993
  * Silaj et al. (2010) Silaj, J., Jones, C. E., Tycner, C., Sigut, T. A. A., & Smith, A. D. 2010, ApJS, 187, 228
  * Simón-Díaz et al. (2020) Simón-Díaz, S., Maíz Apellániz, J., Lennon, D. J., et al. 2020, A&A, 634, L7
  * Sobolev (1960) Sobolev, V. V. 1960, Moving envelopes of stars (Cambridge: Harvard University Press, 1960)
  * Sundqvist & Puls (2018) Sundqvist, J. O. & Puls, J. 2018, A&A, 619, A59
  * Sundqvist et al. (2012) Sundqvist, J. O., ud-Doula, A., Owocki, S. P., et al. 2012, MNRAS, 423, L21
  * Tkachenko (2015) Tkachenko, A. 2015, A&A, 581, A129
  * Tsymbal (1996) Tsymbal, V. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 108, M.A.S.S., Model Atmospheres and Spectrum Synthesis, ed. S. J. Adelman, F. Kupka, & W. W. Weiss, 198

##  Appendix A Coordinate transformations

In this section, we derive the basic transformations of coordinates between
two coordinate systems (see Fig. 4) with orthogonal basis such that
$\Sigma_{1}=\left(\boldsymbol{e}_{x}^{(1)},\boldsymbol{e}_{y}^{(1)},\boldsymbol{e}_{z}^{(1)}\right)$
and
$\Sigma_{2}=\left(\boldsymbol{e}_{x}^{(2)},\boldsymbol{e}_{y}^{(2)},\boldsymbol{e}_{z}^{(2)}\right)$.
Any point $\boldsymbol{r}_{\Sigma_{1}}$ in $\Sigma_{1}$ can then be expressed
by

$\displaystyle\boldsymbol{r}_{\Sigma_{1}}$ $\displaystyle=$ $\displaystyle
x_{1}\boldsymbol{e}_{x}^{(1)}+y_{1}\boldsymbol{e}_{y}^{(1)}+z_{1}\boldsymbol{e}_{z}^{(1)}$
(23) $\displaystyle=$ $\displaystyle
s_{x}x_{2}\boldsymbol{e}_{x}^{(2)}+s_{y}y_{2}\boldsymbol{e}_{y}^{(2)}+s_{z}z_{2}\boldsymbol{e}_{z}^{(2)}+\boldsymbol{t}\,,$

where $(x_{1},y_{1},z_{1})$ and $(x_{2},y_{2},z_{2})$ are the coordinates in
each system, $s_{x,y,z}$ are scale factors allowing us to include different
length scales for both coordinate systems,
$\boldsymbol{t}=t_{x}\boldsymbol{e}_{x}^{(1)}+t_{y}\boldsymbol{e}_{y}^{(1)}+t_{z}\boldsymbol{e}_{z}^{(1)}$
is the translation vector indicating the origin of $\Sigma_{2}$ within
$\Sigma_{1}$, and $\boldsymbol{e}_{x,y,z}^{(2)}$ are the basis vectors of
$\Sigma_{2}$ measured within $\Sigma_{1}$.

Figure 4: Transformation between coordinate systems
$\Sigma_{1}=\left(\boldsymbol{e}_{x}^{(1)},\boldsymbol{e}_{y}^{(1)},\boldsymbol{e}_{z}^{(1)}\right)$
and
$\Sigma_{2}=\left(\boldsymbol{e}_{x}^{(2)},\boldsymbol{e}_{y}^{(2)},\boldsymbol{e}_{z}^{(2)}\right)$,
where $\boldsymbol{t}$ indicates the translation vector from
$\Sigma_{1}\rightarrow\Sigma_{2}$, and $\boldsymbol{r}_{\Sigma_{1}}$ and
$\boldsymbol{r}_{\Sigma_{2}}$ denote the coordinate representations in each
system for a given point.

For our purposes, the length scales of both coordinate systems are constant
(and typically refer to the radius of each individual object). Thus, we define
$s:=L_{\Sigma_{2}}/L_{\Sigma_{1}}=s_{x}=s_{y}=s_{z}$, with $L_{\Sigma_{2}}$
and $L_{\Sigma_{1}}$ describing the length scales of the individual coordinate
systems. Multiplying Eq. (23) with $\boldsymbol{e}_{x}^{(1)}$,
$\boldsymbol{e}_{y}^{(1)}$, and $\boldsymbol{e}_{z}^{(1)}$, we obtain in
matrix form:

$\displaystyle\underbrace{\begin{pmatrix}x_{1}\\\ y_{1}\\\
z_{1}\end{pmatrix}}_{\boldsymbol{r}_{\Sigma_{1}}}$ $\displaystyle=$
$\displaystyle\underbrace{\begin{pmatrix}s&0&0\\\ 0&s&0\\\
0&0&s\end{pmatrix}}_{=:\boldsymbol{S}}\underbrace{\begin{pmatrix}\boldsymbol{e}_{x}^{(2)}\cdot\boldsymbol{e}_{x}^{(1)}&\boldsymbol{e}_{y}^{(2)}\cdot\boldsymbol{e}_{x}^{(1)}&\boldsymbol{e}_{z}^{(2)}\cdot\boldsymbol{e}_{x}^{(1)}\\\
\boldsymbol{e}_{x}^{(2)}\cdot\boldsymbol{e}_{y}^{(1)}&\boldsymbol{e}_{y}^{(2)}\cdot\boldsymbol{e}_{y}^{(1)}&\boldsymbol{e}_{z}^{(2)}\cdot\boldsymbol{e}_{y}^{(1)}\\\
\boldsymbol{e}_{x}^{(2)}\cdot\boldsymbol{e}_{z}^{(1)}&\boldsymbol{e}_{y}^{(2)}\cdot\boldsymbol{e}_{z}^{(1)}&\boldsymbol{e}_{z}^{(2)}\cdot\boldsymbol{e}_{z}^{(1)}\end{pmatrix}}_{=:\boldsymbol{A}}\underbrace{\begin{pmatrix}x_{2}\\\
y_{2}\\\ z_{2}\end{pmatrix}}_{\boldsymbol{r}_{\Sigma_{2}}}$ (24)
$\displaystyle+$ $\displaystyle\underbrace{\begin{pmatrix}t_{x}\\\ t_{y}\\\
t_{z}\end{pmatrix}}_{\boldsymbol{t}}\,.$

Thus, the transformations between the coordinate systems
$\Sigma_{1}\leftrightarrow\Sigma_{2}$ are readily given by:

$\displaystyle\boldsymbol{r}_{\Sigma_{1}}$ $\displaystyle=$
$\displaystyle\left[\boldsymbol{S}\boldsymbol{A}\right]\cdot\boldsymbol{r}_{\Sigma_{2}}+\boldsymbol{t}$
(25) $\displaystyle\boldsymbol{r}_{\Sigma_{2}}$ $\displaystyle=$
$\displaystyle\left[\boldsymbol{S}\boldsymbol{A}\right]^{-1}\cdot\left[\boldsymbol{r}_{\Sigma_{1}}-\boldsymbol{t}\right]\,.$
(26)

##  Appendix B Test calculations

Figure 5: Top panel: Synthetic line profiles as a function of frequency shift
from line centre in units of the terminal velocity for our spherically
symmetric test model (see text). The solid lines indicate the line profiles as
obtained from the binary code at different phases $\varphi$ and a fixed
inclination $i=90^{\circ}$, and the grey crosses correspond to the solution as
obtained from the single-star code. The vertical dashed lines indicate the
transition frequency of the line at the different phases for the assumed
orbital configuration. Bottom panel: Relative error of the obtained line
profiles at different phases when compared to the single-star code (shifted by
the corresponding projected orbital velocities).

In this section, we test the developed BOSS-3D code by calculating synthetic
line profiles of a generic resonance-line transition within a spherically
symmetric wind of a prototypical hot, massive star. In order to compare the
results with the single-object code described in Sect. 2, the secondary object
is assumed to be dark ($T_{\rm eff}^{(2)}=10\,{\rm K}$) and small
($R_{\ast}^{(2)}=0.2\,R_{\odot}$). The coordinate system of the secondary
object, however, has been set to a large extent, in order to test if the
routines for setting up the triangulation and the individual rays are
performing reasonably well. By assigning an artificial circular orbit with
absolute velocity $v_{\rm orb}=1000\,{\rm km}\,{\rm s}^{-1}$ to the primary
object, we further test if the corresponding Doppler-shifts are correctly
implemented for different viewing angles. The primary object is described by
$R_{\ast}^{(1)}=20\,R_{\odot}$, $T_{\rm eff}^{(1)}=40\,{\rm kK}$, and with a
wind stratification given from a $\beta$-velocity law with base velocity
$v_{\rm min}=10\,{\rm km}\,{\rm s}^{-1}$, terminal velocity
$v_{\infty}=2000\,{\rm km}\,{\rm s}^{-1}$, and $\beta=1$. The required density
is obtained from the continuity equation assuming a mass-loss rate
$\dot{M}=10^{-6}M_{\odot}{\rm yr}^{-1}$, and the source function is calculated
from the Sobolev approximation (Sobolev 1960) for a pure scattering line
within the two-level approach. Finally, we specify the boundary condition by
the Planck function, $I_{\nu}=B_{\nu}(T_{\rm eff}^{(1,2)})$ , if a ray hits
the core of one of the objects.

Figure 5 displays the resulting line profiles for edge-on observer directions
at different phases, compared with the result from the single-object code.
Since the relative error remains below few percent, we are highly confident
that the algorithm also works for more complex simulations. Moreover, the
expected velocity shift at the different phases is exactly reproduced. We
finally emphasize that we obtain the same results when interchanging the
primary and secondary object within our algorithm, as required.

##  Appendix C Parameter study

In this section, we perform a parameter study of the disc model introduced in
Sect. 4, and investigate the main effects of the different model parameters on
the various measured quantities (i.e. the line-profile shape, the equivalent
width, the radial-velocity curve and its semi-amplitude). To this end, we
focus on the B-star and BH-disc scenario (model PB_SBH), and on the stripped-
star and Be-star scenario with two discs (model PST_SBE02). The individual
disc parameters are varied one at a time, starting from our best-fit by eye
models presented in Sect. 4 and summarized in Table 3. The radial-velocity
curves have been derived from the Hα line wings defined as those parts of the
synthetic line profile with normalized flux below a threshold of $1/3$ of the
flux peak (see also Liu et al. 2019).

###  C.1 Parameter study for the B-star and BH-disc scenario

Figure 6: Parameter study with base parameters defined in Table 3 for model
PB_SBH. Upper left main-panel: Synthetic line profiles at phase $\varphi=0$
(top left sub-panel), radial-velocity curve (top right sub-panel), equivalent
width with corresponding reduced $\chi^{2}$-value (middle left sub-panels),
semi-amplitude of the radial velocity curve with corresponding reduced
$\chi^{2}$-value (middle right sub-panels), and total reduced $\chi^{2}$-value
(bottom panel), all as a function of the disc’s base density,
$\rho_{0}^{\rm(2)}$. The black dotted line indicates our best-fit by eye (see
Sect. 4.1), and the blue dashed lines display the observed values with
corresponding error-bars (blue-shaded area). Similarly, the top right to
bottom right main-panels indicate the dependence of our measurable quantities
on the disc’s temperature, $T_{\rm disc}^{(2)}$, the slope of the density
stratification, $\beta_{\rm D}^{(2)}$, and the micro-turbulent velocity,
$v_{\rm micro}^{\rm(2)}$.

Figure 6 displays the results of our parameter study for the B-star and BH-
disc scenario, with variations of the four free parameters, namely the disc’s
base density, $\rho_{0}^{\rm(2)}\in 3\cdot[10^{-13},10^{-11}]\,{\rm g}\,{\rm
cm}^{-3}$, the slope of the density stratification, $\beta_{\rm
D}^{(2)}\in[1.5,4]$, the disc’s temperature $T_{\rm
disc}^{\rm(2)}\in[6,15]\,{\rm kK}$, and the micro-turbulent velocity, $v_{\rm
micro}^{\rm(2)}\in[15,105]\,{\rm km}\,{\rm s}^{-1}$. Additionally, we display
the reduced $\chi^{2}$-values of the equivalent width and the radial-velocity
semi-amplitude as compared to the observations, as well as the combined
reduced $\chi^{2}$-value (with an equal weight to both quantities).

#### Impact of $\rho_{0}^{\rm(2)}.$

The strength of the line (and thus also the equivalent width) increases with
increasing base density of the disc. Indeed, the equivalent width follows
roughly a quadratic dependence, following the $\rho^{2}$-dependence of the
opacity for recombination lines. For low densities, only the photospheric line
profile from the primary object (the B star in this scenario) is visible.
Thus, the radial velocity curve when measured from the Hα-line wings follows
the B-star orbit ($K_{1}=53\,{\rm km}\,{\rm s}^{-1}$) for low base densities,
and the BH orbit ($K_{2}=7.4\,{\rm km}\,{\rm s}^{-1}$) for high base
densities, as expected. Additionally, the overestimation of the radial-
velocity amplitude due to the anti-phase motion of the stellar absorption
profile (see Abdul-Masih et al. 2020a) is clearly observed within our models
at intermediate densities.

#### Impact of $T_{\rm disc}^{\rm(2)}.$

The strength of the Hα line decreases with increasing disc temperature. This
effect can be explained as follows. With increasing temperature, the
occupation numbers of the lower level become diminished due to the Boltzmann
factor. Thus, the opacity is considerably reduced, and the disc essentially
becomes transparent in Hα at high temperatures. We emphasize that also at
temperatures $T_{\rm disc}^{\rm(2)}\lesssim 6\,{\rm kK}$ (not shown here), the
opacity drops significantly since most hydrogen will be found in the ground
state. As before, the radial velocity curve measured from the Hα line follows
the BH orbit for low temperatures (and thus a strong disc emission). With
increasing disc temperature, the semi-amplitude of the radial velocity curve
becomes overestimated when measured from the Hα-line wings due to the anti-
phase motion of the stellar absorption profile. In contrast to the low-density
regime described above, the Hα-line wings are mainly controlled by the disc
emission at the upper temperature limit, and the radial-velocity curve is
therefore following the BH orbit for the full parameter space considered here.

#### Impact of $\beta_{\rm D}^{\rm(2)}.$

With increasing slope of the disc’s density stratification, the Hα-line
profile switches rapidly from showing a significant disc emission to the
(almost) pure absorption profile from the B star. Clearly, for $\beta_{\rm
D}^{\rm(2)}\gtrsim 2$, the disc emission vanishes due to decreased densities
(in particular at large distances from the BH and therefore low velocities).
The corresponding radial velocity curves are slightly underestimating the true
orbital velocity of the B star ($K_{1}=53\,{\rm km}\,{\rm s}^{-1}$) due to a
slight emission from the disc in the Hα-line wings (probing the high-velocity
regions and thus the inner parts of the disc).

#### Impact of $v_{\rm micro}^{\rm(2)}$.

Both the strength of the line and the radial velocity curve are only mildly
affected by the micro-turbulent velocity in the disc. The shape of the Hα-line
profile, however, switches from a relatively narrow and strongly pronounced
double-peaked feature at low $v_{\rm micro}^{\rm(2)}$ to a broad and smooth
one at high $v_{\rm micro}^{\rm(2)}$. When compared to observations (see Sect.
4), the latter is to be preferred.

#### Investigation of $\chi^{2}$.

The $\chi^{2}$-values shown in Fig. 6 indicate that both the slope of the
density stratification and the micro-turbulent velocity are relatively well
constrained to the best-fit by eye values found in Sect. 4.1. There exists,
however, a slight degeneracy in the disc’s base density and temperature. For
instance, one might find a similar good model by decreasing both the disc
density and the disc temperature at the same time. With such a model, however,
the observed anti-phase behaviour of line core and line wings (see Fig. 2)
could not be reproduced either. Moreover, at temperatures $T_{\rm
disc}^{(2)}\lesssim\,6{\rm kK}$ (not shown here) and $T_{\rm
disc}^{(2)}\gtrsim\,15{\rm kK}$, the disc’s base density would need to be
increased to obtain a significant equivalent width, translating then to an
X-ray bright disc that has not been observed.

###  C.2 Parameter study for the stripped-B and Be scenario with two discs

Figure 7: As Fig. 6, however for model PST_SBE02 and focussing on the six free
parameters $\rho_{0}^{\rm(1)}$, $\rho_{0}^{\rm(2)}$, $\beta_{\rm D}^{(1)}$,
$\beta_{\rm D}^{(2)}$, $v_{\rm micro}^{(1)}$, and $v_{\rm micro}^{(2)}$.

Figure 7 displays the results of our parameter study for the stripped-B and Be
scenario with two discs, where the discs are described by six free parameters
in total, namely the base density for each disc,
$\rho_{0}^{\rm(1)}\in[10^{-14},10^{-11}]\,{\rm g}\,{\rm cm}^{-3}$ and
$\rho_{0}^{\rm(2)}\in 6\cdot[10^{-13},10^{-11}]\,{\rm g}\,{\rm cm}^{-3}$, the
slope of the density stratification, $\beta_{\rm D}^{(1)}$ and $\beta_{\rm
D}^{(2)}$ (both in a range $\in[1.5,4]$), and the micro-turbulent velocities,
$v_{\rm micro}^{(1)}$, $v_{\rm micro}^{(2)}$ (both in a range
$\in[15,105]\,{\rm km}\,{\rm s}^{-1}$). As before, we additionally display the
reduced $\chi^{2}$-values calculated for the equivalent width, the radial
velocity semi-amplitude, and the equally weighted sum of both.

#### Impact of $\rho_{0}^{\rm(1)}$ and $\rho_{0}^{\rm(2)}$.

The strength of the line increases with increasing base density of the
stripped-star disc, $\rho_{0}^{\rm(1)}$, particularly in the low-velocity
regions (due to the overall low velocities in the stripped-star disc). At the
lower end of our parameter space, the stripped-star disc becomes transparent
and the Hα-line profile consists only of three components, namely the
absorption profiles of both stars in the binary system, and the disc emission
from the secondary object (the Be-star disc). Consequently, the radial-
velocity curves vary significantly. For low $\rho_{0}^{\rm(1)}$, the Hα-line
wings are probing the Be-star orbit (including the effects of all underlying
stellar absorption profiles, of course). With increasing $\rho_{0}^{\rm(1)}$,
the radial-velocity amplitude becomes reduced and switches sign at the upper
limit of our parameter space due to a significant contamination of the Hα-line
emission from the stripped-star disc.

In contrast to the effects of $\rho_{0}^{\rm(1)}$, the strength of the Hα-line
profile increases in both the low and high velocity regions with increasing
$\rho_{0}^{\rm(2)}$. At the lower end of our parameter space, the line
profiles consist only of the emission from the stripped-star disc (and again
the stellar absorption components). Thus, the radial velocity curves as
obtained from the Hα-line wings follow the stripped-star orbit at low
$\rho_{0}^{\rm(2)}$, while the anti-phase behaviour probing the Be-star orbit
is measured at high $\rho_{0}^{\rm(2)}$.

#### Impact of $\beta_{\rm D}^{\rm(1)}$ and $\beta_{\rm D}^{\rm(2)}$.

The strength of the Hα-line profile increases with decreasing slope of the
stripped-star disc’s density stratification, $\beta_{\rm D}^{\rm(1)}$, in
particular in the low-velocity regions. As in Sect. C.1, increasing
$\beta_{\rm D}^{\rm(1)}$ translates to decreasing the density of the stripped-
star disc, particularly in the outer (low-velocity) regions. Thus, the
synthetic line profiles become slightly broader at the upper limit of
$\beta_{\rm D}^{\rm(1)}$. Since the stripped-star disc essentially becomes
transparent at high $\beta_{\rm D}^{\rm(1)}$, the obtained radial-velocity
curve follows the Be-star orbit. Accordingly, the stripped-star orbit is being
probed at low $\beta_{\rm D}^{\rm(1)}$, with a significant contamination from
the Be-star disc though.

Also with decreasing slope of the Be-star disc’s density stratification, the
strength of the Hα-line profile increases, now within both the low and high
velocity regions though. At the upper limit of $\beta_{\rm D}^{\rm(2)}$, the
Be-star disc becomes transparent and only the stripped-star disc is visible
(together with the stellar absorption profiles, again). Thus, a high
$\beta_{\rm D}^{\rm(2)}$ is favourable for measuring the stripped-star orbit,
while a low $\beta_{\rm D}^{\rm(2)}$ would give the Be-star orbit (still
contaminated by the stellar absorption profiles).

#### Impact of $v_{\rm micro}^{\rm(1)}$ and $v_{\rm micro}^{\rm(2)}$.

While the equivalent width of the Hα-line profile is only slightly increasing
with increasing micro-turbulent velocity of the stripped-star disc ($v_{\rm
micro}^{\rm(1)}$), the shape of the line changes significantly from a double-
peaked feature to a relatively smooth and slightly broader appearance. Thus,
the radial-velocity curve is probing primarily the Be-star disc at the lower
limit of $v_{\rm micro}^{\rm(1)}$, and the semi-amplitude is increasing with
$v_{\rm micro}^{\rm(1)}$ since the line-profile wings become increasingly
contaminated from the stripped-star disc emission.

With increasing micro-turbulent velocity of the Be-star disc ($v_{\rm
micro}^{\rm(2)}$) now, the equivalent width of the Hα-line profile is again
only slightly increasing while the shape of the line changes significantly
from a relatively narrow double-peaked line profile to a much smoother and
broader one. Thus, the radial-velocity curves are probing primarily the Be-
star disc at the upper limit of $v_{\rm micro}^{\rm(2)}$. With decreasing
micro-turbulent velocity, the line wings are contaminated by the stripped-star
disc, and the semi-amplitude of the radial-velocity curve therefore becomes
decreased (towards positive values).

#### Investigation of $\chi^{2}$.

As shown by the reduced $\chi^{2}$-values of our stripped-star and Be-star
models (Fig 7), all fit parameters are fairly well constrained to the best-fit
by eye values found in Sect. 4.2. There might be a slight degeneracy of the
base density and the slope of the density stratification within each of the
discs, though. For instance, reducing $\rho_{0}^{\rm(1)}$ and $\beta_{\rm
D}^{(1)}$ (or $\rho_{0}^{\rm(2)}$ and $\beta_{\rm D}^{(2)}$) at the same time
might give similar results as found for our best-fit model.