Autonomous Disentangling for Spectroscopic Surveys

A suite of spectroscopic surveys is producing vast sets of stellar spectra with the goal of advancing stellar physics and Galactic evolution by determining their basic physical properties. A substantial fraction of these stars are in binary systems, but almost all large-survey modeling pipelines treat them as single stars. For sets of multi-epoch spectra, spectral disentangling is a powerful technique to recover or constrain the individual components' spectra of a multiple system. So far, this approach has focused on small samples or individual objects, usually with high resolution ($R \gtrsim 10.000$) spectra and many epochs ($\gtrsim 8$). Here, we present a disentangling implementation that accounts for several aspects of few-epoch spectra from large surveys: that vast sample sizes require automatic determination of starting guesses; that some of the most extensive spectroscopic surveys have a resolution of only $\approx 2,000$; that few epochs preclude unique orbit fitting; that one needs effective regularisation of the disentangled solution to ensure resulting spectra are smooth. We describe the implementation of this code and show with simulated spectra how well spectral recovery can work for hot and cool stars at $R \approx 2000$. Moreover, we verify the code on two established binary systems, the ``Unicorn'' and ``Giraffe''. This code can serve to explore new regimes in survey disentangling in search of massive stars with massive dark companions, e.g. the $\gtrsim 200,000$ hot stars of the SDSS-V survey.


INTRODUCTION
The fact that a significant fraction of all stars or stellar remnants is in multiple stellar systems with a period of less than a few years (e.g.Sana et al. 2012;Moe & Di Stefano 2017) fundamentally affects many aspects of astrophysics.It affects the stellar evolution of both components, in some cases already during their main sequence phases, and more often in the evolved phases that will result in compact objects (such as white dwarfs, neutron stars, and black holes); it affects nucleosynthesis, the formation channels of supernovae, and the interpretation of photometric, astrometric, or spectroscopic sky surveys.And massive binaries -or their descendants -are the most prominent and frequent source of gravitational waves so far (e.g.Abbott et al. 2023).
Most of these systems cannot be spatially resolved, with projected separations that often are ≲ 1 mas.However, their orbital velocities make it possible to separate the constituents of such multiple stellar systems in velocity space, especially if spectra at different orbital phases exist.We commonly categorize spectroscopic binaries into SB1 and SB2.Here, SB1 denotes a single-lined spectroscopic binary, where only one of the component spectra is apparent in the observations, and SB2 describes a double-lined spectroscopic binary, where two sets of lines are visible in the observed spectra.The approach of using multi-epoch observations of spectroscopic binaries ★ E-mail: seeburger@mpia.deto determine the components is called spectral disentangling (e.g.Bagnuolo Jr & Gies 1991;Simon & Sturm 1994;Hadrava 1995).
In broad terms, spectral disentangling assumes that spectra of a presumed multiple stellar system -when observed at different epochs -can be described as the sum of two (or more) spectra that are invariant in their rest-frame, but whose radial velocities (RVs) change as a function time, reflecting orbital motion.The mathematical foundation of spectral disentangling has been established for 30 years (e.g.Bagnuolo Jr & Gies 1991;Simon & Sturm 1994;Hadrava 1995).End-to-end disentangling requires the simultaneous, or iterative, solution to two problems, (a) reconstructing the rest-frame spectra of each component and (b) determining the components' radial velocities at each epoch or, alternatively, the orbital solution of the overall system.If the velocities at all epochs are known, the reconstruction of the disentangled spectra reduces to a linear  2 -optimisation problem, aiming to match the combined spectra at all the different epochs.
However, the application of spectral disentangling to large data sets has some serious practical limitations.First, some literature work has assumed that (a very good guess for) various system parameters can be obtained independently (e.g.Ilĳic (2004)'s code CRES requires input of both the primary's and the secondary's velocities, shift-andadd as described in Shenar et al. (2020) and Shenar et al. (2022b) requires input of a few orbital parameters, see table 1).If the data are of limited resolution, or if the components contribute comparably to the total spectrum's absorption lines, this may not be possible.Second, there are several inherent degeneracies, foremost the fact that any featureless continuum portion of the spectrum can be assigned to either spectral component without consequences in the data match.Third, the approach works manifestly best in the regime of many epochs with data at high spectral resolution (compared to the orbital velocity changes) and at very high signal-to-noise, so that the dimmer component causes distinct changes in the combined spectrum.Finally, solving the full non-linear problem, i.e. optimising the model-data match over all possible primary velocities, mass ratios, and disentangled spectra is very time-consuming.
Over the last decade and for the next decade, vast spectral surveys are driving an exponential growth in the number of high-quality stellar spectra.Current or upcoming surveys include SDSS (York et al. 2000;Kollmeier et al. 2019), LAMOST (Luo et al. 2015), DESI (DESI Collaboration et al. 2023), WEAVE (Dalton et al. 2012) and 4MOST (De Jong et al. 2019).Further, the Gaia Data Release 4 will provide a vest set of spectra.Many of these surveys, in particular Gaia (Gaia Collaboration et al. 2021) and SDSS-V (Kollmeier et al. 2019) have multi-epoch observations scheduled across the entire sky.
These surveys offer vast potential to map the stellar binary population via spectral disentangling.Some surveys, such as SDSS-V have explicit programs to systematically survey stars searching for massive dark companions; there spectral disentangling is crucial to identify "contaminants" with two luminous components (e.g.Shenar et al. 2022b,a;Mahy et al. 2022).
In this context of vast spectral surveys, new requirements -or desiderata -arise for practical approaches to spectral disentangling.The approach must be • fast, so that multi-epoch data for 10 4 to 10 6 systems can be analyzed.
• autonomous, in the sense that initial parameter guesses that permit sensible solutions must be found algorithmically and reliably.
• robust, given that large surveys typically have fewer epochs, less S/N and lower resolution than single-object studies • user-friendly, as a wide community should be in a position to consistenly analyze different data sets, or reanalyze a given data sets with different constraints on acceptable solutions.
In this paper, we propose a new implementation of spectral disentangling, designed to address these issues.It wraps the process of finding starting guesses, solving for the disentangled spectra, and optimising the flux-and mass ratio parameters into one continuous pipeline written in Python.The approach is also fast enough that it can be applied to surveys of many thousands of objects.We include features such as regularisation to ensure desired properties in the spectra, as well as other adaptations on the original method to optimise the code for survey disentangling.
Due to Python's rise to popularity, we have elected to write this implementation of the method in this language.While at the surface level, speed might be a concern, many of Python's modules are partially written in compiled languages (such as C) and merely provide an interface familiar to the average Python user.Thus, this issue remains manageable.Python does have the advantage of being widely known in the scientific community, making an eventual release of the code as a package accessible to many.

SPECTRAL DISENTANGLING METHODOLOGY
We first briefly review the existing method and codes in spectral disentangling in 2.1, then contrast disentangling as an approach with other methods in the literature in 2.2.We lay out the desired characteristics of a code for survey disentangling in 2.3.In subsection 2.4 we describe the prepocessing of the data.We explain the two-step process of optimising for the velocities and mass ratio in 2.5, and how we obtain the systemic velocity and light ratio in 2.6.
These codes employ one of two major methods for solving the linear disentanging problem, either in fourier space (KOREL, fd3) or in wavelength space (CRES, Spectangular).More discussion on this can be found in e.g.Ilĳic et al. (2004), but in summary, both methods come with their own advantages and disadvantages.The Fourier method, most notably, outperforms the wavelength-based alogorithm in terms of speed and thus allows for the implementation of further generalisations (Hadrava 2009).However, this comes at the cost of requiring all spectra to be sampled on the same grid, as well as giving each point the same weight, and implicitly assuming the resulting spectrum to be a periodic function of the wavelength (Sablowski & Weber 2017).Performing the disentangling in wavelength space allows the user to extend the wavelength range on which the component spectra are computed based on the RV shifts of the individual epochs, as well as being less vulnerable to edge effects arising from the Fourier method.

Disentangling vs Spectral Model Fitting
Spectral disentangling entails the "non-parametric" reconstruction of the individual components' rest-frame spectra.Alternatively, one can view the whole problem as a forward-modelling problem, drawing on a set of stellar templates.This has been explored and implemented by several groups (e.g.Traven et al. 2020).An important downside of these methods is the fact that spectra in close binaries often do not look like simple, single-star spectra: they may be "exotic" objects such as stripped stars, they may show exceptionally fast rotation or may show emission lines from decretion disks.It can thus be easy to miss the signatures of "strange" companions (e.g.Jayasinghe et al. 2021Jayasinghe et al. , 2022;;Shenar et al. 2020;El-Badry et al. 2022;Bodensteiner et al. 2020;Frost et al. 2022) when only considering a limited set of stellar templates.
This highlights the strength of template-independent methods, such as direct subtraction (Ferluga et al. 1997) iterative subtraction, also known as shift & add (González & Levato 2006, implemented in Shenar et al. (2020, 2022b)), or, as described here, disentangling (Simon & Sturm 1994).By not having to pre-select a template, we remain flexible to a range of potential outcomes of the procedure.

A Disentangling Approach for Large Spectral Surveys
As mentioned in the Introduction, practical approaches to disentangling for vast spectral surveys call for a code that is fast, autonomous with respect to the initial parameter guesses, astrophysically flexible to accommodate the unusual spectra of close binaries, and robust  1.This table summarises a number of properties of some prominent disentangling codes, which include the author(s) who wrote the codes, the programming language they are written in, whether disentangling takes place in wavelength or fourier space, whether they solve only for the spectra or also attempt to find the RVs and/or Orbital solutions, whether the code is equipped to handle a third component, and which input or guesses are required.Step i) Step ii) Step iii) with respect to suboptimal numbers of epochs, S/N and spectral resolution.
Here we set out to devise, test and verify such an approach.The end-to-end approach can be divided into several steps, which are conceptually illustrated in Figure 1, which serves as a guide and schematic representation of the process.
We want to start with multi-epoch spectra for any given object, and at the end have both an estimate of the (usually two) disentangled spectra of the primary component's velocities at each epoch, the systemic velocity, the component mass ratio, and the mean component flux ratio.To achieve this for any one object, there are essentially three stages: (i) The spectra at all epochs for any one object have to be consistently pre-processed (normalisation, wevelength rebinning); and initial guesses for the velocities of the primary spectral component (defined as the component with the most prominent spectral lines) need to be derived algorithmically.
(ii) The optimisation of the main parameters consists of two parts.
A non-linear part -trying to find the best primary velocities at all epochs, ì   and the mass ratio .And, a linear part, determining the two disentangled component spectra, ì   and ì   that best match the multi-epoch spectra ì for each assumed (ì   , ).(iii) The resulting spectra ì   and ì   have two remaining problems: they are pure mathematical constructs and know nothing about the rest-frame wavelengths of features.And any (modestly small) constant can be subtracted from ì   and added to ì   , leaving the data match to the ì   .The last step then uses template spectra to a) fix the rest-frame of the disentangled spectra (or the systemic velocity of the system) and fix the luminosity ratio of the two components, by requiring physically sensible absorption line equivalent widths.
We now describe these three stages in turn.

Preprocessing of the Spectra
Before disentangling, it proves convenient and useful to pre-process the observed multi-epoch spectra to simplify the subsequent math, labelled as step (i) above.This first entails masking, or interpolating over, bad pixels.Then resampling the spectra to a wavelength grid that is uniformly sampled in ln (hereafter abbreviated as Λ), which linearises velocity shifts.Finally, we normalise the spectra by dividing them by a running median filter, where the filter must be chosen to be much wider than individual spectral features.The point of this is not necessarily to remove the "continuum" of the stellar spectrum, which is often conceptually and practically poorly defined, especially in cool stars with many spectral lines.We are much more concerned with removing the low-frequency variations in the observed spectra, as they may be dominated by instrumental effects.For our disentangling it is most important that we do this consistently across all observed epochs.We also apply the same median filtering to all the templates, i.e. we remove the low-frequency variations in the model spectra fully consistently.Other methods, such as polynomial fitting and clipping were considered, but ultimately median filtering was selected as the method of choice due to its speed, ease of use, and consistency across spectral types.
From this normalised spectrum, we then subtract 1 at all pixels, as this further simplifies the subsequent analysis.The different epochs are then concatenated into one long vector of length    •   (number of epochs times number of pixels).We call this vector ì   .As the next part of pre-processing, we must find initial guesses for the primary component's radial velocities for each epoch  , ì   ), to aid the convergence of the non-linear parameter optimisation of step (ii).We have implemented two ways to obtain these initial guesses, either using template cross-correlation, or using the TIRAVEL algorithm (Zucker & Mazeh 2006, described in the appendix).For the cross-correlation, we first construct a grid of (rest-frame) template spectra from Kurucz (1979).These templates cover a suitable range of effective temperatures (∼ 20, 2.7kK to 25kK), surface gravities (∼ 5, -0.35 to 5.4), and two different rotation velocities sin() (10 and 100 km/s) and are matched in resolution and wavelength sampling to the observed spectra.For any given epoch we perform a cross-correlation (CC) between the observed spectrum and all templates, and take the best template to be the one that yields the highest CC peak, when averaged over all epochs.The position of the CC peak for the best template yields our starting guess for the primary velocity at that epoch, yielding the initial ì   .As starting guesses for the mass ratio, we adopt     ≡  ≡ 1.In the system's center of mass frame, the primary and secondary velocities are related via ì   ≡ − × ì   .But the ì   we derive for each epoch are in the barycentric frame and we do not know the actual systemic velocity,   .However, we can for the subsequent steps simply assume that   ≡ 0. This will then yield disentangled spectra of the correct shape, just in an ill-defined velocity reference frame.However, this can be remedied by correlation with template spectra in a subsequent step, as we show below.

Parameter optimisation and initial disentangling
We can now proceed to step (ii), indicated in the second box (from the top) of the graphical representation in Figure 1.As mentioned, this step consists of two aspects.The linear part estimates the disentangled spectra that best match the observations at all epochs for given ì   and .Each execution of this linear disentangling step yields  2 (ì   , ).This step is wrapped in a nonlinear parameter optimiser (we settled on the Nelder & Mead (1965) method) that then finds most likely parameter estimates ì    ,   as (1)

The linear disentangling step
It is more sensible to describe the linear disentangling step first.We restrict ourselves to the binary (rather than triplet, etc.) case, where we seek the initially disentangled components ì   and ì   , specified on a logarithmic wavelength grid with the same ΔΛ as the observed spectra: ì  / ≡ { / (Λ  )} with  = [0,   ].In this context, non-relativistic velocity changes correspond to constant shifts in ln  and can be specified as In this first step we will retrieve 'initial' disentangled spectra in an ill-determined velocity frame rather than in the physical rest-frame.
We will remedy this in a subsequent step (Section 2.6).For now, we can assume formally that the systemic velocity is not only constant but also zero; then the two component's velocity shifts at epoch  are related via ΔΛ   = −ΔΛ   /.In this case we only need to know the velocities for one component (say, A) and the system's mass ratio to determine all relevant velocities.
The shifted spectra of the components for each epoch  can be written as: The predicted composite spectrum ì  ,   at epoch  is then given by: The predictions for the composite spectra at all epochs  can be cast more elegantly into a matrix of form: where the column vector ì  is the concatenation of ì   and ì   , and ì    the concatenation of all ì  ,   at all different epochs    .Thus, ì  is a column vector of length 2  , and ì    a column vector of length   •    .The matrix M is then a matrix of dimensions (  •    ) × 2  made of 2 ×    blocks (one for each component and each epoch) of dimensions (  ×   ).M is a sparse matrix, with the only nonzero elements being (off-)diagonals, whose position and value are governed by the per-epoch shifts of both components, ΔΛ /  .Figure 2 shows a simplified schematic of Equation 6, and a more detailed description of the matrix structure can be found in Appendix A3.
For a given set of ΔΛ   and q, we then wish to find the ì   for which ì    best matches -in the L2-norm or ℎ 2 sense -the concatenated set of multi-epoch observations ì : As the matrix M is sparse and frequently very large, it makes sense to look for efficient methods for the use case.Simon & Sturm (1994) propose the use of Singular Value Decomposition (SVD) (e.g.Forsythe et al. 1977) to solve for ì .However, we found that the LSMR algorithm by Fong & Saunders (2010) provided a more efficient iterative algorithm that exploits the sparsity of the matrix to arrive at a solution quickly.The resulting ì   can then be separated into ì    and ì    .
[H] x A and x B optimization for a given set of parameters Figure 2. A cartoon of the setup of the linear algebra portion of the disentangling scheme.The top portion of matrix M on the left-hand side consists of (off-) diagonals whose position relates to the per-epoch velocity shifts.The vector in the center is simply the component spectra stacked on top of each other and the desired output of the linear algebra procedure.The vector on the right is the individual epoch composite spectra, stacked as well.The regularisation scheme detailed at the bottom of the cartoon is described in more detail in section 2.5.2 Figure 3 shows an example of the disentangling result for a mainsequence binary consisting of a 1.1  ⊙ primary with a 0.95  ⊙ secondary, i.e. a mass ratio of  = 0.86.The velocity semi-amplitude of the primary is 200 km/s, we have assumed a circular, edge-on orbit and sampled the velocities regularly for 6 epochs at a resolution of  = 2000 and with a signal-to-noise ratio of 30.The systemic velocity is 100 km/s.Disentangling has been performed assuming ΔΛ   and  are known.Both the composite spectra (black and grey lines) and the individual disentangled component spectra are in excellent agreement with the observations (top) and input spectra (two panels below).

Regularisation
While the method described so far can yield acceptable results, the two disentangled spectra are often not smooth.This is due to the nonuniqueness of the solutions: especially for sections of the continuum with no/few lines, it is ill-defined how much each component spectrum is specifically contributing to the composite.
Tikhonov regularisation (Tikhonov 1963;Phillips 1962) introduces an additional term in the least-squares minimisation (Equation 7), which can be used to penalise certain undesired features in the solution, and encourage desired properties.So, instead of minimising the expression in Equation 7, we seek to minimise: where M, ì  and ì   are as previously described.L is a matrix by which we regularise the solution for ì  and  expresses the weight of this regularisation compared to the disentangling problem.There is a trade-off: very small  will lead to very little regularisation (and  = 0 reduces to the original problem posed), while very large  will over-regularise the solution, leading to a loss of the original features.
This can be recast into a form very similar to Eq.7: Here, L is of shape (2  × 2  ), and ì 0 is of length 2  .Thus, this extension to the original problem allows regularisation of ì  by simply appending to the matrix M and the vector ì   .We want matrix L push for smooth ì  solutions but without enforcing a specific shape onto the spectra.We do this by minimising the second derivative or curvature among adjacent elements of ì   .In practice, the matrix L has 1's along the diagonal , and -1/2 along the two first off-diagonals above and below the main, as illustrated in the scheme shown in Figure 2.This matrix leads to a norm that is simply the sum of the second derivatives among all sets of adjacent 3 pixels (expressed as a finite difference).Therefore, it penalizes strong local curvature of the result, and thus "jaggedness" of the disentangled spectra.
In Figure 4 we compare the result of different regularisation strengths, ranging from none (top) to a good amount (middle) to too much (bottom).We see that when no regularisation is applied, the disentangled spectrum is jagged due to the degeneracies.With an appropriate amount of regularisation, this can be remedied.If too much regularisation is applied, features and details will start to be "washed out" -it is thus important to choose the regularisation parameter correctly, e.g. by requiring it not increase the least-squares or  2 significantly.In practice, we seek a regularised solution that fits the data to within Δ 2 ∼ 3 as well as the unregularised solution.We have found a suitable  for each (simulated) dataset by trial and error; automatisation of this is an endeavour for future work.

The non-linear optimisation step
So far, we have shown how to robustly solve the disentangling problem if the RVs ì   for all epochs  and the mass ratio  are known.We now need to implement a robust way to find the optimum parameters ì   and  such that the two constructed component spectra ì   and ì   , shifted by their respective epoch velocities, optimally reconstruct the observed composite spectrum ì   for all epochs.To this end, we define the residual between the observed and reconstructed spectra as: or, more explicitly, Equation 11 explicitly recognises that the residuals are a function of ΔΛ /  .As the next step, we have to find the best values for these ΔΛ /  by minimising the value of  2 as a function of  and ì   , defined as: where cov is the covariance matrix of the observed spectra.
We do this minimisation via the Downhill Simplex Optimisation by Nelder & Mead (1965), specifically the scipy.optimiseimplementation.
Robust optimisation requires a sensible initial starting guess.We initialise it with ì   and  as described in 2.4.Then the optimiser solves at each step the linear problem (Equation 9) for the current values of ì   and .This returns the component spectra ì   and ì   and from them ì    , which yields  2 from Equation 12 from the comparison with the data ì   .Eventually, we expect the optimiser to return the best values for the RVs and mass ratio, those that allowed for the most accurate reconstruction of the observed spectra.

Determining the systemic velocity and light ratio
Mathematically, spectrum disentangling is possible while only knowing the relative per-epoch RVs of the primary and the mass ratio, as has been demonstrated.However, this yields two disentangled spectra that are in an unspecified velocity frame that is neither the center-ofmass frame nor the rest-frame; and the velocity frames for the two components will generally not be the same.But in practice, we are interested in the systemic velocity   and the light ratio of the two components .
We now describe how we find these parameters and move the disentangled spectra to a well-defined rest-frame, as outlined in step (iii) of Figure 1.In doing so, we basically draw on external astrophysical information: the rest-frame wavelengths and (approximate) equivalent widths of prominent stellar absorption lines in the disentangled spectra must resemble -at least roughly -those in a comprehensive set of template spectra.
We start by using such template spectra to find the systemic veloc-ity.We do this by computing the cross-correlation (CC) of the initial solution vectors ì  /  with spectra in the template set and identifing the best template.We denote the velocity of that CC peak as  / , which describes the offset of the velocity found by the optimiser from the true observed velocity of the primary.In this notation we have where ṽ  are the velocities in the initial ill-defined frame,   the true velocities in the observer frame, and   are the centerof-mass-frame velocities.With starting guesses via template crosscorrelation, we expect   to be close to zero.This is because the templates are in the rest frame, so initial velocities found using them should also be in the rest frame, meaning the offset   should be small.TIRAVEL (Zucker & Mazeh 2006) determines initial guesses by cross-correlating the different epochs with each other and solving for the most probable vector of RVs.Then,   could be large, as we have no prior information about absolute RVs, only relative (to other epochs).By contruction from the disentangling, the primary and secondary velocities in our ill-defined frame are related by ṽ = − ṽ  / holds.And, we have the relations with nomenclature analogous to Equation 13; Eq. 14 simply states that the true observed velocity is the secondary's center-of-mass frame velocity plus the systemic velocity, which in turn is related to the true velocity of the primary via the mass ratio .We can now use these offsets  / and Equations 13 and all 14 to determine the systemic velocity of the system: Note that even when   is small,   might be large, especially when the originally assumed systemic velocity is far from the truth.This puts us in a position to solve for the component spectra in the rest-frame using these results and Equation 13.This process also gives us templates ì   and ì   , the templates that had the highest CC peak with the respective disentangled component spectra.
We can now turn to determining the light ratio , using the bestfitting templates.In the end we attribute different portions of the "featureless continuum" to the two disentangled components so that their equivalent widths are physically plausible.We obtain  by minimising arg min which scales the templates for both components, ì   and ì   , such that the scaled templates most resemble the disentangled spectra.
An example of the results of this 3-step (preprocessing, optimising ì   and , finding   and ) process can be seen in Figure 5, where the accuracy of each parameter as found in the optimisation or subsequent step, is explored.The system has the same parameters as the one discussed in section 2.5.1 whose disentangled spectrum can be seen in Figure 3. of the light ratio.We selected these scaling of the different parameters for better visualisation.We see good recovery for all parameters, except two of the RVs.This is, however, not further surprising, as these RVs are fairly close to the systemic velocity, meaning there is only a small doppler shift of the two spectra relative to each other, making determination of the accurate velocities difficult at these low resolutions.

ALGORITHM VALIDATION
We now proceed to validate this autonomous multistep approach to spectral disentangling that we have laid out.In particular, we want to explore under which circumstances it yields sensibly disentangled spectra and reasonable physical parameters.This will depend on both the physical properties of the binary system (velocity semiamplitude , mass and light ratios, effective temperatures) and on the observational set-up (S/N, spectral resolution, number of epochs).
In exploring this, we will particularly focus on the observational parameter regime pertinent to large spectroscopic surveys: modest resolution and S/N along with relatively few epochs.We do this by simulating composite spectra of binaries with both hot and cool primaries, and then running the algorithm 'blindly', or autonomously, on them.These simulations allow us to assess both the robust and the problematic disentangling regimes.After these simulations, we illustrate the approach using a well-established binary system.

Simulated Data
We perform an initial test of the method on simulated data, given that it allows us to precisely control the system parameters and check how well the algorithm recovers them.To simulate co-evolving binary components, we use a 1 Gyr and 250 Myr isochrone from Bressan et al. (2012) for the cool-and hot-star primary simulations, respectively, selecting a 1.1  ⊙ and 3  ⊙ star as the primary.
For resolutions of  ≈ 2, 000 and  ≈ 20, 000, we then explore a grid of different RV semi-amplitudes , linearly spaced from 50 km/s to 250 km/s, and light ratios .For , we consider logarithmically spaced values from 0.01 to 1 to explore the more "extreme" regimes of very faint secondaries, as well as a linearly spaced grid from 0.1 to 0.9.With the primary's parameters (mass, age,  eff and ()) and , we can use the isochrones to get analogous parameters for the secondary; the two components' masses then also yield . We then select spectral templates, using Kurucz (1979) spectra for the two components.Then, we set the orbital parameters of the system (assuming a circular, edge-on orbit), sample the RVs of both components at 6 different epochs uniformly in phase-space over half an orbit, and create composite spectra, adding noise.Here, we assumed a signal-to-noise ratio of 30.
We feed these simulated spectra to the disentangling and optimisation algorithm, assessing how well it is able to obtain the correct system parameters, reconstruct the observed composite spectra ì   and solve disentangled rest-frame spectra ì   0 and ì   0 .We assess the quality of the parameter through this figure of merit: which is bound between 0 and 1, where  ,pred are the estimates returned by our algorithm,  ,true the true (simulation input) values, and   are the normalised weights assigned to each parameter.We have assigned weights of 0.25 to ,   and , and a weight of 0.25    to each of the primary RVs.The weighting avoids the FOM being dominated by the accuracy of the recovered primary RVs (especially in the case of many epochs), which are comparatively easy to find, and places greater importance in correctly determining the parameters associated with the secondary.It is pertinent to note, at this point, that a "good" FOM means that the system parameters have been recovered accurately, not necessarily that the disentangled spectra are equivalent to the "ground truth".It is simply a metric to assess the performance of the optimiser.
The results of this validation are shown in Figure 6 for cool (left panel) and hot primaries (right panel).The different subpanels show a grid in velocity semi-amplitude  and light ration .The background color of each subpanel indicates the FOM, with lighter colors indicating better parameter retrieval.For high velocity semi-amplitudes (compared to the spectral resolution) and distinctly different component temperatures (whenever  ≪ 1) the parameter retrieval is robust and quite precise; this is particularly true if the primary star is already cool.But the Figure also shows that there are two regimes where our algorithm struggles: first, for a mass ratio of unity (twin star spectra) there is a degeneracy between the two components' velocities.Consequently, the optimiser cannot determine the other parameters correctly.For the very faint secondary regime ( = 0.01), the optimiser finds the correct primary velocities, but has some difficulties with the mass and light ratio, as well as q.This is not particularly surprising, as in this regime the signal of the secondary is on the level of the noise, and thus finding parameters that pertain not just to the primary (as the RVs do, in this case) is problematic.
Second, Fig. 6 also shows that lower velocity semi-amplitudes (below the spectral resolution) lead to a worse determination of the individual component's velocities, and consequently of the other parameters.For the higher-mass and hotter primary, the broader spectral lines and comparatively similar spectra between the primary and secondary also lead to an underestimation of the RVs, as the CFF peak is substantially affected by the secondary.
Figure 7 is analogous to Figure 6, with  varying from 0.1 to 0.9 in steps of 0.2.Again we see that for both the cool and the hotter primary, our disentangler has difficulties correctly recovering the velocities for the lowest RV semi-amplitudes.As we are now excluding the more extreme cases of twin stars and extreme light ratios, where the secondary contribution is on the level of the noise, we observe generally higher (and thus better) figure-of-merit values.With the linear grid in  we can see a more gradual trend in the robustness of our disentangling approach, with accuracies getting worse as we get "too close" to an equal mass binary, especially for the 3M ⊙ primary, as the spectra of the primary and secondary show rather similar features.
Figure 8 shows individual disentangled spectra for different simulated systems.We show a smaller wavelength range for the cooler primary, as there are more and narrower lines, whereas for the hotter primary, the lines are wider and fewer.We see that the disentangler successfully recovers the primary component in all cases, using both the "true" parameters as well as the ones found by the optimiser.For the smallest light ratio of 0.01, both in the hot-and the cool-star case, the disentangler struggles to recover the secondary, due to its small contribution to the spectrum.There are still some issues with the secondary spectrum recovery for the hotter primary even at a more moderate light ratio of 0.1.We believe this to be due to the very similar lines between the primary and secondary -for moderate velocity shifts, two similar spectra being shifted against each other look like the lines are "widening" and "narrowing", rather than fully seperating, which causes some issues in the disentangling process.

Disentangling SB1 Systems with Dark Companions
Part of the motivation for this work is to find systems where the secondary is dark, i.e. does not show up in the disentangling.To explore how our code responds to such systems, we have simulated a binary consisting of a 1.1  ⊙ primary with a dark companion (a so-called SB1 system).Like before, we set a velocity semi-amplitude of 200 km/s for the bright component and assume a circular, edge-on orbit.RVs are sampled uniformly in time at a spectral resolution of  = 2000 and a signal-to-noise ratio of 30.
We then applied our optimiser and disentangling procedure to this system in the same way as the previous simulations, presuming we have no prior knowledge of its SB1 nature.Inevitably, the optimiser will find some   and , which must be spurious.And then it determines an  from the disentangled spectra that were computed using these   and .We would expect the secondary spectrum to be essentially noise, and the resulting  to be small, presumably an upper limit.
Figure 9 shows the results of the optimiser in one realization of this scenario.We see that the primary velocities are found near perfectly.In the right panels, we show the formal results for the three problematic parameters.As expected, the best-fit   and  vary greatly among different noise-realizations of the mock data, after being initialized as before.
However,  is always found to be small.This can be traced back to the fact that the disentangled spectrum of the (non-existing) "secondary" is close to noise, as illustrated in Figure 10, which shows the disentangling results using the parameters as determined by the optimiser.The primary spectrum is of course recovered well.The computed secondary spectrum essentially appears to be noise, which then leads to the small recovered .
This shows that our approach degenerates gracefully for an SB1 The colour of each parameter indicates how well it was recovered within its specific system, with pink indicating a comparatively large distance between truth and fitted parameter, and cyan indicating a small distance, and thus good agreement.Lastly, the size of each symbol is related to its variance in the bootstrapping process, with a large symbol indicating a large variance (and thus a low precision) and a small symbol indicating a small variance and higher precision.For each panel, the truth is plotted in black, with the disentangling result achieved by using the "correct" input parameters in green, and the actual result cyan the autonomous disentangler in magenta.The disentangler solutions are scaled by the light ratio (true value for cyan, value recovered by the optimiser for magenta) to match them to the true input.Thus, a "correctly shaped" spectrum that is scaled incorrectly indicates an incorrect light ratio, while a spectrum that is offset relative to the truth indicates that the systemic velocity was not found correctly.For the most extreme light ratios (top row), even with "correct" input parameters, the disentangler struggles to recover the correct spectrum of the secondary.The optimiser also has difficulties recovering the correct spectrum in the higher-mass case, owing to the wide lines found in hot stars.system towards an ill-defined   and , with high variance or uncertainty, and towards a small estimate of .Further, inspection of the actual disentangled component spectra reveals a lack of features in the secondary, in line with the input of a zero secondary.The physically most sensible upper limit on  in this scenario depends on the template and its spectral features, effectively   ( eff ,v sin().If one has external priors e.g. on the temperature, this can be easily incorporated into the   estimate.

Sample Application: the "Unicorn" & "Giraffe"
The two systems, V723 Mon ("The Unicorn") & 2M04123153+6738486 ("The Giraffe"), present excellent examples of the power of disentangling compared to cross-correlation and other template-based methods.In an initial study, Jayasinghe et al. (2021) and Jayasinghe et al. (2022) had reported to have found SB1 binaries with a primary mass and orbital parameters requiring the secondary to be a dormant black hole.Further study by El-Badry et al. (2022) (hereafter E22) found that the authors of the first study had been led to incorrect inferences about the nature and mass of the primary, whihc then led to an incorrectly inferred a mass limit for the secondary.In E22, spectral disentangling revealed a non-degenerate stellar secondary in both cases.Both systems were products of previously occurring mass transfer, leading to a stripped, luminous, but very low-mass primary "masquerading" as a much higher mass star, and a higher mass secondary.In the case of the Unicorn, the secondary had also been significantly spun up, smearing out its spectral lines.
Disentangling was able to reconstruct the component spectra and help identify the stars, but took significant attention to detail in the analysis, as well as prior knowledge of the velocities of the primary, and solid guesses about further system parameters.
In this work, we present "blind" disentangling of the Unicorn (with the Giraffe included in the Appendix), assuming no prior knowledge of the systems other than the assumption that they are binaries, and thus disentangling is expected to yield sensible results.We also explore the performance of the disentangler assuming the spectra had been more "large survey style", i.e. at lower resolution.We apply our method to the Unicorn here, as it is the more "complicated" of the two systems to identify, due to the very rapid rotation of the secondary.
Both in E22 and in this work, we use data from the Keck/HIRES spectrograph (Vogt et al. 1994) to identify the components of the systems.Natively, the data have resolution of  ≈ 60, 000, with 7 epochs for the Unicorn, and 8 for the Giraffe.The data cover a wavelength range from 3900 to 8000 Å, and have a typical SNR of 20 per pixel at 5000 Å.More details can be found in Jayasinghe et al. (2021) and Jayasinghe et al. (2022) for the Unicorn and Giraffe data, respectively.
Figure 11 shows how well the optimiser was able to recover the parameters of the Unicorn system at the native resolution of the data,  ≈ 60, 000.This was performed on the wavelength range from 5275 to 5370 Å, as this section contains enough lines to successfully disentangle while still being narrow enough for the assumption of a constant light ratio not to break down.We see the optimiser seems to recover the individual velocities of the primary very well, owing to the clear and dominant lines in the composite spectrum.It struggles slightly more with the systemic velocity, and light ratio, likely because of the very rapid rotation and thus broadened lines of the secondary.Despite this, the light ratio is still found to an acceptable accuracy.This is likely due to a combination of the light ratio finder using both the primary and secondary spectra, and the set of templates used in the process containing non-rotating and rapidly rotating (v = 100km/s) templates.
In contrast to this, we see Figure 12, which shows the performance of the optimiser for the same data, now artificially resampled to a lower resolution of  ≈ 2, 000.A larger wavelength window of 5200 to 5450 Å was chosen here for disentangling to account for the more broadened and thus fewer lines.We see here that the RVs tend to be underestimated, with the disentangler finding generally smaller velocities than E22.This is unsurprising, as the lower resolution leads to broader CC peaks, where eventually the peaks of the primary and secondary velocities are so broad that the secondary peak "drags" the primary down.Because of this underestimation in velocities and the broadened lines of the secondary, the optimiser also cannot correctly assess the mass ratio.However, we do see that the optimiser has pushed q away from its initial starting guess of () = 0 towards higher values, meaning the data do tell us that the secondary is heavier than the primary, but struggles to perfectly determine how much heavier.As a consequence of this, the disentangled spectra are less accurate, also making it difficult to find the correct light ratio.
In Figure 13, we see the disentangled Unicorn spectra for the native resolution of  ≈ 60, 000, as well as the models from E22.Here, we see generally good agreement between the two, with one of the more apparent differences stemming from the slightly different light ratios between their work and this one.For example, we see that the E22 model for the secondary generally shows slightly "stronger" lines, i.e. a bigger difference between maxima and minima.There does not seem to be a significant offset along the wavelength axis, suggesting that the systemic velocity was recovered (mostly) accurately, at least to within the resolution limit.
Figure 14 shows the disentangling result for the  ≈ 2, 000 case, zoomed out more to cover a larger wavelength range than Figure 13.We elected to do this due to the "smearing out" at lower resolutions that removes many of the finer features and broadens the dominant lines.To still show a good number of features, a wider wavelength window is necessary.Despite the underestimations of the velocities and mass ratio, the disentangled component spectra still resemble the (downsampled) models from E22.Many of the visible features are recovered, even though the underestimated velocities likely lead to some of the features of the primary being more "washed out".Sim-Figure 10.Formal disentangling result shown for one epoch for a simulated SB1 system, computed using the best-fit parameters found by the optimiser, analogous to figure 3. The spectrum of the primary is recovered, as it is the only luminous component in the system.The code also constructs of course a second component, which is close to just noise, as it should be for a dark secondary.
ilarly, these in conjunction with the incorrect  have an analogous effect on the recovered spectrum of the secondary.The issues with the light ratio then arise naturally as a consequence of the recovered spectra being more "smeared out" than the templates we are comparing to.Additionally, for larger wavelength window, the assumption of a constant light ratio across the spectrum holds less true, which is likely a contributing factor here.

DISCUSSION AND SUMMARY
In this work, we have set out to implement spectral disentangling as described by Simon & Sturm (1994) and Hadrava (1995), among others, in an algorithm tailored for million-stars surveys.This meant adapting the process to be suitable for the regime of few epochs, modest S/N, and moderate resolution ( ≈ 2, 000).At the same time, we required the code to be robust in a range of scenarios (e.g. both hot and cool stars with a range of companions), as well as able to autonomously find the parameters necessary for the disentangling procedure precisely and accurately.Given the goal of applying it to large volumes of data, the code also needed to be fast.
Our approach has been to combine a number of known aspects and strategies to create an end-to-end autonomous pipeline that fulfills all of the above criteria, while consisting of a number of internal steps.We have implemented wavelength-space disentangling based on Simon & Sturm (1994); Tikhonov regularisation (Tikhonov 1963;Phillips 1962) to ensure smoothness of the disentangled spectra; downhill-simplex optimisation (Nelder & Mead 1965) to optimise parameters; template cross-correlation and TIRAVEL (Zucker & Mazeh 2006) to get initial velocity guesses; as well as templatebased methods to find the systemic velocity and light ratio of the system, all in one.We have explored the efficacy of this process using a variety of synthetic SB2 spectra that emulate those provided by the large surveys we developed this pipeline for.
We have found overall satisfactory performance on our simulated systems, as well as a real-life example (the Unicorn), which had previously been misidentified (Jayasinghe et al. 2021), and required  great care to analyse correctly (E22).We have shown that in the target regime of few epochs, modest S/N and moderate resolution, the algorithm can still recover the relevant system parameters and the component spectra for many of the systems explored.However, there are still areas of the parameterspace where our method fails or encounters issues: • Extreme light ratio, relative to S/N.For light ratios below  ∼ 0.1, while the algorithm generally still recovers the primary and associated velocities correctly, the spectrum of the secondary is overwhelmed by the noise, and thus both the spectrum of the sec-ondary and the mass ratio  cannot be determined correctly.In this case, the disentangler reacts similarly as if the system were an SB1, meaning we can set a light ratio "cut-off", below which we can no longer separate an SB1 from an SB2.
• Low velocity semi-amplitude, relative to resolution.If the velocity semi-amplitude (especially of the primary) falls close to or below the minimum resolvable velocity, the first step of the optimisation, finding starting guesses for the velocities of the primary, fails, and subsequent steps cannot rectify the issue.Here, the majority of parameters are determined incorrectly, and thus the disentangled spectra are also spurious.
• Equal-mass binaries.For systems where the primary and secondary are of equal mass, or close to it, the associated spectra look too similar, introducing an ambiguity as to which set of lines belongs to which object at which epoch.In this case, the velocities for the primary and secondary are degenerate for any one epoch (as the method cannot tell "which is which"), causing further issues with the disentangler.One method to rectify this might be to invoke orbital fitting, however this would only be possible in the regime of sufficient epochs.Careful, individual treatment of these systems might provide another path to a successful solution, but goes against the spirit of this work, and autonomy of the process.However, despite the issues with finding the correct parameters, the disentangled spectra often still resemble the truth quite well.This is because a simple switch of the velocities in the matrix solver step simply means that the spectra are swapped for the affected epoch.If the spectra are the same or very similar, this has only a minor effect on the disentangling result.
• Hot stars.Due to the relatively few and broad lines of hot stars compared to cooler stars, the disentangler encounters more issues for these objects.The broader lines lead to a less precise determination of radial velocities, which creates problems when trying to find the mass ratio, spectra, and other parameters.In many cases, the disentangler does still arrive at a satisfactory solution, but less consistently than for a similar system in terms of light ratio and velocity semi-amplitude but with a cooler primary.
One of the novelties of this work is the inclusion of Tikhonov regularisation (Tikhonov 1963;Phillips 1962) to ensure that spectra are smooth.This process allows us to remove much of the highfrequency noise from the solution, much like a truncated SVD would, but it is independent of the solver method employed and allows us to use fast, iterative algorithms to solve for the spectra while retaining desired characteristics.It also does not enforce physical constraints on the solutions, allowing a great deal of "freedom" when determining the component spectra.
An important aspect of survey disentangling is the sheer volume of data to be processed, which requires short optimisation times for any one object and consistent optimisation.We have parallelised the code using Python's multiprocessing library.On 72 CPUs, the computation of the parameters and spectra of 120 simulated systems (2 masses, 2 resolution regimes, 5 different light ratios, and 6 different velocity semi-amplitudes, see section 3.1 for details) took about 3 minutes.This includes running each system through the optimiser 6 times (once for each epoch that is removed from the data) for bootstrapping.As the systems are all independent of each other, the process scales without much difficulty to larger clusters with more cores.This bodes well for the application to even a million stellar systems.
The examples provided in this work have demonstrated that even in the regime of few epochs, moderate resolution, and signal-tonoise ratio, disentangling poses a very viable option for analyzing spectroscopic binaries.This, as well as the fast runtime of the algo- The other two panels display the disentangled rest-frame solutions for the primary and secondary, respectively, in magenta and cyan.The dashed, darker lines in these two panels indicate the model spectrum found by E22 to most closely fit the disentangled solutions found there.The bottom panel shows the residual between the reconstructed and observed composite spectra.We see generally good agreement between the solutions from this work and E22 models.rithm, allows disentangling to be performed on large surveys, such as LAMOST (Cui et al. 2012) and SDSS-V (York et al. 2000).The SDSS-V catalog will contain multi-epoch spectra of ∼ 200, 000 hot and massive stars, which are predicted (Sana et al. 2012;Moe & Di Stefano 2017) to have a high multiple fraction, making them great candidates for disentangling.Thus, the method developed here will allow us to weed out contaminants in our continued search for dark companions, and select ideal targets for higher-resolution, better-S/N and more-epoch follow up investigation.

A3 Matrix Structure: Details
The index of the first non-zero element in the first row of each block of M depends on the size of the shift.
Let    be the shift of the spectrum of component  at epoch  in units of the resolution element : We then make use of the   () operation, which returns the next lowest integer for any number ; for example,   (3.2) returns 3,   (−1.5) returns −2,   (0.7) returns 0, etc.For component  at epoch , then: or Λ  ,+1 in Equation A7) gives the index of the second non-zero element of the first row of each block of the matrix, or the shift of its associated (off-)diagonal, with values    .For a shift in the negative Λ direction that is smaller than 1 integer pixel, the block for component  and epoch  would look as follows: The task of disentangling the composite spectra from all epochs into two rest-frame spectra can then be expressed as a minimisation of the sum over epochs: This can be written, more explicitly, as the sum over epochs and pixels: (A12) Then using the interpolation described, as well as the notation in Equation A10, we get:

A4 Selection and scaling of optimisation parameters
For most optimisers it is advantageous to have normalised parameters.In this case, we are performing some rescaling by taking the common logarithm of , and dividing all epoch velocities (ì   ) by 100.This ensures that the optimiser can take similarly large steps in different directions with relative ease, as well as changing the mass ratio to a more sensible scale -a step from a  = 1 to  = 2 is relatively large, meaning a doubling of the secondary mass, while a step from  = 50 to  = 51 will cause very little difference in the overall solution.In log-space however, steps of similar size have a more consistent effect along the parameter range.
Even with these considerations finding these parameters is not trivial, considering the dimensionality of the parameter space is frequently high.(Simon & Sturm 1994), Hadrava (1995), Sablowski & Weber (2017) and others have suggested optimising over the orbital parameters instead, which can reduce the dimensionality of the problem to 7, which, in the case of many epochs, can be significantly less than fitting each epoch velocity individually.Here, however, we wish to demonstrate that even for a smaller number of epochs, spectral disentangling is viable, and thus optimise on the velocities and mass ratio.

A5 TIRAVEL
If we wish to forgo using templates to find initial guesses for the velocities of the primary, or do not have suitable templates available, we can employ the TIRAVEL algorithm (Zucker & Mazeh 2006) in order to find the relative inter-epoch velocities.This algorithm is most suitable to SB1s, or SB2s with comparatively faint secondaries.TIRAVEL achieves this by computing the cross-correlations between all epochs and then determining the velocitiy vector, with an offset, that solves for all the inter-epoch velocities by maximising the value of the CC at the inter-epoch velocity shifts.The offset arises from the fact that this method employs no templates, and thus no knowledge of what a rest-frame spectrum looks like is included.In this implementation, we set the observed RV of the first epoch to 0, and then calculate the other RVs relative to this.Once TIRAVEL has computed a candidate velocity vector, we create an estimation of the primary's spectrum by shifting all epoch spectra by their relative RVs and co-adding the spectra.This should, in theory, minimise the contribution from the secondary (in an SB2) and boost the signal of the primary -being more efficient the more epochs are available.
The gross offset of all velocities from their true values might seem problematic at first, but we remind ourselves that the systemic velocity is corrected for in a later step in our disentangling pipeline.This also accounts for any constant offset among all velocities that arises due to TIRAVEL, and thus this should cause no further issues.

B1 High Resolution Simulations
We repeat the simulations from section 3 for the same systems, but now with a higher resolution of  ≈ 20, 000.Here, the disentangler's difficulty with (close to) equal mass binaries becomes even more apparent, demonstrating that this is not a resolution issue.We also note that there are now fewer problems arising for small velocities (such as for the K=50 km/s systems), as the resolving power is now high enough to accurately determine even the smaller velocity shifts.Further, we see that the higher mass (3   ) primary still generally poses a bigger challenge than the lower mass (1.1   ) primary, again, due to the wider lines and fewer features overall.

B2 The "Giraffe"
We perform "blind" disentangling on the Giraffe system, both in the native resolution regime of  ≈ 60, 000, and on spectra artificially reduced to  ≈ 2, 000, as with the Unicorn in section 3.2.
In Figure B3 we see the performance of the optimiser on the parameters of the Giraffe in the higher resolution regime, with the wavelength window the same as for the Unicorn.Compared to the Unicorn, we note a markedly higher difficulty in recovering the light ratio correctly.This is likely due to the fact that the primary and secondary spectra of the Giraffe look more similar to each other than those of the Unicorn.Due to the rapid rotation of the Unicorn's secondary, the lines of that component are significantly broadened, leading to a somewhat non-typical spectrum that is quite distinct in shape from the spectrum of the primary.Here, however, the secondary is rotating much more slowly, thus this effect is minimal.Looking at Figure B5 lends credence to this idea, as we can see the two spectra are fairly similar in shape.We also see the effect of the underestimation of the light ratio: relative to the E22 models, the lines of the primary spectrum appear deeper (meaning a larger fraction of total light is attributed to the primary, and thus the light ratio is lower).
Figure B4 shows the performance of the optimiser on the data at the artificially reduced resolution, again on a larger wavelength window analogous to the Unicorn.We note the same difficulty with accurately assessing the light ratio, as well as some small issues with the velocities, and, resulting from this, the mass ratio, owing to the lower resolving power.The effects of the incorrect light ratio in particular become very apparent in Figure B6.Many lines of the primary appear deeper in the E22 models than in the spectrum found here, due to the too-low light ratio found by the optimiser.
c r o s sc o r r e l a t e c r o s s -c o r r e l a t e

Figure 3 .
Figure 3.An example of a disentangled spectrum, displaying one epoch.The top panel shows the observed spectrum at that epoch in black, and the sum of the two disentangled and velocity-shifted component spectra in grey.The second and third panel show the rest-frame spectra for the disentangled components, with the true spectra shown as dotted lines.The bottom panel displays the difference between the observed and the reconstructed spectrum (shown separately in the top panel).

Figure 4 .
Figure 4.The spectrum of the primary as obtained from disentangling with varying degrees of regularisation.The first panel shows the result without regularisation, the second panel with  = 5 and the last with  = 500.In each panel, the true spectrum is shown with a black dashed line.

Figure 5 .
Figure 5.An example of the results for the different parameters obtained by the optimiser.The correct value is plotted along the x-axis, with the value recovered by the optimiser along the y-axis.The dashed line indicates agreement between the two, and the further away from this line a value falls, the less accurately it was determined.The colour of each symbol also indicates how well the parameter was fit, with pink indicating a bad, and cyan a good agreement of the fit value with the truth relative to the other parameters.The symbol used for each point represents the parameter, with numbers 0-5 indicating the per-epoch RVs for each epoch (divided by 100),  the  10 of the mass ratio,  the systemic velocity (divided by 100) and  the  10of the light ratio.We selected these scaling of the different parameters for better visualisation.We see good recovery for all parameters, except two of the RVs.This is, however, not further surprising, as these RVs are fairly close to the systemic velocity, meaning there is only a small doppler shift of the two spectra relative to each other, making determination of the accurate velocities difficult at these low resolutions.

Figure 6 .
Figure6.A grid of the results of the autonomous disentangler on a range of different systems.The left panel shows the case of a low mass primary (1.1  ⊙ ), while the right panel displays the case of a higher mass primary (3  ⊙ ).Each panel consists of multiple cells, varying in the semi-amplitude of the radial velocities, , along the y-axis, and in light ratio , and consequently, the mass ratio  and effective temperature of the secondary,  eff,2 , along the x-axis.The colour of each cell indicates the figure-of-merit value the disentangler achieved for the recovery of the parameters, with white indicating a higher, and thus better, figure-of-merit, and black indicating a lower, worse FOM.Each cell also contains the accuracy and precision achieved for each individual parameter, analogous to Figure5.The x-axis here displays the true value of the parameter, and the y-axis the one recovered by the optimiser.If the optimiser determined the correct parameter, we expect it to lie along the dashed grey diagonal line.The colour of each parameter indicates how well it was recovered within its specific system, with pink indicating a comparatively large distance between truth and fitted parameter, and cyan indicating a small distance, and thus good agreement.Lastly, the size of each symbol is related to its variance in the bootstrapping process, with a large symbol indicating a large variance (and thus a low precision) and a small symbol indicating a small variance and higher precision.

Figure 7 .
Figure 7. Analogous to Figure6, with the light ratio now varying between 0.1 and 0.9 in steps of 0.2.

Figure 8 .
Figure8.A selection of individual results of the autonomous disentangler, here for a RV semi-amplitude of 200 km/s.The left column shows the low-mass primary case, and the right the higher-mass primary.The different rows show different light ratios  and thus mass ratios  and effective temperature of the secondary,  eff,2 .The top panel (white background) in each row shows the result of the disentangling for the primary, and the bottom (grey background) for the secondary.For each panel, the truth is plotted in black, with the disentangling result achieved by using the "correct" input parameters in green, and the actual result cyan the autonomous disentangler in magenta.The disentangler solutions are scaled by the light ratio (true value for cyan, value recovered by the optimiser for magenta) to match them to the true input.Thus, a "correctly shaped" spectrum that is scaled incorrectly indicates an incorrect light ratio, while a spectrum that is offset relative to the truth indicates that the systemic velocity was not found correctly.For the most extreme light ratios (top row), even with "correct" input parameters, the disentangler struggles to recover the correct spectrum of the secondary.The optimiser also has difficulties recovering the correct spectrum in the higher-mass case, owing to the wide lines found in hot stars.

Figure 9 .
Figure 9. Formal disentangling results for a simulated SB1 system (dark companion), analogous to figure 5.The single-epoch velocities of the luminous component are well-determined.The parameters   and  are ill-determined (because there is no second spectrum); the 1  range for one noise-realization is shown by the orange bands in the right panels.The value for  in such cases (implicitly an upper limit) is always found to be small.

Figure 11 .
Figure11.The parameters of the Unicorn as recovered by the optimiser (yaxis) compared to the ones found by E22 (x-axis).In the case of agreement, we expect the points for each parameter to lie on the grey, dashed, diagonal line.

Figure 13 .
Figure13.The results of the autonomous disentangling applied to the Unicorn for one epoch.The top panel shows the observed spectrum (black, solid line) as well as the reconstruction from the velocity shifted disentangled component spectra (grey, dashed line).The other two panels display the disentangled rest-frame solutions for the primary and secondary, respectively, in magenta and cyan.The dashed, darker lines in these two panels indicate the model spectrum found by E22 to most closely fit the disentangled solutions found there.The bottom panel shows the residual between the reconstructed and observed composite spectra.We see generally good agreement between the solutions from this work and E22 models.
For a shift in the negative Λ direction that is between 1 and 2 integer pixels, the block for component  and epoch  would look as followsFor a shift in the positive Λ direction that is smaller than 1 integer pixel, the block for component  and epoch  would look as follows:

Figure B3 .
Figure B3.The parameters of the giraffe as recovered by the optimiser (yaxis) compared to the ones found by E22 (x-axis).in the case of agreement, we expect the points for each parameter to lie on the grey, dashed, diagonal line.

Figure B5 .
Figure B5.The results of the autonomous disentangling applied to the Giraffe for one epoch.The top panel shows the observed spectrum (black, solid line) as well as the reconstruction from the velocity shifted disentangled component spectra (grey, dashed line).The other two panels display the disentangled rest-frame solutions for the primary and secondary, respectively, in magenta and cyan.The dashed, darker lines in these two panels indicate the model spectrum found by E22 to most closely fit the disentangled solutions found there.The bottom panel shows the residual between the reconstructed and observed composite spectra.We see generally good agreement between the solutions from this work and E22 models.
A10) Here,    gives distance in the positive Λ direction from    for each    -specifically, this    is equal to the    or    in Equation A7.