Issue 
A&A
Volume 660, April 2022



Article Number  A17  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202141831  
Published online  01 April 2022 
Binaryobject spectralsynthesis in 3D (BOSS3D)
Modelling H_{α} emission in the enigmatic multiple system LB1
^{1}
Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium
email: levin.hennicker@kuleuven.be
^{2}
National Solar Observatory, 22 Ohi’a Ku Street, Makawao, HI 96768, USA
^{3}
Anton Pannekoek Institute for Astronomy, University of Amsterdam, Postbus 94249, 1090 GE Amsterdam, The Netherlands
^{4}
European Southern Observatory, KarlSchwarzschildStrasse 2, 85748 Garching bei München, Germany
^{5}
European Southern Observatory, Alonso de Cordova 3107, Vitacura, Casilla, 19001 Santiago de Chile, Chile
^{6}
Institut de Planétologie et d’Astrophysique de Grenoble, 414 Rue de la Piscine, 38400 SaintMartind’Hères, France
Received:
20
July
2021
Accepted:
23
November
2021
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 LB1 constitutes a perfect example of such a complex multiD problem. Thus far, the LB1 system has been indirectly investigated by 1D stellaratmosphere codes and by spectral disentangling techniques, yielding differing conclusions about the nature of the system (e.g., a Bstar and blackhole binary with an accretion disc around the black hole or a strippedstar and Bestar binary system have been proposed).
Aims. To improve our understanding of the LB1 system, we directly modelled the phasedependent H_{α} line profiles of this system. To this end, we developed and present a multipurpose binaryobject spectralsynthesis code in 3D (BOSS3D).
Methods. BOSS3D calculates synthetic line profiles for a given state of the circumstellar material. The standard pzgeometry 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 LB1 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 phasedependent 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 Bestar disc or the blackhole 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 postasymptotic giant branch binaries including circumstellar and circumbinary discs and massivestar binaries with stellar winds and disc systems. For the LB1 system, our modelling provides strong evidence that each object in the system contains a disclike structure.
Key words: radiative transfer / methods: numerical / stars: emissionline, Be / stars: black holes / binaries: spectroscopic
© ESO 2022
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_{*}, T_{eff}, R_{*}), 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 OBstar regime (on which this paper focusses), stateoftheart atmospheric modelling and spectral synthesis codes typically assume spherical symmetry (e.g., PHOENIX: Hauschildt 1992; CMFGEN: Hillier & Miller 1998, Hillier 2012; WMBASIC: Pauldrach et al. 2001; PoWR (see also Sander et al. 2019, Sect. 1): 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 nonLTE occupation numbers within multiD calculations, corresponding spectral synthesis codes are only gradually developed by applying different solution methods from finitevolume methods (FVM) to short and longcharacteristics methods (SC and LC) to MonteCarlo methods (MC). A nonexhaustive list of 3D radiative transfer codes accounting for supersonic velocity fields includes Wind3D (see also Lobel & Blomme 2008, Sect. 1; 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 twolevelatom (TLA) approach, whereas the last two include a multilevel 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 multiD alternative to calculating synthetic spectra from massivestar binary systems with circumstellar material around them. For example, the SPAMMScode by AbdulMasih 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 (BOSS3D, i.e. binaryobject spectralsynthesis in 3D)^{1} 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 largescale, 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 singlestar pzgeometry 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 LB1. The spectrum of this system clearly shows an antiphase 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 ≈ 79 d, consisting of a B star (visible component) and a 70 M_{⊙} black hole (BH) with associated disc yielding the H_{α} emission. Since it is difficult to explain the origin of a 70 M_{⊙} BH in the Milky Way from stellar evolution theory (though, see Belczynski et al. 2020), there has been vivid discussion about LB1 in the literature. While SimónDíaz et al. (2020) reanalysed the visible component (B star, primary object) and found a significantly lower mass for the B star thus also reducing the BH mass, AbdulMasih et al. (2020a) and ElBadry & 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 LB1 and proposed a binary system consisting of a stripped Hestar (the primary) in a 79 d orbit with a rapidlyrotating Bestar, presumably the product of a recent masstransfer event. Alternatively, Rivinius et al. (2020) proposed that the Bestar is a tertiary object, leaving the secondary component unidentified. Recently, Lennon et al. (2021) modelled the spectral energy distribution (SED) of LB1 in the optical and ultraviolet regime using 1D planeparallel modelatmospheres, and compared several models (including the aforementioned strippedstar/Bestar and Bstar/BH hypotheses) with observations from the Hubble Space Telescope. While their models slightly favour the Bstar and BH scenario, the strippedstar and Bestar 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 LB1 by applying the BOSS3D 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 phasedependent H_{α} line profile provided the LB1 system contains an additional disc attached to the Bstar 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 LB1 system. Finally, we summarize our findings in Sect. 5.
2. Singlestar spectral synthesis
In this section, we briefly summarize the basic solution method for obtaining synthetic spectra from a given singlestar model in 3D. Following, for instance, Busche & Hillier (2005), the radiation flux at frequency ν emerging from any unresolved emitting object as seen by an observer at distance d can be formulated as
where R_{max} is the maximum size of the emitting object, (p,ζ) are polar coordinates describing the projected disc of the emitting object perpendicular to the observer’s direction n, and I_{ν} is the emergent specific intensity into that direction evaluated at a distance .
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}, ζ_{j}) domain of a cylindrical coordinate system with corresponding Cartesian reference frame , and with the axis being aligned with the direction to the observer (see also Fig. 1). The corresponding (timeindependent) equation of radiative transfer reads:
Fig. 1. Geometry used for the singlestar algorithm (top panels) and for the binary code (bottom panels). Top left panel: cylindrical coordinate system (blue) with corresponding Cartesian reference frame (red) for an observer’s direction used to obtain the synthetic spectrum of a given (single) object described in a spherical coordinate system (magenta) with corresponding Cartesian reference frame (black). Top middle panel: p − ζ (or plane) of the cylindrical coordinate system (perpendicular to the observer’s direction). Top right panel: an arbitrary slice, with 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 coordinates indicated by the red dots. Bottom left panel: two local coordinate systems Σ_{1} and Σ_{2} describing two individual objects indicated by the cyan and magenta quartercircles, respectively. Both coordinate systems are embedded within a global coordinate system Σ_{0} (typically but not necessarily chosen to be the centre of mass of both objects). The vectors t_{1} and 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 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 rays propagating through the system of object 1 (blue) and object 2 (red) potentially need to be merged. 
with the opacity χ_{ν} and source function S_{ν} 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 ↔ u with transition frequency ν_{lu}. Further, we omit the notation for explicit spatial dependencies. Then, the opacity and source functions are given by:
where n_{l} and n_{u} are the occupation numbers, g_{l} and g_{u} are the statistical weights, and
is the profile function (approximated here by a Doppler profile) at observers frame frequency ν including nonrelativistic Doppler shifts from the observer’s to the comoving frame via the projected velocity n ⋅ v. The width of the profile, , is given from the thermal velocity of the considered species with mass number m_{A}, accounting for the local temperature T and a microturbulent velocity v_{micro}:
For a given atmospheric model (with known temperatures, velocity vectors, occupation numbers, and microturbulent velocities), the equation of radiative transfer, Eq. (2), becomes an ordinary differential equation with exact solution from any point to
We note that is the opticaldepth 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 opticaldepth 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 with corresponding Cartesian reference frame , and if we assume that all structures are well resolved by the corresponding discretized grid, the discretization of is easily constructed following the standard pzgeometry approach (e.g., Mihalas 1978, see also Fig. 1):
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}, ζ_{j}) pair are determined in the cylindrical coordinate system . 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 :
with e_{x, y, z} and the unit vectors of the Cartesian reference frames for the spherical and cylindrical systems, respectively, and A the transformation matrix from the Cartesian ray system to the object’s Cartesian reference frame. The unit vectors and can be chosen arbitrarily under the constraint that and . With the obtained coordinates in the spherical system of the object, all quantities on the discretized ray are calculated by trilinear interpolations. To ensure that the lineprofile function along a ray is resolved by the grid, we refine the grid if the projected velocity steps are above v_{th}/3.
For a single object, synthetic line profiles are thus obtained for a specific observer’s direction by performing the following steps:

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

Calculation of the grid for each (p_{i}, ζ_{j}) following Eq. (8).

Transformation of the coordinates to the object’s spherical system by Eqs. (9)–(11), and interpolation of all required quantities (i.e. opacities, source functions, and velocity components).

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 userspecified subroutine to couple either line profiles as obtained from Kurucz atmosphere models, e.g., Castelli & Kurucz 2003, or from FASTWIND)^{2}. For noncore rays, the incident intensity at the outer boundary is set to zero.

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

Restarting the procedure at step 2 for the next frequency bin.
3. Multiobject 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 meshrefinement 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 and the individual systems (see Appendix A):
where and 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, t_{q} is the translation vector to the origin of the individual systems measured in the global system, and the transformation matrices Q_{q} are given by:
The unit vectors of the individual coordinate systems are described in the global system , 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 n_{0} measured in the global frame then, we compute the corresponding representations in the individual coordinate systems, , and assign cylindrical coordinate systems to each of the objects following the procedure outlined in Sect. 2 by setting . Again, to keep track of the different coordinate systems, we define the transformation from the local raycoordinates to the local object coordinates:
with transformation matrix
and e_{x} = (1, 0, 0), e_{y} = (0, 1, 0), 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 , 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 Delaunaytriangulation of all (p_{i}, ζ_{j})_{q} points from all individual cylindrical coordinate systems by using the GEOMPACK2 library^{3} (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:

Definition of the individual coordinate systems for all objects q, that is lengthscales L_{q}, translation vectors t_{q}, basis vectors (described in the global coordinate system) , global velocities v^{(q)} (corresponding to the orbital velocities of the objects), and calculation of the transformation matrices Q_{q} via Eq. (14).

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

Transformation of all (p, ζ)_{q} coordinates to the global raycoordinate system using Eqs. (12) and (15) to obtain the corresponding coordinates in the global ray coordinate system.

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

Backtransformation of each coordinate to all individual coordinate systems via Eqs. (13) and (15).

Step 2 of the singlestar algorithm: Setting the grid along a ray within each individual coordinate system following Eq. (8).

Step 3 of the singlestar algorithm: Transformation of all coordinates to the individual spherical systems and interpolation of all required physical quantities (i.e. opacities, source functions, velocity fields) onto the ray.

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

Merging all individual rays to one single ray for a given point (see Fig. 1).

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

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

Step 6 of the singlestar 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 postasymptotic 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. LB1
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 singlestar models). In this section we model the H_{α} line of the binary (or multiple) system LB1 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ónDí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ónDí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 sphericallysymmetric NLTE code PoWR (Hamann & Gräfener 2003; Sander et al. 2017) to the disentangled spectrum of the Bestar in the strB+Be scenario, while relying on the Grid Search in Stellar Parameters tool (GSSP, see also Shenar et al. 2020, Appendix B.1; Tkachenko 2015) with synthetic spectra derived from 1D planeparallel LTE modelling using the SYNTHV (see also Shenar et al. 2020, Appendix B.1) radiative transfer code (Tsymbal 1996) and a grid of LLMODEL (see also Shenar et al. 2020, Appendix B.1) atmospheres (Shulyak et al. 2004) when analysing the stripped star. The updated stellar parameters by Lennon et al. (2021) using the 1D planeparallel 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 Bestar disc to calculate synthetic H_{α} line profiles.
Stellar parameters for the B+BH and strB+Be scenario used for our simulations of the LB1 system, and as found in the literature.
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 powerlaw, and assumes hydrostatic equilibrium in the vertical direction. Then, the density in the disc is:
with r, Θ the radial and colatitudinal coordinates, M_{*} the stellar or BH mass, R_{*} the stellar radius or the minimum radius of the BH disc, and ρ_{0} and β_{D} the base density and powerlaw index for the density stratification. When β_{D} = 15/8, this formulation is equivalent to the standard αaccretion disc model by Shakura & Sunyaev (1973) at large distances from the BH, with the αprescription and accretion rate hidden^{4} within the parameter ρ_{o}. As such, ρ_{0} is related to the mass flux through the disc, and β_{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 Bestar decretion disc). For Bestars, typical values are of the order of β_{D} ∈ [1.5, 4] and ρ_{0} ∈ 5 ⋅ [10^{−10}, 10^{−13}] g 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, , evaluated here for a typical temperature T = 10 kK and a mean molecular weight μ = 0.6, the disc density is specified for given stellar parameters and the input parameters ρ_{0}, β_{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
Since the temperature stratification crucially depends on the type of the disc (BHaccretion disc or Bestar 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 SahaBoltzmann 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 semiamplitude for the primary and secondary of K_{1} ≈ 53 km s^{−1} and K_{2} ≈ [6, 13] km 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 semiamplitude 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 . However, as shown by AbdulMasih 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 for LB1, in good agreement with a detailed investigation of nearinfrared emission lines by Liu et al. (2020) estimating . 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 semiamplitude when applying the barycentric method to the synthetic H_{α} lines, while also yielding from the actual orbit of the system.
Fig. 2. H_{α} line profiles at different phases φ (top row) with corresponding meansubtracted dynamical spectra (middle row) and radial velocity curves (bottom row) of the BH or strB companion. The x_{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: Bestar, model 1), PST_SBE02 (primary: stripped star, secondary: Bestar, 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, meansubtracted dynamical spectrum and radial velocity curves of LB1, 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 km 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 (phaseaveraged) equivalent width evaluated in velocity space in units of 100 km s^{−1} and the inclination used within the syntheticspectra calculations are indicated. 
Since the eccentricity is small (e = 0.03 ± 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 ≈ 79 d then, the orbit is obtained from the twobody problem:
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 semiamplitude of the radial velocity curve for the primary object by v_{1} ⋅ sin(i) = K_{1}. For our best models, the orbital parameters are summarized in Table 2.
Orbital parameters for the LB1 system as calculated from our models and found in the literature.
4.1. Bstar and BHdisc (B+BH) scenario
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 given, the mass of the black hole is constrained via Eqs. (20) and (21) yielding , where we have chosen the blackhole mass as M_{2} = 30 M_{⊙} and inclination i = 22° in order to reproduce the semiamplitude of the radial velocity K_{1}. With this inclination, we adopted a rotational velocity for the B star, (with the observed value v sin i = 8 km s^{−1}, Lennon et al. 2021). To model the disc, we have set the minimum radius to , 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 α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 Bstar companion completely dominates the total radiation flux. In the B+BH binary system, however, the disc will also be heated by the Bstar companion. In order to avoid complex multiD radiationhydrodynamic 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_{disc}, the microturbulent velocity in the disc, , the base density of the disc, , and the slope of the density stratification, .
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_{disc} = 8 kK, microturbulent velocity , base density , and radial density stratification β_{D} = 3/2 are shown in Fig. 2, with corresponding density, temperature, and velocity profiles in the equatorial plane shown in Fig. 3.
Fig. 3. Disc models in the equatorial plane corresponding to our bestfit parameters of the LB1 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 Bestar disc within the singledisc model (model PST_SBE01), respectively. For model PST_SBE02 (containing two discs), the Bestar disc is indicated by the green dotted line (partially overlapping with the red solid line from model PST_SBE01), and the strippedstar disc is indicated by the green solid line. 
The slope of the density stratification is in fairly good agreement with accretion discs described within the α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 lowvelocity 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 microturbulent 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 doublepeaked emission feature of accretion discs even at such inclinations (i = 22°). High turbulent velocities, however, are frequently found when modelling dynamical systems with simplified stationary models (which by assumption neglect smallscale 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 km s^{−1} Gaussian ‘macroturbulence’ 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 phaseaveraged equivalent width measured in velocity space in units of 100 km s^{−1}, , in good agreement with observations (). We emphasize, however, that the radial velocity curves show some discrepancies when compared to observations (see also Table 2). While the actual orbital velocities, , are slightly below the observed range (), we find a higher semiamplitude when applying the barycentric method (). Indeed, this discrepancy can be explained by the absorption profile of the B star (see discussion in AbdulMasih 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 (meansubtracted) dynamical spectra of the H_{α} line. While the observed dynamical spectrum clearly displays an antiphase 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_{disc} = 8 kK in order to reproduce the observations. This temperature region can indeed be considered as a sweet spot, since lower temperatures (T_{disc} ≲ 5 kK) would require much higher densities to provide enough emission in H_{α}. The corresponding densities would translate to a high massaccretion rate of the BH within the standard αdisc prescription, thus predicting an Xray bright accretion disc that has not been observed. Similarly, too high temperatures (T_{disc} ≳ 15 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 Xray bright accretion disc. If the sweetspot temperature as found in our model is a reasonable substitute for detailed numerical simulations still needs to be investigated.
4.2. Strippedstar and Bestar (strB+Be) scenario
For the stellar parameters of the LB1 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°. Motivated by these authors, we have set the rotational velocity of the stripped star to (with v sin i = 7 km s^{−1} following Shenar et al. 2020 and Lennon et al. 2021), and assumed the Bestar to rotate at its critical velocity, . 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 radiativetransfer calculations in Carciofi & Bjorkman (2006, Sect. 5.2). In regions close to the stellar surface (where the vertical electronscattering optical depth τ_{z}(r sin Θ) > 0.1), we adopted an analytical temperature stratification described by a flat blackbody reprocessing disc (their Eq. (12)):
with T_{0} the base temperature of the disc. Since the base densities of our disc models are typically relatively low, we can neglect backwarming effects of the photosphere by the disc, and therefore have set the base temperature to the stellar effective temperature, T_{0} = T_{eff}. To mimick the heating of the disc in (vertically) optically thin regions by the star, we adopted a constant disc temperature T_{disc} = 0.6 ⋅ T_{eff} where τ_{z}(r sin Θ) < 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 semiamplitude of the radial velocity curve of the primary, K_{1}, the model is completely specified by six input parameters, namely the microturbulent velocities, , the base density , and the slope of the density stratification of both discs, . When considering only the Bestar 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 BOSS3D 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.
Disc parameters for the three different (bestfit) models as obtained from our simulations of the LB1 system (see text).
Model PST_SBE01. This model describes a stripped star (primary object) in orbit with a Bestar (secondary object) hosting a circumstellar disc. The obtained disc parameters, and 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 continuum^{5} 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 microturbulent velocity to obtain the broad H_{α} emission wings and to smooth out the otherwise doublepeaked 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 (, again measured in 100 km s^{−1}) is somewhat smaller than expected from observations (). Moreover, the obtained semiamplitude of the radial velocity of the secondary object, , shows a clear deviation from the actually observed values (). As in the B+BH scenario, the semiamplitude is larger than the orbital velocity of the secondary object (v_{2}sin(i) = 9.9 km 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 lowvelocity emission, we additionally implemented a circumstellar disc attached to the stripped star. Although speculative, this disc might have formed by reaccretion from the Bestar (e.g., Shenar et al. 2020) or from earlier masstransfer phases.
While the best fit (by eye) gives a similar disc for the Bestar as found above, we predict a slightly higher base density of the strippedstar disc, , and a relatively steep slope of the density stratification, . Since the orbital velocities of the strippedstar 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 microturbulent velocities. The qualitative behaviour of the observed H_{α} line profile then is very well reproduced, with corresponding equivalent width, , in very good agreement with observations. Since the H_{α} line wings are slightly contaminated by the (antiphase) emission from the strippedstar disc, the obtained radial velocity curve () agrees fairly well with the corresponding observations ( 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 AbdulMasih 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 (βvelocitytype) 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 LB1 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 superlinearly) 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 lengthscales of each object. The presented algorithm and associated code (BOSS3D) 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 postasymptotic 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 Bstar/BHdisc (B+BH) and strippedstar/Bestar (strB+Be) scenarios for the enigmatic (multiple) system LB1 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 Bestar disc. The level populations for the H_{α} line transition have then been computed from SahaBoltzmann statistics assuming full ionization.
Under these assumptions none of the models can reproduce the detailed phasedependent shape of the observed line profiles, particularly due to differences in the dynamical spectrum (and missing lowvelocity 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 AbdulMasih et al. 2020a). It remains open whether NLTE calculations can help to explain the detailed phasedependent 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 masstransfer phases or could have originated from reaccretion of material from the Bestar 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 LB1 system. In any case, our findings provide strong evidence that LB1 contains a discdisc 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).
Following 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 , with Ṁ_{16} the massaccretion rate in 10^{16} g s^{−1}, m_{1} the mass of the central object in M_{⊙}, and R_{*, 10} the inner disc radius in units of 10^{10} cm.
Acknowledgments
We thank our referee, Dr. Maria Bergemann, for many helpful comments and suggestions. L. H. and J. O. S. gratefully acknowledge support from the Odysseus program of the Belgian Research Foundation Flanders (FWO) under grant G0H9218N. J. B. 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. T.S. 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łodowskaCurie grant agreement No 101024605. I.E.M. 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
 AbdulMasih, M., Sana, H., Conroy, K. E., et al. 2020a, A&A, 636, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AbdulMasih, M., Banyard, G., Bodensteiner, J., et al. 2020b, Nature, 580, E11 [Google Scholar]
 Belczynski, K., Hirschi, R., Kaiser, E. A., et al. 2020, ApJ, 890, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Bollen, D., Kamath, D., De Marco, O., Van Winckel, H., & Wardle, M. 2020, A&A, 641, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Busche, J. R., & Hillier, D. J. 2005, AJ, 129, 454 [Google Scholar]
 Carciofi, A. C., & Bjorkman, J. E. 2006, ApJ, 639, 1081 [Google Scholar]
 Castelli, F., & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, eds. N. Piskunov, W. W. Weiss, & D. F. Gray, 210, A20 [Google Scholar]
 ElBadry, K., & Quataert, E. 2020, MNRAS, 493, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition [Google Scholar]
 Hamann, W. R., & Gräfener, G. 2003, A&A, 410, 993 [CrossRef] [EDP Sciences] [Google Scholar]
 Hauschildt, P. H. 1992, J. Quant. Spectr. Rad. Transf., 47, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Hauschildt, P. H., & Baron, E. 2006, A&A, 451, 273 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2018, A&A, 616, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2020, A&A, 633, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hillier, D. J. 2012, in From Interacting Binaries to Exoplanets: Essential Modeling Tools, eds. M. T. Richards, & I. Hubeny, 282, 229 [NASA ADS] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I. 1988, Comput. Phys. Commun., 52, 103 [Google Scholar]
 Hubeny, I., & Lanz, T. 1995, ApJ, 439, 875 [Google Scholar]
 Joe, B. 1991, Adv. Eng. Software Workstations, 13, 325 [CrossRef] [Google Scholar]
 Kee, N. D., Owocki, S., & Kuiper, R. 2018, MNRAS, 474, 847 [NASA ADS] [CrossRef] [Google Scholar]
 Kunasz, P., & Auer, L. H. 1988, J. Quant. Spectr. Rad. Transf., 39, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Lamers, H. J. G. L. M., CerrutiSola, M., & Perinotto, M. 1987, ApJ, 314, 726 [NASA ADS] [CrossRef] [Google Scholar]
 Lennon, D. J., Maíz Apellániz, J., Irrgang, A., et al. 2021, A&A, 649, A167 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Liu, J., Zhang, H., Howard, A. W., et al. 2019, Nature, 575, 618 [Google Scholar]
 Liu, J., Zheng, Z., Soria, R., et al. 2020, ApJ, 900, 42 [Google Scholar]
 Lobel, A., & Blomme, R. 2008, ApJ, 678, 408 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (San Francisco: W. H. Freeman and Co.) [Google Scholar]
 Owocki, S. P., udDoula, A., Sundqvist, J. O., et al. 2016, MNRAS, 462, 3830 [NASA ADS] [CrossRef] [Google Scholar]
 Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Najarro, F., Sundqvist, J. O., & Sen, K. 2020, A&A, 642, A172 [EDP Sciences] [Google Scholar]
 Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&ARv, 21, 69 [Google Scholar]
 Rivinius, T., Baade, D., Hadrava, P., Heida, M., & Klement, R. 2020, A&A, 637, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sander, A. A. C., Hamann, W. R., Todt, H., Hainich, R., & Shenar, T. 2017, A&A, 603, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sander, A. A. C., Hamann, W. R., Todt, H., et al. 2019, A&A, 621, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
 Shenar, T., Bodensteiner, J., AbdulMasih, M., et al. 2020, A&A, 639, L6 [EDP Sciences] [Google Scholar]
 Shulyak, D., Tsymbal, V., Ryabchikova, T., Stütz, C., & Weiss, W. W. 2004, A&A, 428, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Silaj, J., Jones, C. E., Tycner, C., Sigut, T. A. A., & Smith, A. D. 2010, ApJS, 187, 228 [NASA ADS] [CrossRef] [Google Scholar]
 SimónDíaz, S., Maíz Apellániz, J., Lennon, D. J., et al. 2020, A&A, 634, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sobolev, V. V. 1960, Moving Envelopes of Stars (Cambridge: Harvard University Press) [Google Scholar]
 Sundqvist, J. O., & Puls, J. 2018, A&A, 619, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., udDoula, A., Owocki, S. P., et al. 2012, MNRAS, 423, L21 [Google Scholar]
 Tkachenko, A. 2015, A&A, 581, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tsymbal, V. 1996, in M.A.S.S., Model Atmospheres and Spectrum Synthesis, eds. S. J. Adelman, F. Kupka, & W. W. Weiss, ASP Conf. Ser., 108, 198 [Google Scholar]
Appendix A: Coordinate transformations
In this section, we derive the basic transformations of coordinates between two coordinate systems (see Fig. A.1) with orthogonal basis such that and . Any point r_{Σ1} in Σ_{1} can then be expressed by
Fig. A.1. Transformation between coordinate systems and , where t indicates the translation vector from Σ_{1} → Σ_{2}, and r_{Σ1} and r_{Σ2} denote the coordinate representations in each system for a given point. 
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, is the translation vector indicating the origin of Σ_{2} within Σ_{1}, and are the basis vectors of Σ_{2} measured within Σ_{1}.
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_{Σ2}/L_{Σ1} = s_{x} = s_{y} = s_{z}, with L_{Σ2} and L_{Σ1} describing the length scales of the individual coordinate systems. Multiplying Eq. (A.1) with , , and , we obtain in matrix form:
Thus, the transformations between the coordinate systems Σ_{1} ↔ Σ_{2} are readily given by:
Appendix B: Test calculations
In this section, we test the developed BOSS3D code by calculating synthetic line profiles of a generic resonanceline transition within a spherically symmetric wind of a prototypical hot, massive star. In order to compare the results with the singleobject code described in Sect. 2, the secondary object is assumed to be dark () and small (). 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_{orb} = 1000 km s^{−1} to the primary object, we further test if the corresponding Dopplershifts are correctly implemented for different viewing angles. The primary object is described by , , and with a wind stratification given from a βvelocity law with base velocity v_{min} = 10 km s^{−1}, terminal velocity v_{∞} = 2000 km s^{−1}, and β = 1. The required density is obtained from the continuity equation assuming a massloss rate Ṁ = 10^{−6} M_{⊙} yr^{−1}, and the source function is calculated from the Sobolev approximation (Sobolev 1960) for a pure scattering line within the twolevel approach. Finally, we specify the boundary condition by the Planck function, , if a ray hits the core of one of the objects.
Figure B.1 displays the resulting line profiles for edgeon observer directions at different phases, compared with the result from the singleobject 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.
Fig. B.1. 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 φ and a fixed inclination i = 90°, and the grey crosses correspond to the solution as obtained from the singlestar 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 singlestar code (shifted by the corresponding projected orbital velocities). 
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 lineprofile shape, the equivalent width, the radialvelocity curve and its semiamplitude). To this end, we focus on the Bstar and BHdisc scenario (model PB_SBH), and on the strippedstar and Bestar scenario with two discs (model PST_SBE02). The individual disc parameters are varied one at a time, starting from our bestfit by eye models presented in Sect. 4 and summarized in Table 3. The radialvelocity 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 Bstar and BHdisc scenario
Figure C.1 displays the results of our parameter study for the Bstar and BHdisc scenario, with variations of the four free parameters, namely the disc’s base density, , the slope of the density stratification, , the disc’s temperature , and the microturbulent velocity, . Additionally, we display the reduced χ^{2}values of the equivalent width and the radialvelocity semiamplitude as compared to the observations, as well as the combined reduced χ^{2}value (with an equal weight to both quantities).
Fig. C.1. Parameter study with base parameters defined in Table 3 for model PB_SBH. Upper left mainpanel: Synthetic line profiles at phase φ = 0 (top left subpanel), radialvelocity curve (top right subpanel), equivalent width with corresponding reduced χ^{2}value (middle left subpanels), semiamplitude of the radial velocity curve with corresponding reduced χ^{2}value (middle right subpanels), and total reduced χ^{2}value (bottom panel), all as a function of the disc’s base density, . The black dotted line indicates our bestfit by eye (see Sect. 4.1), and the blue dashed lines display the observed values with corresponding errorbars (blueshaded area). Similarly, the top right to bottom right mainpanels indicate the dependence of our measurable quantities on the disc’s temperature, , the slope of the density stratification, , and the microturbulent velocity, . 
Impact of . 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 ρ^{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 Bstar orbit (K_{1} = 53 km s^{−1}) for low base densities, and the BH orbit (K_{2} = 7.4 km s^{−1}) for high base densities, as expected. Additionally, the overestimation of the radialvelocity amplitude due to the antiphase motion of the stellar absorption profile (see AbdulMasih et al. 2020a) is clearly observed within our models at intermediate densities.
Impact of 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 (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 semiamplitude of the radial velocity curve becomes overestimated when measured from the H_{α}line wings due to the antiphase motion of the stellar absorption profile. In contrast to the lowdensity regime described above, the H_{α}line wings are mainly controlled by the disc emission at the upper temperature limit, and the radialvelocity curve is therefore following the BH orbit for the full parameter space considered here.
Impact of 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 , 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 km s^{−1}) due to a slight emission from the disc in the H_{α}line wings (probing the highvelocity regions and thus the inner parts of the disc).
Impact of . Both the strength of the line and the radial velocity curve are only mildly affected by the microturbulent velocity in the disc. The shape of the H_{α}line profile, however, switches from a relatively narrow and strongly pronounced doublepeaked feature at low to a broad and smooth one at high . When compared to observations (see Sect. 4), the latter is to be preferred.
Investigation of χ^{2}. The χ^{2}values shown in Fig. C.1 indicate that both the slope of the density stratification and the microturbulent velocity are relatively well constrained to the bestfit 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 antiphase behaviour of line core and line wings (see Fig. 2) could not be reproduced either. Moreover, at temperatures (not shown here) and , the disc’s base density would need to be increased to obtain a significant equivalent width, translating then to an Xray bright disc that has not been observed.
C.2. Parameter study for the strippedB and Be scenario with two discs
Figure C.2 displays the results of our parameter study for the strippedB and Be scenario with two discs, where the discs are described by six free parameters in total, namely the base density for each disc, and , the slope of the density stratification, and (both in a range ∈[1.5, 4]), and the microturbulent velocities, , (both in a range ∈[15, 105] km s^{−1}). As before, we additionally display the reduced χ^{2}values calculated for the equivalent width, the radial velocity semiamplitude, and the equally weighted sum of both.
Fig. C.2. As Fig. C.1, however for model PST_SBE02 and focussing on the six free parameters , , , , , and . 
Impact of and . The strength of the line increases with increasing base density of the strippedstar disc, , particularly in the lowvelocity regions (due to the overall low velocities in the strippedstar disc). At the lower end of our parameter space, the strippedstar 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 Bestar disc). Consequently, the radialvelocity curves vary significantly. For low , the H_{α}line wings are probing the Bestar orbit (including the effects of all underlying stellar absorption profiles, of course). With increasing , the radialvelocity 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 strippedstar disc.
In contrast to the effects of , the strength of the H_{α}line profile increases in both the low and high velocity regions with increasing . At the lower end of our parameter space, the line profiles consist only of the emission from the strippedstar disc (and again the stellar absorption components). Thus, the radial velocity curves as obtained from the H_{α}line wings follow the strippedstar orbit at low , while the antiphase behaviour probing the Bestar orbit is measured at high .
Impact of and . The strength of the H_{α}line profile increases with decreasing slope of the strippedstar disc’s density stratification, , in particular in the lowvelocity regions. As in Sect. C.1, increasing translates to decreasing the density of the strippedstar disc, particularly in the outer (lowvelocity) regions. Thus, the synthetic line profiles become slightly broader at the upper limit of . Since the strippedstar disc essentially becomes transparent at high , the obtained radialvelocity curve follows the Bestar orbit. Accordingly, the strippedstar orbit is being probed at low , with a significant contamination from the Bestar disc though.
Also with decreasing slope of the Bestar 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 , the Bestar disc becomes transparent and only the strippedstar disc is visible (together with the stellar absorption profiles, again). Thus, a high is favourable for measuring the strippedstar orbit, while a low would give the Bestar orbit (still contaminated by the stellar absorption profiles).
Impact of and . While the equivalent width of the H_{α}line profile is only slightly increasing with increasing microturbulent velocity of the strippedstar disc (), the shape of the line changes significantly from a doublepeaked feature to a relatively smooth and slightly broader appearance. Thus, the radialvelocity curve is probing primarily the Bestar disc at the lower limit of , and the semiamplitude is increasing with since the lineprofile wings become increasingly contaminated from the strippedstar disc emission.
With increasing microturbulent velocity of the Bestar disc () 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 doublepeaked line profile to a much smoother and broader one. Thus, the radialvelocity curves are probing primarily the Bestar disc at the upper limit of . With decreasing microturbulent velocity, the line wings are contaminated by the strippedstar disc, and the semiamplitude of the radialvelocity curve therefore becomes decreased (towards positive values).
Investigation of χ^{2}. As shown by the reduced χ^{2}values of our strippedstar and Bestar models (Fig C.2), all fit parameters are fairly well constrained to the bestfit 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 and (or and ) at the same time might give similar results as found for our bestfit model.
All Tables
Stellar parameters for the B+BH and strB+Be scenario used for our simulations of the LB1 system, and as found in the literature.
Orbital parameters for the LB1 system as calculated from our models and found in the literature.
Disc parameters for the three different (bestfit) models as obtained from our simulations of the LB1 system (see text).
All Figures
Fig. 1. Geometry used for the singlestar algorithm (top panels) and for the binary code (bottom panels). Top left panel: cylindrical coordinate system (blue) with corresponding Cartesian reference frame (red) for an observer’s direction used to obtain the synthetic spectrum of a given (single) object described in a spherical coordinate system (magenta) with corresponding Cartesian reference frame (black). Top middle panel: p − ζ (or plane) of the cylindrical coordinate system (perpendicular to the observer’s direction). Top right panel: an arbitrary slice, with 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 coordinates indicated by the red dots. Bottom left panel: two local coordinate systems Σ_{1} and Σ_{2} describing two individual objects indicated by the cyan and magenta quartercircles, respectively. Both coordinate systems are embedded within a global coordinate system Σ_{0} (typically but not necessarily chosen to be the centre of mass of both objects). The vectors t_{1} and 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 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 rays propagating through the system of object 1 (blue) and object 2 (red) potentially need to be merged. 

In the text 
Fig. 2. H_{α} line profiles at different phases φ (top row) with corresponding meansubtracted dynamical spectra (middle row) and radial velocity curves (bottom row) of the BH or strB companion. The x_{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: Bestar, model 1), PST_SBE02 (primary: stripped star, secondary: Bestar, 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, meansubtracted dynamical spectrum and radial velocity curves of LB1, 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 km 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 (phaseaveraged) equivalent width evaluated in velocity space in units of 100 km s^{−1} and the inclination used within the syntheticspectra calculations are indicated. 

In the text 
Fig. 3. Disc models in the equatorial plane corresponding to our bestfit parameters of the LB1 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 Bestar disc within the singledisc model (model PST_SBE01), respectively. For model PST_SBE02 (containing two discs), the Bestar disc is indicated by the green dotted line (partially overlapping with the red solid line from model PST_SBE01), and the strippedstar disc is indicated by the green solid line. 

In the text 
Fig. A.1. Transformation between coordinate systems and , where t indicates the translation vector from Σ_{1} → Σ_{2}, and r_{Σ1} and r_{Σ2} denote the coordinate representations in each system for a given point. 

In the text 
Fig. B.1. 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 φ and a fixed inclination i = 90°, and the grey crosses correspond to the solution as obtained from the singlestar 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 singlestar code (shifted by the corresponding projected orbital velocities). 

In the text 
Fig. C.1. Parameter study with base parameters defined in Table 3 for model PB_SBH. Upper left mainpanel: Synthetic line profiles at phase φ = 0 (top left subpanel), radialvelocity curve (top right subpanel), equivalent width with corresponding reduced χ^{2}value (middle left subpanels), semiamplitude of the radial velocity curve with corresponding reduced χ^{2}value (middle right subpanels), and total reduced χ^{2}value (bottom panel), all as a function of the disc’s base density, . The black dotted line indicates our bestfit by eye (see Sect. 4.1), and the blue dashed lines display the observed values with corresponding errorbars (blueshaded area). Similarly, the top right to bottom right mainpanels indicate the dependence of our measurable quantities on the disc’s temperature, , the slope of the density stratification, , and the microturbulent velocity, . 

In the text 
Fig. C.2. As Fig. C.1, however for model PST_SBE02 and focussing on the six free parameters , , , , , and . 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.