Spectrum of the secondary component and new orbital elements of the massive triple star δ Ori A

δ Orionis is the closest massive multiple stellar system and one of the brightest members of the Orion OB association. The primary (Aa1) is a unique evolved O star. In this work, we applied a two-step disentangling method to a series of spectra in the blue region (430 to 450nm), and we detected spectral lines of the secondary (Aa2). For the ﬁrst time, we were able to constrain the orbit of the tertiary (Ab) – to 55450d or 152yr – using variable γ velocities and new speckle interferometric measurements, which have been published in the Washington Double Star Catalogue. In addition, the Gaia DR3 parallax of the faint component (Ca+Cb) constrains the distance of the system to (381 ± 8)pc, which is just in the centre of the Orion OB1b association, at (382 ± 1)pc. Consequently, we found that the component masses according to the three-body model are 17 . 8, 8 . 5, and 8 . 7M (cid:12) , for Aa1, Aa2, and Ab, respectively, with the uncertainties of the order of 1M (cid:12) . We used new photometry from the BRITE satellites together with astrometry, radial velocities, eclipse timings, eclipse duration, spectral line proﬁles, and spectral energy distribution to reﬁne radiative properties. The components, classiﬁed as O9.5II + B2V + B0IV, have radii of 13 . 1, 4 . 1, and 12 . 0R (cid:12) , which means that δ Ori A is a pre-mass-transfer object. The frequency of 0 . 478cycles per day, known from the Fourier analysis of the residual light curve and X-ray observations, was identiﬁed as the rotation frequency of the tertiary. δ Ori could be related to other bright stars in Orion, in particular, ζ Ori, which has a similar architecture, or ε Ori, which is a single supergiant, and possibly a post-mass-transfer object.


Introduction
The bright star δ Ori (HR 1852, HD 36486, HIP 25930, ADS 4134) is a multiple stellar system consisting of six Tables B.1 and B.2 are also available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https:// cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/672/A31 Based on spectroscopic CCD observations with a coudé spectrograph attached to the 2m reflector of the Astronomical Institute AS ČR at Ondřejov, archival Haute-Provence and ESO La Silla spectra, groundbased UBV photometry from Hvar, and data collected by the BRITE Constellation satellite mission, designed, built, launched, operated, and supported by the Austrian Research Promotion Agency (FFG), the University of Vienna, the Technical University of Graz, the University of Innsbruck, the Canadian Space Agency (CSA), the University of Toronto Institute for Aerospace Studies (UTIAS), the Foundation for Polish Science & Technology (FNiTP MNiSW), and National Science Centre (NCN).
components: Aa1, Aa2, Ab, B, Ca, and Cb, more specifically, the eclipsing binary Aa1+Aa2, the interferometric binary (Aa1+Aa2)+Ab, the faint visual companion B, and the spectroscopic binary Ca+Cb (see Fig. 1).Their properties can be summarised as follows: -Aa1+Aa2 (V Aa1 = 2.55 mag, V Aa2 5.5 mag) 1 is a detached eclipsing binary with a negligible mass transfer, the orbital period P 1 = 5.732436 d (Mayer et al. 2010), a slightly eccentric orbit (0.08), and apsidal motion (1.45 • yr −1 ) (Pablo et al. 2015).-Ab (V Ab = 3.7 mag) is a nearby companion, which forms an interferometric pair with Aa1+Aa2.It was discovered by Heintz (1980), confirmed by speckle interferometry (Mason et al. 1999) and by Hipparcos astrometry (Perryman & ESA 1997) of the (Aa1+Aa2)+Ab system.The corresponding orbital period P 2 must be of the order of tens of thousands of days.-B (V B 14 mag) is a very faint distant companion 2 that is probably not associated with the system.Assuming that the component is a main-sequence star, its absolute magnitude of 6.7 mag corresponds to the spectral type K.
Many researchers have studied the system since the end of the 19th century.For a detailed summary of the early investigation of δ Ori, we refer readers to our earlier study of the system (Mayer et al. 2010).As far as studies of the 21st century are concerned, Harvin et al. (2002) carried out a tomographic separation of the ultraviolet and optical spectra into two systems of spectral lines, interpreted them as the lines of the primary and secondary of the eclipsing subsystem, and concluded that the components have unexpectedly low masses (m 1 = 11.2M and m 2 = 5.6 M ).However, Mayer et al. (2010), showed that the optical spectra are dominated by the spectral lines of the O9.5 II primary (Aa1; Walborn 1972) and the similarly hot tertiary (Ab), and that the system has normal masses for O and early-B stars (Harmanec 1988).The previous solution of the light curves (LCs) led Mayer et al. (2010) to the conclusion that the faint secondary (Aa2) contributes only a few percent to the total flux.Although they carried out disentangling of the spectra, they were unable to find its spectral lines convincingly, and could only rely on an indirect estimate of the mass ratio m 2 /m 1 .
Five in-depth studies of δ Ori were published in 2015 (the first four are a series): Corcoran et al. (2015) presented an overview of deep Chandra HETGS X-ray observations that covered nearly the entire binary (Aa1+Aa2) orbit.The observed X-ray emission was dominated by wind shocks from the primary (Aa1).Nichols et al. (2015) discussed the time-resolved and phase-resolved variability seen in the Chandra spectra.For the first time, they found phase-dependent variability in the Xray emission line widths.They identified two periods in the total X-ray flux: 4.76±0.30and 2.04±0.05days.Pablo et al. (2015) carried out a detailed analysis of space-based photometry from Microvariability and Oscillations of STars (MOST) and simultaneously secured ground-based spectroscopy in the residuals of the orbital LC, with periods ranging from 0.7 to 29 days.Shenar et al. (2015) carried out a multi-wavelength non-local thermodynamic equilibrium (NLTE) analysis of spectra.The determined parameters led to a O9.5 II, B1 V, and B0 IV spectral classification for Aa1, Aa2, and Ab, respectively, with evolved primary (Aa1) and tertiary (Ab) components.They also found wind-driven mass loss by the Aa1 component at 4×10 −7 M yr −1 .Richardson et al. (2015) used cross-correlation of the ultraviolet spectra from HST to obtain stellar parameter estimates for the primary, secondary, and the tertiary that was angularly resolved in the observations. 2α J2000 = 5 h 31 m 58.745 s and δ J2000 = −00 • 18 18.65 . 3α J2000 = 5 h 32 m 00.406 s and δ J2000 = −00 • 17 04.38 .In this work, we continue our earlier analysis (Harmanec et al. 2013), which was devoted to the detection of very weak He i 6678 Å lines of the secondary in the red spectral region.Hereinafter, we focus on the blue spectral region.This study was also motivated by the tentative evidence of the secondary reported by Richardson et al. (2015), namely in the ultraviolet region, observed by the Hubble Space Telescope (Space Telescope Imaging Spectrograph).
However, a robust detection of the secondary (Aa2) spectrum is still lacking.Now, we have a larger set of spectra in the blue part of the optical spectrum and procedures to successfully detect the secondary's spectrum.Moreover, new Gaia DR3 parallax measurements have been published.This provides the possibility to estimate the distance of bright stars, saturated in the Gaia images, from the measured distances of their faint companions.We also have new high-resolution astrometric measurements at our disposal, which enables us to constrain the long-period orbit of (Aa1+Aa2)+Ab.

Observational data
In this section, only the spectroscopic and photometric data sets are described as these data sets are new and fundamental to our analysis.Details of other data sets (astrometry, spectral energy distribution SED, speckle interferometry, etc.) are described in the following sections (Sects. 4,8,and 9).

Spectroscopy
We used digital spectra covering the blue spectral region secured at the coudé focus of the Ondřejov 2 m reflector (Škoda et al. 2002).We supplemented these data sets with spectra from the public archives of the ELODIE echelle spectrograph (Moultaka et al. 2004) at the Haute Provence Observatory, and the FEROS echelle spectrograph (Kaufer et al. 1999) at the ESO La Silla Observatory.The journal of the observations is presented in Table 1 (see Table B.1 for more details).The coverage of orbital phase ϕ 1 is illustrated in Fig. 2. The short period P 1 of 5.732436 d is well covered.The mean signal-to-noise ratio (S/N) is 208.5 (S/N values of individual spectra are given in Tables B.1 and B.2), which was sufficient for spectral disentangling.We normalised the spectra using polynomials of degree at least 4, with the program reSPEFO24 (written by Adam Harmanec).

Photometry
We used space-based photometric data from instruments on board the BRITE (BRIght Target Explorer; Pablo et al. 2016) and the MOST (Carroll et al. 1998) satellites and ground-based photometric data obtained at the Hvar Observatory with the 0.65m telescope.The time coverage is illustrated in Fig. 3.We did not use the saturated photometry from the Transiting Exoplanet Survey Satellite (TESS).
Each BRITE nanosatellite hosts a telescope, which has a 3 cm aperture.The BTr, BHr, and UBr satellites are equipped with a red filter (with effective wavelength 620 nm); BAb and BLb have a blue filter (420 nm).We have eliminated instrumental effects from the raw BRITE data by removing outliers and worst orbits, and by decorrelations.For more information on BRITE data processing, see Pigulski (2018).
The MOST passband covers the visible range of the spectrum (350-750 nm).The satellite performs high-precision optical photometry of single bright stars.It is equipped with a Maksutov telescope with an aperture of 15 cm and a custom broadband filter.It can point with an error of less than 1 arcsec.Other information can be found in Table 2.
The δ Ori LC from MOST continuously covers 3 weeks of observation.During calibration, we numerically shifted the measured magnitude to the V magnitude from the differential photometry at the Hvar Observatory.Then, we constructed normal points by centring the errors on the satellite orbital periods from Table 2, omitting the points with larger than the average uncertainty (0.5 mmag).
The Cassegrain 0.65m f /11 telescope at the Hvar observatory is equipped with a photoelectric detector (Božić 1998).This telescope was constructed at the Ondřejov Observatory of the Czechoslovak Academy of Sciences and brought to the Hvar Observatory at the beginning of 1972.A monitoring programme of bright variable stars has continued until today.The Hvar all-sky photometry provides accurate UBVR magnitudes in the Johnson system.For δ Ori A, we used UBV differential magnitudes obtained between October 2006 and October 2008 and UBVR between January 2019 and March 2021.(2006).
Hereinafter, we assume that the components (Aa1+Aa2)+Ab as well as Ca+Cb are located at the same distance.Statistically, they are located close to each other.The number of stars brighter than Ca+Cb (6.62 mag) is limited, there is only 15 of them within 7200 .Given the separation of 52 , the probability that stars are physically unrelated is low, p < 10 −3 .To the contrary, δ Ori B, which is also a formal member of the multiple visual system ADS 4134, is located at a substantially smaller distance (by almost 100 pc).It is therefore not physically related to δ Ori A. Either way, it is too faint (14 mag) to affect our results.
The Orion OB1 stellar association is usually divided into four subgroups, OB1a, OB1b, OB1c, and OB1d (Brown et al. 1994).The system δ Ori belongs to OB1b.We used the distances of 131 members from the Gaia DR3 catalogue and estimated the median distance to be (382±1) pc, using a cumulative distribution function that is sensitive to the local number density of stars (see Fig. 4).We obtained the same distance as the distance of δ Ori Ca+Cb, within the respective intervals.We consider this to be an independent estimate for the δ Ori A system since massive stars are often located in the centre of the given association.
Other bright stars in the Orion belt are also located at very similar distances (Table 3).For instance, the faint components of ζ Ori C, σ Ori C, D, and E all have precise parallaxes.Moreover, the single star ε Ori has a similar spectroscopic distance modulus.Again, this is an independent confirmation for the δ Ori system.
For comparison, the dispersion of distance in the radial direction (1σ) of the OB1b subgroup is only 15 pc, as seen in Fig. 4, while the angular dispersion (1σ) is about 0.5 • , which corresponds to 3 pc, at the distance of 382 pc.In other words, 1 A31, page 3 of 22 Table 3. Information about bright stars and their companions in Orion.(grey) associations (Brown et al. 1994).Parallax data were taken from the Gaia DR3 catalogue (Gaia Collaboration 2023).The median values of distance for OB1b and OB1a are (382±1) pc (brown) and (345±1 pc) (black), respectively.The distance of δ Ori C (red), (381±8) pc, is very close to the median distance of the OB1b distribution.

HD Name
corresponds to 0.11 pc, and 1 to 0.0018 pc; this is a range of separations for the faint components discussed above.The age of the OB1b association is estimated between 4 and 5 Myr (Maucó et al. 2018)  5 .The OB1a subgroup (northwest) is older and at a smaller distance (by approximately 37 pc), while the OB1c and OB1d subgroups (north-east), including the Trapezium, are younger and at larger distances.
5 Some of the outliers seen in Fig. 4 might actually be former members of the OB1b association.If they were ejected at the typical speed of 10 km s −1 , they may travel 50 pc or 7.5 • in the radial or tangential directions.The same is true for δ Ori B.
Another data set incorporated into the model included values of the systemic velocities γ 1 of δ Ori (Harvey et al. 1987;Harvin et al. 2002), which vary between approximately 12 and 23 km s −1 (see Table 5).This should correspond to the radial velocity of the (Aa1+Aa2) component.We did not take into account data points with possible systematic errors in γ 1 , that is, blending with Ab (1910Ab ( , 1948)), low amplitude of RV curve K 1 (1951,1969,1981), and different γ 1 for Aa2 (1987,1997).In some cases, also the RV of Ab was measured.
In total, we had N = 88 data points and M = 8 free parameters, which means N − M = 80 degrees of freedom.The model resulted in the best-fit with χ 2 = 95, with contributions χ 2 SKY = 60 for astrometry and χ 2 RV = 35 for RVs.Although χ 2 > N − M, the fit is still acceptable.The RV amplitude is in agreement, as well as directly measured RV values of Ab, which is lower than γ 1 .
The resulting parameters and parameters that were fixed are shown in Table 4.We fixed the mass of the (Aa1+Aa2) components based on the phoebe2 model (Sect.7) and the distance d based on the parallax (Sect.3).The orbit is illustrated in Fig. 5.The fit of RVs is shown in Fig. 6.
We estimated the uncertainties of parameters using χ 2 mapping and verified the MCMC method.According to the χ 2 statistics, a 1σ level corresponds to χ 2 101.Since the distance was fixed, the uncertainties are relatively small (1% for the period, 10% for the mass).We determined the mass of the Ab component to be 11.0 M .Therefore, the total mass of the (Aa1+Aa2)+Ab system is around 37.5 M .
Mirror solution.We are aware of the existence of a mirror solution, with the opposite sign of inclination i 2 .It exhibits higher total mass (up to 52 M ), higher eccentricity (0.95), shorter period (40 000 d), closer periastron passage.The RV curve of (Aa1+Aa2) component is also opposite, with a 'spike' due to the eccentricity.According to our throughout testing, it always has a worse best-fit χ 2 , especially the χ 2 rv contribution.Moreover, in the mirror solution, the RVs of Ab are larger than γ 1 , which is incorrect.A more complex model is needed to test other constraints (see Sect. 9).

Spectral disentangling of residuals
After obtaining the reliable value of the long-period (P 2 ) of δ Ori Aa and Ab, we searched for the secondary's lines in the spectra.Our experience with the disentangling technique is that the result is often sensitive to the choice of initial values of the parameters.5 were also used for the fit, extending the time span to 44 000 days.Bottom: detail of the observed arc.The most precise astrometric measurements from 2013 and 2019 constrain the orbital period P 2 .
This is understandable since the χ 2 sum based on all data points of all spectra is a complicated function of the orbital elements, and it is easy to end up in a local minimum.Moreover, the rotationally broadened spectral lines of the primary (Aa1) and tertiary (Ab) blend with each other at all orbital phases and altogether dominate the spectrum.Consequently, the contribution of the faint secondary (Aa2) to the χ 2 sum is almost comparable to the noise.The mass ratio q 1 of the Aa2 and Aa1 components is therefore poorly constrained.Nevertheless, the lines of the secondary can be detected in the residuals by a procedure called a two-step disentangling.
To disentangle the spectra, we used the KOREL program developed by Hadrava (1995Hadrava ( , 1997Hadrava ( , 2004Hadrava ( , 2005)).Rebinning of the spectra to a linear scale in RV, needed as input for KOREL, was carried out using the HEC35D program written by P.H. 8 The relative fluxes for the new wavelength points were derived using the INTEP program (Hill 1982), which is a modification of the Hermite interpolation formula.It is possible to choose both boundaries of the desired spectral region, and the program smoothly interpolates the rebinned spectra with the continuum values of 1.0 at both edges.1899 1913 1927 1941 1954 1968 1982 1995 2009 RV [km/s] JD − 2400000 synthetic γ 1 synthetic RV 3 observed γ 1 observed RV 3 omitted residuals γ 1 residuals RV 3 Fig. 6.Synthetic RVs of the Aa1+Aa2 (green) and Ab (black) components, relative to the barycentre of the (Aa1+Aa2)+Ab system.We used a simplified two-body model and fitted data from Table 5, plotted with blue symbols.The residuals are plotted in red; the value of χ 2 = 95.The last most precise point is from this work.
Notes.It is variable due to the presence of the third (Ab) component.
In some cases, also the RV of Ab was measured.If the reference is not provided, the value is taken from the list of Harvey et al. (1987) or Mayer et al. (2010), where more information about RV observations is provided (in their App.A). ( * ) Denotes the data that were not included in the fit due to systematic errors (see text).
To account for the variable quality of the individual spectra, we measured their S/N ratios in the line-free regions and assigned each spectrum a weight w according to the formula: where (S /N) mean denotes the root mean square of S/N ratio of all spectra.
Fitting with KOREL was performed with the following measurement equation: where j denotes the component; i, the spectrum; I, the normalised intensity; F I, its Fourier transform; y, the Fouriertransformed quantity x ≡ ln λ/λ 0 , related to the wavelength λ; s i j , the intensity factors (constant or variable); v i j , the radial velocity.
Two-step method.We used the period derived by Mayer et al. (2010), pericentre rate derived by Pablo et al. (2015), and parameters from Table 4 as the initial conditions.With the method of spectral disentangling, we needed to detect a line spectrum of the secondary in the blue spectral region 4275-4509 Å.
In the first step, we fitted the orbit of the close pair (Aa1+Aa2) and converged q 1 , e 1 , ω 1 , K 1 , T 0 , while P 1 , ω were fixed as well as the outer orbit (Sect.4).We set the same and constant intensity of lines of Aa1 and Ab (s 1 = 1, s 3 = 1) constant and assigned zero intensity to Aa2 (s 2 = 0).The result of the first step was the disentangled spectra of only the primary (Aa1) and tertiary (Ab), and the residuals for all individual spectra after disentangling (O−C).
In the second step, we added a value of 1.0 to the residuals and reran KOREL on this 'residual' data set.Now, the intensity factors of Aa1, Ab were zero (s 1 = 0, s 3 = 0) and the intensity factor of Aa2 was constant s 2 = 1.We fitted the spectrum of the Aa2 component by converging the mass ratio q and fixing T 0 , e 1 , ω 1 , K 1 , P 1 , and ω1 .We successfully detected the desired spectrum of the Aa2 component.The determined parameters are summarised in Table 6.This method gave higher e 1 , lower q 1 , precise K 1 and ω 1 , which were well constrained.All disentangled spectra in these two steps have a flat continuum, not wavy.To confirm the detection, we created a pseudo-χ 2 map (see Fig. 7).9 ; uncertainties as in Table 6.
Three-step method.A more precise orbital solution can be obtained by using KOREL in a sequence (three-step method).We started the process by fitting the primary and tertiary (Aa1+Ab) with variable intensities s 1 , s 3 .We fixed P 1 and ω1 of the close orbit and converged T 0 , e 1 , ω 1 , K 1 .The outer parameters of the orbit were fixed.
We continued by fitting Aa1, Ab with the constant s-factors and Aa2 with the variable one.Except for q 1 , all parameters were fixed.Finally, we found the solution for all three components with constant s-factors, free T 0 , e 1 , ω 1 , K 1 , q 1 , and fixed P 1 , ω1 .
The resulting s 1 (t) is variable with time and should correspond to the LC; however, the amplitude of the eclipses (without reflection) is too low (0.04 mag) to be seen.The resulting value of q 1 = 0.4517 is higher, compared to the two-step disentangling, while e 1 = 0.0761 is close to that found from the LC, and K 2 = 239.7 km s −1 .
In the three-step method, which we considered to be more reliable, we also computed the radial velocities of all three Notes.We prefer the solution from the three-step method (bold).The anomalistic period P anom,1 and the pericentre rate ω1 were fixed.Free parameters were the time of periastron passage T 0 , eccentricity e 1 , argument of periastron ω 1 , semi-amplitude of the primary K 1 , mass ratio q 1 , standard deviations of the intensity factors for the primary, secondary, and tertiary, σ s 1 , σ s 2 , σ s 3 , respectively.The dependent parameter is the semi-amplitude of the secondary K 2 .The models are quantified by pseudo-χ 2 in Fourier space.Fig. 7. Pseudo-χ 2 in Fourier space vs. mass ratio q 1 = m 2 /m 1 .This is related to the two-step disentangling method, to its second step, when the signal of the secondary component (Aa2) was sought for in the residuals.
components (see Tables B.1 and B.2).We estimated the uncertainties as a standard deviation weighted by S/N.

Atmospheric parameters of Aa1, Aa2, and Ab
We used the disentangled blue spectra to estimate the atmospheric parameters, namely, T eff , log g, v sin i, and the relative luminosities of the three components using the program PYTERPOL (Nemravová et al. 2016) 10 .The program uses the simplex minimisation technique to fit the synthetic spectra to the observed ones.As model spectra, OSTAR and BSTAR grids (Lanz & Hubený 2003, 2007) were used.
The results are summarised in Table 7, where the uncertainties were estimated from several independent trials.The fitted spectral line profiles of all components are shown in Fig. 8.Most of the lines are fitted reasonably well, except for He i 4471.The value of log g was determined primarily from the Hγ wings.The metallicity Z is not well constrained.
The sum of relative luminosities that were fitted independently (0.692 + 0.035 + 0.194 = 0.921) is close to 1, which is  (f) indicates the fixed parameter.For comparison with previous results, the values from Shenar et al. (2015) are shown in grey.In the case of Aa1 and Ab components, the results are usually in agreement within the uncertainties.More significant differences are for the Aa2 component; however, our values for the secondary are constrained by the disentangled spectra and mass ratio from KOREL.The uncertainties of the parameters are given in concise form in brackets.
an independent verification of the correctness of the KOREL disentangling.The effective temperatures agree with the spectral classifications of Aa1 (O9.5 II) and Ab (B0 IV), although the effective temperature of the Aa2 component is significantly lower (around 25 000 K) than that of the other components, corresponding to B1 V, according to Harmanec (1988) calibrations.The metallicities were fixed to the solar value since they are not well constrained by the blue spectra containing only one strong magnesium line.The values from Table 7 are the initial parameters for the phoebe2 model.

Orbit of eclipsing binary Aa1+Aa2
For the eclipsing binary Aa1+Aa2, we solved the inverse problem using PHOEBE2 (Conroy et al. 2020), obtaining a more precise model than with PHOEBE1 (Prša & Zwitter 2005) in our preliminary analysis (Oplištilová et al. 2020).The initial values of parameters for PHOEBE2 were inferred from the analysis performed with PHOEBE1.We had three photometric data sets available for analysis, SMEI, MOST, and BRITE (Oplištilová et al. 2020), but we preferred to use only BRITE data to have a homogeneous data set spanning nine seasons.MOST was used in Sect.9. We did not use SMEI data since they suffer from a contamination problem because stellar images in SMEI image have angular sizes of the order of 1 degree.
PHOEBE211 , a Python module, is software for modelling eclipsing binaries.To achieve the smallest possible discretisation error, the software uses a mesh of triangular elements.Each element of the mesh is assigned local properties (e.g., temperature, intensity), and the eclipse algorithm determines which elements are visible, which are partially visible, and which are not visible at all.The total flux is obtained by integrating over all visible elements.
We implemented a custom object-oriented Python wrapper to construct a model of the eclipsing binary and combine different data types.Each model was quantified by the χ 2  7).
value.First, we fitted the stellar parameters using the simplex method (Nelder & Mead 1965) or the supblex method (Rowan 1990).Second, we used the Markov chain Monte Carlo method (MCMC; Robert & Casella 2011;Tierney 1994), which was originally invented by Slanisław Ulam alongside with the atomic bomb.This method uses a sequence of random samples and provides a straightforward algorithm for the numerical estimation of parameters and their uncertainties.In other words, it describes the topology of the parameter space in the vicinity of the local/global minimum.The MCMC method was run using the API OpenMP (application program interface), which allows our code to run in parallel on multiple CPUs.
The MCMC method relies on Bayes' Theorem, which relates four probabilities as follows: where D denotes the vector of data; Θ M , the vector of parameters of our model; P( D), the probability of obtaining the data (normalisation); P(Θ M ), the prior, a priori knowledge of parameters (we used uniform, uninformative priors); P(D|Θ M ), the likelihood function, which is equivalent to the forward model or χ 2 ; and P(Θ M | D), the posterior distribution, which quantifies our belief in the parameters after combining our prior distribution with the current observations and normalising by the overall evidence.
The input data for the script are the RV curves of the primary and secondary, and the LCs in the blue and red filters.The synthetic fluxes were normalised by two free parameters S red and S blue satisfying: where c denotes the colour of the filter (blue or red), and i is the point number.
We set the algorithm parameters as follows: for the spatial discretisation, we used 1500 triangles covering the surface of the primary and 500 triangles for the small secondary surface.As a sampler, we used emcee (Foreman-Mackey et al. 2013) with 30 walkers and 2000 iterations.After some initial tests, we set the number of initial steps (burn-in) to 300.These are not taken into account as they are irrelevant and randomly distributed within the prior.The program ran on 30 CPUs.
In our modelling, we fixed the orbital sidereal period to 5.732436 d following Mayer et al. (2010), the pericentre rate ω = 1.45 • y −1 (Pablo et al. 2015), and in some models also the effective temperature of the primary T 1 = 31 000 K, and the third light, additional to the components Aa1 and Aa2, l 3 = 0.26685 calculated from Table 7 (Sect.6).
We used the following parameters in our model: atmosphere, black-body (approximation), limb darkening, linear, limb darkening coefficients, interpolated based on van Hamme (1993), gravity brightening, 1.0 (corresponding to the β coefficient for gravity darkening corrections), reflection and heating fraction, 1.0, distortion method, Roche, irradiation method, Wilson (1990), Wilson's original reflection effect scheme incorporates all irradiation effects, including reflection and redistribution, radial velocity method, flux-weighted (i.e.radial velocities are determined by the radial velocity of each element of visible surface, weighted by its intensity).Consequently, the RV curve includes the Rossiter-McLaughlin effect.We did not take into account either the effects of light travel time or gravitational redshift.This setting of the phoebe2 model is used for all models in Sects.7.1 and 7.2.

Model for season 2016
We had the BRITE LCs from 9 seasons at our disposal (Table 8).First, we selected well-covered season 2016 and fitted several A31, page 8 of 22 models with some parameters free or fixed, namely the effective temperature of the primary and the third light (for both blue and red filters).The results are presented in Table 9.We prefer the model with the fixed effective temperature of the primary, which also has the lowest value of χ 2 .The data, the model, and the residuals are shown in Fig. 9.Then, we used the MCMC method to estimate the uncertainties of the parameters.Figures 10 and 11 show the corner plot and the paths of walkers.In particular, masses M i and radii R i show strong positive correlations.In contrast, the inclination i 1 and R i show negative correlations due to geometrical reasons.In binaries, the sum of masses is inversely proportional to the third power of sin i; thus, i 1 and m i show negative correlations.The value of the systemic velocity γ is a little problematic as the value of 21.96 km s −1 was assumed and subtracted, then our model drifted to about −2.5 km s −1 , so that the resulting value is 18.5 km s −1 .The detached binary system Aa1+Aa2 is shown in Fig. 12.In addition, we derived several parameters from the nominal phoebe2 model (χ 2 = 604); see Table 10.We estimated the synthetic apparent brightness of δ Ori A as follows.The passband flux in Johnson V at the observer location is (in [W m −2 ]): where ∆ eff [m −1 ] stands for the effective wavelength range; ω = 1 m 2 /d 2 [sr −1 ], for the solid angle; d, for the distance of the system; k , for the summation over the triangular elements (grid); I λ [W m −3 sr −1 ], the monochromatic intensity on the stellar surface; S [m 2 ], the surface area of the element; µ = cos θ, where θ is the angle between the normal and the line of sight; η, the visibility in the range from 0 (visible element) to 1 (hidden or eclipsed element).We assumed the monochromatic calibration flux of Bessell et al. (2000): and the Johnson V passband flux: Notes.We fixed or released effective temperature T 1 of the primary and the third light l 3 .We prefer the model with fixed T 1 and l 3 (bold).The uncertainty σ is estimated the same for all models.The following numbers of data points were used: N total = 321, N LCB = N LCR = 100, N RV1 = 71, and N RV2 = 50.T 0 denotes the time of periastron passage ( ( * ) means HJD−2457733); T 1 and T 2 , the effective temperatures of the primary and secondary, respectively; R 1 and R 2 , the radii of the primary and secondary, respectively; i, inclination; S B and S R , the coefficients adjusting normalisation of the flux defined in Eq. ( 4); m 1 and m 2 , the masses of the primary and secondary, respectively; e 1 , eccentricity; ω 1 , the argument of periastron; γ, the systemic velocity; l 3R and l 3B , the third light in the blue and red filters, respectively; χ 2 sum , the total value of χ 2 ; χ 2 LCB , χ 2 LCR , χ 2 RV1 , and χ 2 RV2 , the contributions to χ 2 from the LCs in the red and blue filters and radial velocities of the primary and secondary; N total , the total number of data points; N LCB , N LCR , N RV1 , and N RV2 , the corresponding numbers of data points.Uncertainties were estimated from Figs. 11 and 13.
The apparent magnitude V 0 (without absorption) of the primary Aa1 and Aa2 is then: For the tertiary, we used the value of the third light: Comparing V 0 of the Aa1 + Aa2 + Ab system with Table 3, we get synthetic values 2.65 + 5.91 + 4.02 = 2.34 mag and observed values 2.42 + 5.4 + 3.70 = 2.08 mag.Thus, the total synthetic magnitude V 0 is about 0.26 mag fainter than observed.This result is acceptable, especially because the phoebe2 model was constrained only by the relative BRITE photometry (see also Sects.8 and 9).In other words, the result can be considered to be an independent confirmation of the distance.

Model for all observing seasons
Then, we made fits for all seasons.We kept the effective temperature of the primary T 1 and the third light l 3 fixed (as for the preferred model).The detailed results are presented in Table B.3.We took season 2016 as a reference.In Fig. 13, we show the deviations from the solution for season 2016.
In season 2018, the red filter measurements had greater uncertainties, however, it did not lead to any artefacts.We also omitted season 2019 since only measurements in the blue filter are available.We cannot confirm that variations of the parameters are intrinsic.Since the masses of the components must be constant, the mean values over all seasons should be preferred.The variations are most likely due to the oscillations.

Spectral energy distribution (SED)
The absolute flux is an additional observational constraint.In the case of δ Ori, the absorption is low because the star is not located behind the Orion molecular clouds.According to the reddening maps of Lallement et al. (2019) 12 , E(B − V) = 0.042 mag, with a substantial scatter of individual samples (due to the Orion

mag
Notes.log g denotes the logarithm of surface gravity, Φ V the passband flux at Earth, V 0 the corresponding apparent magnitude (without absorption, at a distance of 382 pc).
clouds).Therefore, the total absorption A V 3.1 E(B − V) = 0.130 mag, if extinction is not anomalous.There are not enough line-of-sight samples in the maps of Green et al. (2019)  13 .
We have removed clear outliers, multiple entries, and points without uncertainties.The IRAS 60 and 100 µm measurements show an excess, probably due to the far-infrared emission behind δ Ori; thus, they have been removed too.In the end, our photometric data set contained 31 data points (see Fig. 14).
For the Hvar photometry (Božić 1998), we performed a removal of eclipse phases (around 0.0, 0.45, 1.0), and computed average values at the maximum light.In this case, the absolute photometry is more reliable.The magnitudes transformed to the Johnson system are as follows: U = 0.940 mag, B = 1.977 mag, V = 2.221 mag, R = 2.334 mag, with uncertainties less than 0.010 mag.The comparison star used was HD 36591 (HR 1861): V = 5.341 mag, B − V = −0.190mag, U − B = −0.927mag, V − R = −0.050mag.In order to compare with the absolute flux, we used the calibrations from Bessell et al. (2000) (see also Fig. 14).

Three-body model with all observables
In order to account for additional observables in Xitau, we replaced the two-body model (Aa1+Aa2)+Ab with a three-body model Aa1+Aa2+Ab.Thus, the equations of motion were: where the first term is Newtonian gravitational interactions; the second, oblateness; and the third, relativistic effects.This model includes all relevant N-body perturbations (e.g., the radial velocities with respect to the common centre of mass, the light-time effect, precession of Ω 1 , 1 , Ω 2 , 2 , variation, evection; see also Appendix A), even though some of them are of minor importance for δ Ori.We included the oblateness of the bodies, parametrised with the Love number k F10 13 http://argonaut.skymaps.info/0.015 (Fabrycky 2010)14 , which results in the observed value of precession ω 1.45 • y −1 .Finally, we also included the parametrised post-Newtonian (PPN) approximation of relativistic effects (Standish & Williams 2006;Brož et al. 2022b) since the stars are both massive and close.The motion was integrated numerically using a Bulirsch-Stoer integrator, with a precision parameter ε = 10 −8 , and output every 0.2 d (plus exact times of observations).
Our model was constrained by astrometry (as in Sect.4), RVs of all components (Aa1, Aa2, Ab), eclipse timings, eclipse duration, LCs, synthetic spectra, and the SED.Individual contributions to the χ 2 metrics were multiplied by weights: where subscripts denote the data sets mentioned above, respectively.We used w SKY = w ECL = 10, due to the limited number of points, and w SYN = w SED = 0.1, due to the remaining systematics of rectification of spectra and absolute flux measurements.
Our model had 27 free parameters (see Table 11).The osculating elements are referenced to the epoch 2458773.188651(TDB), corresponding to the most precise speckle interferometry measurement.They are defined in the Jacobi coordinates, suitable for a system with hierarchical geometry.In this particular case, the distance d pc (Sect.3) was fixed.
We used the MOST LC (Pablo et al. 2015) to derive three times of the primary eclipse timings: 2456283.521, 2456289.277, and 2456294.994 (TDB, barycentric).Additional timings were obtained from the TESS (Ricker et al. 2014): 2458473.344, 2458479.080, 2458484.830 (TDB, barycentric).Due to large-amplitude oscillations, the uncertainty is degraded to 0.005 d.The primary eclipse duration is 0.667 d, with an uncertainty of 0.010 d, again due to the oscillations.We used a simplified eclipse algorithm for spherical stars.
At the same time, we computed the synthetic LC with the modified version of the Wilson-Devinney program (Wilson & Devinney 1971;Wilson 1979;Van Hamme & Wilson 2007;Wilson et al. 2010;Brož 2017); similarly as in Sect.7. In this case, however, the instantaneous true phase and distance were determined by the N-body integration.The third light is no longer an independent parameter; instead, it is determined by the third component (m 3 , T 3 , log g 3 ).This allowed us to constrain our model by eclipse depths.Other improvements included: a correction of computations for highly eccentric binaries, precise computations of the Roche potential from the volume-equivalent radius (Leahy & Leahy 2015), and more photometric filters (Prša & Zwitter 2005), including MOST.As the oscillations were not accounted for in the synthetic LC, uncertainties of 0.01 mag were assigned to all data points (see Fig. B.1).The observed spectra cover the blue region (430 to 450 nm).The synthetic spectra were interpolated by Pyterpol (Nemravová et al. 2016) from the BSTAR and OSTAR grids (Lanz & Hubený 2003, 2007).
We used the Planck (black body) approximation for the whole range of wavelengths, or absolute synthetic spectra for the limited range of the respective grids.The fit was performed with the simplex algorithm (see Fig. 15).We consider the best-fit model to be a compromise because it exhibits a tension between  9).Each sub-diagonal panel shows a 2D histogram, the isolines corresponding to the confidence intervals, and the correlations between parameters.
i) the synthetic spectra (in particular, log g 2 or R 2 ) and the duration of eclipses, ii) the minima timings, RVs, and oblateness (see also Fig. B.2).The best-fit parameters are summarised in Table 11, and the derived parameters in Table 12.Uncertainties were estimated by the χ 2 mapping and by the MCMC method.
Actually, for the Aa1, Aa2 components, they seem to be comparable to the phoebe2 model (Sect.7), but here we use a different and more extensive set of observations, in order to constrain all components at the same time.
The observed and synthetic SEDs are compared in Fig. 14.Even though the corresponding contribution χ 2 LC is larger than the number of data points N SED , we consider the fit to be acceptable as there are several multiple (but independent) measurements of the same band that are not consistent with each other.
At the same time, there is no systematic offset of the SED.In other words, our model provides independent confirmation of the parallax distance.
All blue spectra are shown in Fig. 16.There were remaining systematics between the observations and the model related to the rectification procedure, especially close to the He i 4387 line.While the Hγ was fitted without problems, the synthetic He i 4471 line was much shallower than the observed one.These problematic regions were removed from the fitting.These spectra constrain not only the RVs but also the relative luminosities L, log g, or radii R of all components.
Contributions of individual components are demonstrated in Fig. 17.Indeed, the secondary (Aa2) is faint (relative L 2 = 0.038).Unfortunately, its contribution is comparable to the A31, page 12 of 22 30 walkers, 1600 iterations, without burn-in Fig. 11.Parameter values versus iterations during the MCMC analysis (as in Fig. 10).We used 30 walkers (one walker corresponds to one colour) for the computation and 300 steps for the burn-in.For most parameters, the distribution of walkers is already stationary.The systemic velocity γ is stationary only after 1000 iterations.The value of 21.96 km s −1 should be added to γ.  systematics mentioned above.Consequently, some of the parameters are not very stable (in particular, log g 2 , v rot2 ).Nevertheless, our fitting in the direct space is independent and complementary to the fitting in Fourier space (Sect.5).Moreover, the secondary is constrained by other observables (e.g., eclipse duration, eclipse depth, RVs of the primary, the 3rd-body orbit, or the total mass m 1 + m 2 + m 3 ).
Mirror solution.Eventually, we explored the mirror solution (Sect.4).We fixed the total mass m 1 + m 2 + m 3 = 52.0M and performed a similar testing as in  B.3).In all models, the effective temperature T 1 of the primary (Aa1) and the third light l 3 were fixed.Season 2016 was taken as a reference.of the parameters can be found in Table 9.
χ 2 = 25468, which is worse than the nominal model.It exhibits a strong tension between the synthetic spectra and the SED.Especially the Hγ line profiles are fitted poorly (χ 2 syn = 69 441 vs. 44 795).This is directly related to the log g 3 value, which is very low (around 3.2) according to our model, as well as modelling of the disentangled spectrum of the tertiary component (Sect.6).Consequently, we exclude the mirror solution and prefer the nominal model presented above.

Pulsations
After obtaining a well-constrained model, we analysed the residual MOST LC (Fig. 18), in order to address the large-amplitude oscillations (BRITE data were not used because of instrumental issues).Our analysis is similar to that in Pablo et al. (2015), A31, page 13 of 22 synthetic observed residua w. reddening Planck 10 -3 10 -2 300 400 500 600 700 800  albeit it is different in several aspects: i) our model (from Sect.9) is constrained by all observables, ii) not only did we subtract the synthetic LC but we also removed the eclipse intervals (both primary and secondary), to suppress the binary signal.Consequently, the remaining frequencies should be preferentially related to rotation or pulsations, even though 'gaps' also create spurious signals.Notes.The best-fit model has non-reduced χ 2 = 23739.m tot denotes the total mass; q 1 = m 2 /m 1 , q 2 = m 3 /(m 1 + m 2 ) the respective mass ratios; P, osclulating period; e, eccentricity; i, inclination; Ω, longitude of node; , longitude of pericentre; λ, true longitude; T , effective temperature; g, gravitational acceleration; v rot , rotational velocity; C 20 , quadrupole moment; z 0 , magnitude zero point; γ, systemic velocity; and d pc , distance. (f) indicates the respective parameter was fixed.Orbital elements are osculating for the epoch T 0 = 2458773.188651(TDB).
We used the Period04 program (Lenz & Breger 2004) to compute the Fourier spectrum (Fig. 19), subtract the dominant term (prewhitening), recompute the spectrum again, and repeat these steps ten times.Our result is shown in Table 13.Notes.m denotes the mass; a, the semi-major axis; e, the eccentricity; R, the stellar radius; L, the relative luminosity (in V).
For reference, the minimum frequency (also spacing) is given by the time span of observations, f ∆ = 1/∆ = 0.045 c d   frequency of the primary (Aa1) by Nichols et al. (2015).Here, we interpret the second one as the rotation of the tertiary (Ab).This is confirmed by the rotational broadening.For parameters from Tables 11, 12, where we assumed an alignment.Interestingly, f rot1 = v rot1 /(2πR 1 sin i 1 ) = 0.174 c d −1 derived from rotational broadening is not equal to f 1 , which differs by f ∆ .It corresponds to the orbital frequency, which would indicate a synchronous binary.
The remaining frequencies ( f 3 , f 4 , f 10 ) are likely associated with pulsations, namely low-order modes = 0, 1, 2, or 3, typical for β Cep or SPB stars (Paxton et al. 2015).They can be present either on the primary (Aa1), or the tertiary (Ab), which contributes up to 40% of the light.

Conclusions
In this paper, we studied the triple star Aa1+Aa2+Ab in the multiple system δ Ori.The close eclipsing binary Aa1+Aa2 contains an O star, is non-interacting, and has a negligible mass transfer.Consequently, it represents a target suitable for defining the intrinsic parameters of evolved O stars.Our main results are as follows: 1.The distance of the system was estimated from the new Gaia DR3 parallax of the faint component δ Ori (Ca+Cb).2. The outer orbit (Aa1+Aa2)+Ab was constrained by new speckle interferometric measurements from the WDS (the period of approximately 152 years and eccentricity 0.58) and by γ velocities.3. The secondary (Aa2) spectrum in the blue spectral region was detected by the two-and three-step disentangling.4. The RV curve of the secondary was obtained by crosscorrelation with a disentangled template spectrum.5.The two-body model of the eclipsing binary was constructed in PHOEBE2.6.The three-body model in Xitau was constrained by all observables.Compared to previous studies, we obtained significantly lower masses than Pablo et al. (2015), in their low-mass model (23.81+ 8.54) M .In contrast to the study of Shenar et al. (2015), where models were calculated for two distances, 212 and 380 pc, we adopted the latter.Our results give lower radii of Aa1, Aa2; Shenar et al. (2015) have [(16.5±1)+ (6.5±2)] R .Nevertheless, the radius of Ab is in agreement with Shenar et al. (2015) who reported (10.4±2)R .We also obtained low log g 3 similarly as Shenar et al. (2015).
Hertzprung-Russel diagram.Given the spectral types of Aa1 + Aa2 + Ab, O9.5 II + B2 V + B0 IV (Pablo et al. 2015 and this work), the primary has evolved from the main sequence; however, it has not reached the overflow yet.The Hertzsprung-Russel diagram (Fig. 20) with the positions of Aa1, Aa2, Ab components indicates an interesting problem -the Ab component is very offset from a normal position.This offset is either related to its log g 3 value (3.2) or to its m 3 value (8.7 M ).However, it is not easy to modify these values, because they are well constrained by observations.To put Ab on the evolutionary track, either log g 3 3.7, or m 3 18 M .If all A31, page 16 of 22  (2015).Numbers next to the main sequence indicate the theoretical masses, small numbers the theoretical gravitational acceleration log g.According to the 3-body model, the masses are 17.1, 8.5, 8.7 M .For the tertiary (Ab), it is in disagreement with the evolutionary track, but in agreement with the value of log g 3 inferred from normalised spectra (Hγ).The instability region of pulsations (β Cep type, order = 0) is indicated as grey area.Other modes ( = 1, 2, 3) can be found in a very similar region (Paxton et al. 2015).Both Aa1, Ab components are located here, and they can exhibit photometric variability attributed to pulsations.
components were normal, the sum of masses should be about (24+18+10) M = 52 M .Interestingly, this is similar to the mirror solution, which was excluded (see the discussion in Sects.4 and 9).Consequently, we are left with an unusual stellar component.Actually, it is not unusual -see for example δ Ori Ca, or σ Ori E which are both helium-rich, with Hα emission (Table 3).Detailed stellar-evolution models with possible mass transfer between (some of) the components shall be computed in future work.Additionally, long-baseline optical or near-infrared interferometry may be able to measure precisely the angular diameters of the component stars (e.g., Shabun et al. 2008) and the separation of the inner orbit in the sky, giving direct constraints on the size of the orbit, helping to resolve any discrepancies in the masses measured between this study and other similar studies of δ Ori.
A comparison of δ Ori with other bright stars in the Orion belt (see Table 3) shows that σ Ori has a similar architecture (((Aa+Ab)+B)+C)+D+E, and even a very similar angular scale (Simón-Díaz et al. 2015).All of its components seem to be less evolved15 .On the other hand, ζ Ori exhibits an angular scale about 10 times larger and has the primary evolved in an O supergiant (Hummel et al. 2000).In this sense, ε Ori, which seems to be a single variable B supergiant (Puebla et al. 2016), may represent an even more evolved object.
Given the fact that all these stars (δ, ε, ζ, σ) are the most massive within the Orion OB1b association, they might have encountered and perturbed (destabilised) each other.Again, a possible convergence of their proper motions will be analysed in future work.because the respective radius (|r| = a 1 ) is the same as the reference radius (R = a 1 ).Given that J 2 = 0.109, and the ratio of a 1 /a 2 = 0.002, the precession rates should be of the order of 10 −6 deg y −1 .
However, for the massive triple system (Aa1+Aa2)+Ab, we used the theory of Breiter & Vokrouhlický (2015); for the longitudes of pericentre and node: where 2 , the angular momenta; cos J = cos i 1 cos i 2 +sin i 1 sin i 2 cos(Ω 1 −Ω 2 ), the mutual inclination; and γ = L 1 /(L 2 η 2 ).Assuming parameters from Table 11, ˙ 2 = 3.4 • 10 −6 deg y −1 , Ω2 = 8.5 • 10 −5 deg y −1 , which is negligible on the observational time span.All these effects were nevertheless included in our numerical N-body model (Sect.9).Every simplex was initialised with a different combination of the gravitational accelerations log g 1 , log g 3 , which were kept fixed to obtain a regular grid.All other parameters were free.The number of convergence steps was limited to 2000, consequently, 200000 models were computed in total.Axes correspond to the values of log g 1 , log g 3 , colours to χ 2 (see also tiny numbers).The colour scale was adjusted as follows: cyan the very best fit for a given data set, blue acceptable fits (<3.0 min χ 2 ), orange poor fits (≥3.0 min χ 2 ).The factor was 1.5 for the SKY, total; and 30.0 for the ECL data set.'Forbidden regions' can be seen, in particular, large log g 1 , log g 3 due to the SYN, SED, or large log g 1 due to the ECL, LC.The weighted very best fit is denoted by the red circle.

Appendix B: Supplementary figures and tables
A31, page 20 of 22  8) and all RVs were used to constrain the models.We assumed a fixed value of the temperature T 1 and the third light l 3 .The uncertainties are the same as in Notes.* −2457000 HJD.The explanation of variables is the same as in Table 9. f indicates the respective parameter was fixed.var denotes variable values for each season (they were between 1.004 and 1.010).For each season, the numbers of data points were: N total = 321, N LCB = N LCR = 100, N RV1 = 71, and N RV2 = 50.

Fig. 5 .
Fig.5.Orbit of 2-body model.Top: orbit of (Aa1+Aa2)+Ab components in the (u-v) plane (green), calculated using the two-body model and the simplex method.Observations are shown with blue symbols and uncertainty ellipses (orange).The residuals are plotted in red; the value of χ 2 = 95.The astrometry used for the fit is from the WDS.The radial velocities of Aa1, from Table5were also used for the fit, extending the time span to 44 000 days.Bottom: detail of the observed arc.The most precise astrometric measurements from 2013 and 2019 constrain the orbital period P 2 .

Fig. 8 .
Fig.8.Comparison of the disentangled spectra (blue) of the Aa1, Aa2, and Ab components with the best-fit synthetic spectra (orange) found by PYTERPOL.The range from 428 to 449 nm was used for disentangling.The flux is normalised to the local continuum.The small panels show the residuals (red).The relative luminosities of the Aa1 and Ab components significantly exceed that of the Aa2 component (see also Table7).

Fig. 9 .
Fig.9.Comparison of observations and the phoebe2 model of δ Ori with χ 2 = 604.The values of the effective temperature of the primary T 1 and the third light l 3 were fixed.The upper panel shows the phased LCs in the blue and red BRITE filters.The lower panel shows the RV curves for the primary Aa1 (green) and the secondary Aa2 (purple).The grey points correspond to our model, the red lines to the residuals, or contributions to χ 2 .

Fig. 10 .
Fig. 10.Corner plot (a full covariance matrix) for the δ Ori model, as derived by the MCMC analysis.The model is the same as in Fig. 9.Each diagonal panel shows a 1D histogram (posterior distribution) for one parameter (explained in Table9).Each sub-diagonal panel shows a 2D histogram, the isolines corresponding to the confidence intervals, and the correlations between parameters.

Fig. 12 .
Fig. 12. Mesh plot of δ Ori Aa1+Aa2 binary from the phoebe2 model (with χ 2 = 604).This is a (u, v) plane projection with scale in R , at phase ϕ = 0.75.The grey scale corresponds to the monochromatic intensity I λ [W m −3 sr −1 ] for the effective wavelength of the BRITE blue filter (420 nm).

FigFig. 13 .
Fig. 13.Relative differences of parameters derived for δ Ori for eight seasons (TableB.3).In all models, the effective temperature T 1 of the primary (Aa1) and the third light l 3 were fixed.Season 2016 was taken as a reference.of the parameters can be found in Table9.

Fig. 14 .
Fig. 14.Model from SED data.Top: comparison of the observed (blue) and synthetic (orange) SED of δ Ori.The residuals are plotted in red.The wavelength range is from 350 nm (ultraviolet) to 25 µm (farinfrared).Bottom: the same for the limited wavelength range of synthetic spectra.The Hvar differential UBVR photometry with removed eclipses is plotted in green.

Fig. 15 .
Fig. 15.Example of a convergence of the three-body model of δ Ori.The individual contributions to χ 2 correspond to astrometry (SKY), radial velocities (RV), eclipse timings (ETV), eclipse duration (ECL), synthetic spectra (SYN), and the spectral energy distribution (SED).The total χ 2 is summed with non-unit weights (see values in the main text).

Fig. 16 .
Fig. 16.Comparison of the observed (blue) and synthetic (orange) rectified spectra.The residuals are plotted in red.All components Aa1, Aa2, Ab of δ Ori contribute to the total flux.The wavelength range included the spectral lines: Hγ 4341, He i 4378, and numerous weaker lines.

Fig. 19 .
Fig. 19.Periodogram of the residual MOST LC with the 10 principal frequencies indicated on top.The first periodogram is shown; it is subsequently modified by subtraction (prewhitening).The broad peak between 10 to 20 c d −1 corresponds to the satellite frequency f most and its combinations with f 1 , etc.

Fig. 20 .
Fig.20.Hertzsprung-Russel diagram with the positions of the Aa1, Aa2, and Ab components of δ Ori and evolutionary tracks fromPaxton et al. (2015).Numbers next to the main sequence indicate the theoretical masses, small numbers the theoretical gravitational acceleration log g.According to the 3-body model, the masses are 17.1, 8.5, 8.7 M .For the tertiary (Ab), it is in disagreement with the evolutionary track, but in agreement with the value of log g 3 inferred from normalised spectra (Hγ).The instability region of pulsations (β Cep type, order = 0) is indicated as grey area.Other modes ( = 1, 2, 3) can be found in a very similar region(Paxton et al. 2015).Both Aa1, Ab components are located here, and they can exhibit photometric variability attributed to pulsations.

Fig. B. 1 .Fig
Fig. B.1.Comparison of the observed MOST (blue) and synthetic (orange) light curve.The residuals are plotted in red.Apart from the eclipses, the light curve contains large-amplitude oscillations (not included in our model); uncertainties 0.01 mag were thus assigned to all data points.In Tables B.1 and B.2, we present more details on the spectral data sets discussed in Sect.2.1.In Table B.3, we report parameters derived for the 8 seasons observed by the BRITE satellites.The LC from our three-body model (Sect.9) is shown in Fig. B.1.Individual contributions to χ 2 computed for an extensive set models of δ Ori is shown in Fig. B.2.

Table 1 .
Journal of digital spectra covering the blue spectral region.

Table 2 .
Information on satellites.

Table 6 .
Solution to the disentangling of 75 blue spectra in KOREL.

Table 7 .
Atmospheric parameters derived with PYTERPOL from the blue spectral region 4271-4513 Å, with the mean resolution of 0.0144 Å.Notes.T eff denotes the effective temperature; log g, logarithm of surface gravity; v sin i, projected rotational velocity; Z, metallicity; L R , relative luminosity; χ 2 R , the reduced value of χ 2 (divided by the degrees of freedom), and

Table 8 .
Time intervals of the BRITE LCs (red and blue filters) when δ Ori was observed.Notes.The whole time span was divided into 9 seasons.Time is given in HJD − 2 400 000.Season 2019 denoted by ( * ) was omitted since redfilter data were not available.In total, BRITE nanosatellites measured 499 656 raw data points.The numbers of points are shown for each season and filter.

Table 9 .
Results of fitting three phoebe2 models for δ Ori Aa1+Aa2.LC from season 2016 and all RVs were used to constrain the models.

Table 11 .
Free parameters of the three-body model of δ Ori.

Table 12 .
Derived parameters of the three-body model of δ Ori.

Table 13 .
Nichols et al. 2015)the residual MOST LC (from Fig.18).The synthetic LC was subtracted and eclipses were removed to suppress the binary signal.The first two frequencies are then identified as the rotation frequency of the primary (Aa1;Nichols et al. 2015)and the tertiary (Ab; this work).The uncertainty is determined by the time span, f ∆ = 0.045 c d −1 .For reference, the orbital frequency of the Aa1+Aa2 binary f orb = 0.174 c d −1 , and the orbital frequency of the satellite f most = 14.19 c d −1 .
Table B.2. Details on ELODIE and FEROS spectra in the blue region.ELODIE is at the upper part of the table, FEROS at the lower part.The RVs were determined during the spectral disentangling (three-step method) in KOREL.Table B.3. Results of eight phoebe2 models of δ Ori.LCs from eight individual seasons S. (Table