Issue 
A&A
Volume 672, April 2023



Article Number  A31  
Number of page(s)  22  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202245272  
Published online  27 March 2023 
Spectrum of the secondary component and new orbital elements of the massive triple star δ Ori A^{⋆,}^{⋆⋆}
^{1}
Charles University, Faculty of Mathematics and Physics, Astronomical Institute, V Holešovičkách 2, 180 00 Praha 8Trója, Czech Republic
email: betsimsim@seznam.cz
^{2}
Uniwersytet Wrocławski, Instytut Astronomiczny, Kopernika 11, 51622 Wrocław, Poland
^{3}
Hvar Observatory, Faculty of Geodesy, Zagreb University, Kačicéva 26, 10000 Zagreb, Croatia
^{4}
Czech Academy of Sciences, Astronomical Institute, 25165 Ondřejov, Czech Republic
^{5}
AAVSO, 185 Alewife Brook Pkwy, Cambridge, MA, USA
^{6}
Department of Physics, Mount Allison University, Sackville, NB E4L1E6, Canada
^{7}
Department of Physics and Space Science, Royal Military College of Canada, Kingston, Ontario K7K 7B4, Canada
^{8}
Institut für Astro und Teilchenphysik, Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria
^{9}
Department of Electronics, Electrical Engineering and Microelectronics, Silesian University of Technology, Akademicka 16, 44100 Gliwice, Poland
^{10}
University of Vienna, Institute for Astrophysics, Türkenschanzstraße 17, 1180 Vienna, Austria
Received:
23
October
2022
Accepted:
9
January
2023
δ 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 twostep disentangling method to a series of spectra in the blue region (430–450 nm), and we detected spectral lines of the secondary (Aa2). For the first time, we were able to constrain the orbit of the tertiary (Ab) – to 55 450 d or 152 yr – 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 threebody model are 17.8, 8.5, and 8.7 M_{⊙}, for Aa1, Aa2, and Ab, respectively, with the uncertainties of the order of 1 M_{⊙}. We used new photometry from the BRITE satellites together with astrometry, radial velocities, eclipse timings, eclipse duration, spectral line profiles, and spectral energy distribution to refine radiative properties. The components, classified as O9.5 II + B2 V + B0 IV, have radii of 13.1, 4.1, and 12.0 R_{⊙}, which means that δ Ori A is a premasstransfer object. The frequency of 0.478 cycles per day, known from the Fourier analysis of the residual light curve and Xray observations, was identified 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 postmasstransfer object.
Key words: binaries: close / stars: massive / stars: individual: δ Ori / binaries: eclipsing / stars: fundamental parameters / techniques: spectroscopic
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/vizbin/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 HauteProvence 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).
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The bright star δ Ori (HR 1852, HD 36486, HIP 25930, ADS 4134) is a multiple stellar system consisting of six 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 mainsequence star, its absolute magnitude of 6.7 mag corresponds to the spectral type K.

Ca+Cb (V_{Ca + Cb} = 6.85 mag) is another distant companion^{3} to (Aa1+Aa2)+Ab that is a spectroscopic, noneclipsing binary with a period of 29.96 d and of spectral types B3 V + A0 V (Leone et al. 2010).
Fig. 1. Scheme of the multiple system δ Ori (HD 36486, ADS 4134, Mintaka). Orbital periods were taken from ^{*}Leone et al. (2010), ^{**}Mayer et al. (2010), and ^{* * *}this paper. 
In the present paper, we focus on the triple subsystem δ Ori (Aa1+Aa2)+Ab, with V_{Aa1 + Aa2 + Ab} = 2.223 mag (from the differential photometry at the Hvar Observatory), α_{J2000} = 5^{h} 32^{m} 00.400^{s}, and δ_{J2000} = −00° 17′ 56.74″. Hereinafter, the parameters corresponding to the inner orbit Aa1+Aa2 and to the outer orbit (Aa1+Aa2)+Ab are denoted by indices 1 and 2, respectively. The parameters of the components Aa1, Aa2, and Ab are denoted by indices 1, 2, and 3, respectively.
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.2 M_{⊙} 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 earlyB 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 indepth studies of δ Ori were published in 2015 (the first four are a series): Corcoran et al. (2015) presented an overview of deep Chandra HETGS Xray observations that covered nearly the entire binary (Aa1+Aa2) orbit. The observed Xray emission was dominated by wind shocks from the primary (Aa1). Nichols et al. (2015) discussed the timeresolved and phaseresolved variability seen in the Chandra spectra. For the first time, they found phasedependent variability in the Xray emission line widths. They identified two periods in the total Xray flux: 4.76 ± 0.30 and 2.04 ± 0.05 days. Pablo et al. (2015) carried out a detailed analysis of spacebased photometry from Microvariability and Oscillations of STars (MOST) and simultaneously secured groundbased spectroscopy in the residuals of the orbital LC, with periods ranging from 0.7 to 29 days. Shenar et al. (2015) carried out a multiwavelength nonlocal 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 winddriven mass loss by the Aa1 component at 4 × 10^{−7} M_{⊙} yr^{−1}. Richardson et al. (2015) used crosscorrelation 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.
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 highresolution astrometric measurements at our disposal, which enables us to constrain the longperiod orbit of (Aa1+Aa2)+Ab.
2. 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).
2.1. 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 signaltonoise 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 reSPEFO2^{4} (written by Adam Harmanec).
Fig. 2. Coverage of the orbital period P_{1} by blue spectra. Phases are determined with respect to time T_{0} = HJD 2454002.8735 (time of periastron passage determined by KOREL) for the eclipsing binary. 
Journal of digital spectra covering the blue spectral region.
2.2. Photometry
We used spacebased photometric data from instruments on board the BRITE (BRIght Target Explorer; Pablo et al. 2016) and the MOST (Carroll et al. 1998) satellites and groundbased 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).
Fig. 3. Photometric data from MOST and BRITE displayed with respect to time. The BRITE data covers six consecutive seasons between 2013 and 2021. 
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 highprecision 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.
Information on satellites.
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 allsky 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.
3. Parallax and distance of δ Ori
In Gaia DR3 (Gaia Collaboration 2016, 2021, 2023), the parallaxes of the faint components of bright stars in the Orion OB1 association were measured (see Table 3). The parallax of δ Ori Ca+Cb, π = (2.6244 ± 0.0538) mas implies a distance d = (381 ± 8) pc and a distance modulus μ = (7.90 ± 0.04) mag. 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}.
Information about bright stars and their companions in Orion.
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.
Fig. 4. Sorted distances of objects from the Orion OB1b (blue) and OB1a (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. 
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′ 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 (northeast), including the Trapezium, are younger and at larger distances.
4. Visual orbit of (Aa1+Aa2)+Ab
To determine the parameters of the longperiod (P_{2}) orbit of the (Aa1+Aa2)+Ab system, we constructed a simplified twobody model in the Xitau program (Brož 2017; Brož et al. 2021, 2022a,b)^{6}. We used astrometric data from the Washington Double Star (WDS) Catalogue (Mason et al. 2001). If uncertainties were not available, we assumed the uncertainties of the separation ρ and the position angle measured from v, +DE direction as follows: σ_{θ} = 1.0°, σ_{ρ} = 0.01 mas, or σ_{θ} = 0.2°, σ_{ρ} = 0.005 mas for measurements before and after 2013, respectively. After removing 8 outliers due to poor resolution, incorrect plate scale, or calibration (from 1879.12^{7}, 1978.10, 1979.06, 1980.02, 1981.01, 1985.74, 1995.05, 1999.78), we used N = 74 data points (both ρ, θ).
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 (1910, 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 bestfit with χ^{2} = 95, with contributions for astrometry and 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.
Fig. 5. Orbit of 2body model. Top: orbit of (Aa1+Aa2)+Ab components in the (u–v) plane (green), calculated using the twobody 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 Table 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}. 
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 twobody 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. 
Parameters of the orbit (Aa1+Aa2)+Ab, for the model with unreduced χ^{2} = 95.
Observed values of γ_{1} velocity of the Aa1+Aa2 components.
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_{⊙}.
4.1. 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 bestfit χ^{2}, especially the 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).
5. Spectral disentangling of residuals
After obtaining the reliable value of the longperiod (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. 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 twostep disentangling.
To disentangle the spectra, we used the KOREL program developed by Hadrava (1995, 1997, 2004, 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.
To account for the variable quality of the individual spectra, we measured their S/N ratios in the linefree 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; ℱI, its Fourier transform; y, the Fouriertransformed quantity x ≡ lnλ/λ_{0}, related to the wavelength λ; s_{ij}, the intensity factors (constant or variable); v_{ij}, the radial velocity.
5.1. Twostep 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 . 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.
Fig. 7. Pseudoχ^{2} in Fourier space vs. mass ratio q_{1} = m_{2}/m_{1}. This is related to the twostep disentangling method, to its second step, when the signal of the secondary component (Aa2) was sought for in the residuals. 
Solution to the disentangling of 75 blue spectra in KOREL.
5.2. Threestep method.
A more precise orbital solution can be obtained by using KOREL in a sequence (threestep 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 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 sfactors 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 sfactors, free T_{0}, e_{1}, ω_{1}, K_{1}, q_{1}, and fixed P_{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 twostep disentangling, while e_{1} = 0.0761 is close to that found from the LC, and K_{2} = 239.7 km s^{−1}.
In the threestep method, which we considered to be more reliable, we also computed the radial velocities of all three components (see Tables B.1 and B.2). We estimated the uncertainties as a standard deviation weighted by S/N.
6. 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.
Fig. 8. Comparison of the disentangled spectra (blue) of the Aa1, Aa2, and Ab components with the bestfit 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 Table 7). 
Atmospheric parameters derived with PYTERPOL from the blue spectral region 4271–4513 Å, with the mean resolution of 0.0144 Å.
The sum of relative luminosities that were fitted independently (0.692 + 0.035 + 0.194 = 0.921) is close to 1, which is 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.
7. 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.
PHOEBE2^{11}, 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 objectoriented Python wrapper to construct a model of the eclipsing binary and combine different data types. Each model was quantified by the χ^{2} 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 (ForemanMackey et al. 2013) with 30 walkers and 2000 iterations. After some initial tests, we set the number of initial steps (burnin) 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 (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, blackbody (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, fluxweighted (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 RossiterMcLaughlin 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.
7.1. Model for season 2016
We had the BRITE LCs from 9 seasons at our disposal (Table 8). First, we selected wellcovered season 2016 and fitted several 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.
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 (forestgreen greenblack) and the secondary Aa2 (pypurple purpleblack). The grey points correspond to our model, the red lines to the residuals, or contributions to χ^{2}. 
Time intervals of the BRITE LCs (red and blue filters) when δ Ori was observed.
Results of fitting three phoebe2 models for δ Ori Aa1+Aa2. LC from season 2016 and all RVs were used to constrain the models.
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}.
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 Table 9). Each subdiagonal panel shows a 2D histogram, the isolines corresponding to the confidence intervals, and the correlations between parameters. 
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 burnin. 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 γ. 
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}]):
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). 
Derived parameters corresponding to the δ Ori Aa1+Aa2 model with χ^{2} = 604.
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:
where f(λ) denotes the filter transmission; λ_{eff} = 545 nm, the effective wavelength; and Δ_{eff} = 85 nm, the effective range.
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.
7.2. 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.
Fig. 13. Relative differences of parameters derived for δ Ori for eight seasons (Table 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. Explanations of the parameters can be found in Table 9. 
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.
8. 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 clouds). Therefore, the total absorption A_{V} ≐ 3.1 E(B − V) = 0.130 mag, if extinction is not anomalous. There are not enough lineofsight samples in the maps of Green et al. (2019)^{13}.
From the photometric catalogues in http://vizier.ustrasbg.fr/vizier/sed/doc/ (Allen et al. 2014), we took the standard Johnson system photometry (Ducati 2002) and measurements from HIPPARCOS (Anderson & Francis 2012), Gaia DR3 (Gaia Collaboration 2020), 2MASS (Cutri et al. 2003), WISE (Cutri et al. 2012), MSX (Egan et al. 2003), Akari (Ishihara et al. 2010), and IRAS (Neugebauer et al. 1984). The data covered the spectral range from 0.35 to 100 μm.
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 farinfrared emission behind δ Ori; thus, they have been removed too. In the end, our photometric data set contained 31 data points (see 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 redredblack. 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. 
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.190 mag, U − B = −0.927 mag, V − R = −0.050 mag. In order to compare with the absolute flux, we used the calibrations from Bessell et al. (2000) (see also Fig. 14).
9. Threebody model with all observables
In order to account for additional observables in Xitau, we replaced the twobody model (Aa1+Aa2)+Ab with a threebody 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 Nbody perturbations (e.g., the radial velocities with respect to the common centre of mass, the lighttime 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} ≃ 0.015 (Fabrycky 2010)^{14}, which results in the observed value of precession . Finally, we also included the parametrised postNewtonian (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 largeamplitude 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 WilsonDevinney 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 Nbody 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 volumeequivalent 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 bestfit model to be a compromise because it exhibits a tension between 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 bestfit 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.
Fig. 15. Example of a convergence of the threebody 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 nonunit weights (see values in the main text). 
Free parameters of the threebody model of δ Ori.
Derived parameters of the threebody model of δ Ori.
The observed and synthetic SEDs are compared in Fig. 14. Even though the corresponding contribution 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.
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. 
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 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 3rdbody orbit, or the total mass m_{1} + m_{2} + m_{3}).
Fig. 17. Same as Fig. 16, with the distinguished contributions of individual components Aa1, Aa2, Ab of δ Ori. The wavelength range is limited to 430–440 nm. Only the sixth spectrum, 2455836.582700 is shown. 
9.1. Mirror solution.
Eventually, we explored the mirror solution (Sect. 4). We fixed the total mass m_{1} + m_{2} + m_{3} = 52.0 M_{⊙} and performed a similar testing as in Fig. B.2. The overall bestfit has χ^{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 ( = 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.
10. Pulsations
After obtaining a wellconstrained model, we analysed the residual MOST LC (Fig. 18), in order to address the largeamplitude oscillations (BRITE data were not used because of instrumental issues). Our analysis is similar to that in Pablo et al. (2015), 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.
Fig. 18. Residuals of the MOST LC (red), and the Fourier series up to the order of ℓ = 10 (black; Table 13). Individual orders as they are summed up are shown in grey. The times of primary (P) and secondary (S) eclipses are indicated on top. 
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.
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. 
For reference, the minimum frequency (also spacing) is given by the time span of observations, f_{Δ} = 1/Δ = 0.045 c d^{−1}, which means that frequencies differing by f_{Δ} are certainly possible to distinguish. The maximum frequency (also Nyquist) is given by the sampling, f_{Ny} = 1/(2δ ≐ 1000 c d^{−1}. The orbital frequency of the binary Aa1+Aa2 is f_{orb} = 1/P_{1} = 0.174 c d^{−1}.
In the MOST data, the first two frequencies are f_{1} = 0.218 c d^{−1}, f_{2} = 0.478 c d^{−1}. These correspond to the Chandra data (Nichols et al. 2015), f_{X1} = 0.210 c d^{−1}, f_{X2} = 0.490 c d^{−1}, within uncertainties. The first one was identified as the rotation 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, f_{rot3} = v_{rot3}/(2πR_{3} sin i_{2}) = 0.474 c d^{−1}, 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.
Alternatively, when we assume an asynchronous primary, the synchronicity is F_{1} = f_{1}/f_{orb} = 1.250. For comparison, the periastronsynchronised value is F_{1} = [(1 + e_{1})/(1 − e_{1})^{3}]^{1/2} = 1.181. In both cases, it means a minor modification of the LC.
The remaining frequencies (f_{3}, f_{4}, f_{10}) are likely associated with pulsations, namely loworder 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.
11. 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 noninteracting, 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:

The distance of the system was estimated from the new Gaia DR3 parallax of the faint component δ Ori (Ca+Cb).

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.

The secondary (Aa2) spectrum in the blue spectral region was detected by the two and threestep disentangling.

The RV curve of the secondary was obtained by crosscorrelation with a disentangled template spectrum.

The twobody model of the eclipsing binary was constructed in PHOEBE2.

The threebody model in Xitau was constrained by all observables.
Compared to previous studies, we obtained significantly lower masses than Pablo et al. (2015), in their lowmass 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).
11.1. 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 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 heliumrich, with Hα emission (Table 3). Detailed stellarevolution models with possible mass transfer between (some of) the components shall be computed in future work. Additionally, longbaseline optical or nearinfrared 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.
Fig. 20. Hertzsprung–Russel diagram with the positions of the Aa1, Aa2, and Ab components of δ Ori and evolutionary tracks from Paxton et al. (2015). Numbers next to the main sequence indicate the theoretical masses, small numbers the theoretical gravitational acceleration log g. According to the 3body 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. 
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ónDíaz et al. 2015). All of its components seem to be less evolved^{15}. 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.
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.
Niesten (1904) reported micrométriques measures by 15cm Merz refractor, with two position angles 162° and a note ‘en contact’. However, it is unlikely that it corresponds to δ Ori Ab, because its separation at that epoch was only 0.11″; separations in other binaries were 1″ or more. This observation is not compatible with our model, which indicates θ ≐ 127°.
The program HEC35D with User’s Manual is available at https://astro.troja.mff.cuni.cz/ftp/hec/HEC35/
We had some success using only a onestep method and setting Aa1, Aa2 to have fixed intensity factors s_{1} = 1, s_{2} = 1 and Ab to have free intensity factor s_{3}. Parameters T_{0}, e_{1}, ω_{1}, K_{1}, q_{1} were converged and P_{1}, were fixed. This setting prevents fluctuations related to variable s_{1}, s_{2} factors, making the model much more stable. Otherwise, the continuum is wavy when s_{1} is free. The results from this approach are e_{1} = 0.0804, ω_{1} = 153.9°, K_{1} = 108.3 km s^{−1}, q_{1} = 0.4893.
In the Fabrycky (2010) model, only the radial force component is included. In the multipole model (Brož et al. 2021), containing all components, the value of J_{2} = −C_{20} ≃ 1.5 × 10^{−4}, or equivalently k_{2} = J_{2}(Ω_{0}/n_{0})^{2} ≃ 2.5 × 10^{−3}.
Acknowledgments
This research was supported by the grants P209/10/0715, GA1502112S, and GA1901995S of the Czech Science Foundation. M.B. was supported by the Czech Science Foundation grant GA2111058S. We acknowledge the support from the Research Programmes MSM0021620860 “Physical study of objects and processes in the solar system and in astrophysics” of the Ministry of Education of the Czech Republic and the grant AV0Z10030501 of the Academy of Sciences of the Czech Republic. Thanks to Dr. Rachel Matson for providing the data from Washington Double Star Catalogue. H.B. acknowledges financial support from the Croatian Science Foundation under the project 6212 “Solar and Stellar Variability”. We gratefully acknowledge the use of the program KOREL written by Petr Hadrava, PYTERPOL by Jana Nemravová, and reSPEFO2 by Adam Harmanec. This research has made use of the AMBER data reduction package of the JeanMarie Mariotti Center (http://www.jmmc.fr/amberdrs). Our thanks are due to Š. Dvořáková, D. Korčáková, J. Kubát, J. Nemravová, M. Oksala, P. Škoda, and V. Votruba who obtained some of the Ondřejov blue spectra used in this study. We also acknowledge the use of the electronic database from CDS Strasbourg and the electronic bibliography maintained by the NASA/ADS system. CCL and GAW acknowledge Discovery Grant support by the Natural Science and Engineering Research Council (NSERC) of Canada. Last but not least, we would like to thank our anonymous referee for the constructive useful comments that helped to improve the paper.
References
 Allen, M. G., Ochsenbein, F., Derriere, S., et al. 2014, in Astronomical Data Analysis Software and Systems XXIII, eds. N. Manset, & P. Forshay, ASP Conf. Ser., 485, 219 [NASA ADS] [Google Scholar]
 Anderson, E., & Francis, C. 2012, Astron. Lett., 38, 331 [Google Scholar]
 Bagnuolo, W. G. Jr, Riddle, R. L., Gies, D. R., & Barry, D. J. 2001, ApJ, 554, 362 [NASA ADS] [CrossRef] [Google Scholar]
 Bessell, M. 2000, in Encyclopedia of Astronomy and Astrophysics, ed. P. Murdin, 1939 [Google Scholar]
 Božić, H. 1998, Hvar Obs. Bull., 22, 1 [Google Scholar]
 Breiter, S., & Vokrouhlický, D. 2015, MNRAS, 449, 1691 [NASA ADS] [CrossRef] [Google Scholar]
 Brož, M. 2017, ApJS, 230, 19 [Google Scholar]
 Brož, M., Marchis, F., Jorda, L., et al. 2021, A&A, 653, A56 [Google Scholar]
 Brož, M., Ďurech, J., Carry, B., et al. 2022a, A&A, 657, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brož, M., Harmanec, P., Zasche, P., et al. 2022b, A&A, 666, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown, A. G. A., de Geus, E. J., & de Zeeuw, P. T. 1994, A&A, 289, 101 [NASA ADS] [Google Scholar]
 Burssens, S., SimónDíaz, S., Bowman, D. M., et al. 2020, A&A, 639, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carroll, K., Zee, R., & Matthews, J. 1998, The MOST Microsatellite Mission: Canada’s First Space Telescope [Google Scholar]
 Conroy, K. E., Kochoska, A., Hey, D., et al. 2020, ApJS, 250, 34 [Google Scholar]
 Corcoran, M. F., Nichols, J. S., Pablo, H., et al. 2015, ApJ, 809, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
 Cutri, R. M., Wright, E. L., & Conrow, T. 2012, VizieR Online Data Catalog: II/311 [Google Scholar]
 Ducati, J. R. 2002, VizieR Online Data Catalog: Catalogue of Stellar Photometry in Johnson’s 11color system [Google Scholar]
 Egan, M. P., Price, S. D., Kraemer, K. E., et al. 2003, VizieR Online Data Catalog: V/114 [Google Scholar]
 Fabrycky, D. C. 2010, in Exoplanets, ed. S. Seager, 217 [Google Scholar]
 Fitzpatrick, R. 2012, An Introduction to Celestial Mechanics (UK: Cambridge University Press) [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Gaia Collaboration 2020, VizieR Online Data Catalog: I/350 [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, in press, https://doi.org/10.1051/00046361/202243940 [Google Scholar]
 Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Hadrava, P. 1995, A&AS, 114, 393 [NASA ADS] [Google Scholar]
 Hadrava, P. 1997, A&AS, 122, 581 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hadrava, P. 2004, Publ. Astron. Inst. Acad. Sci. Czech Rep., 92, 15 [Google Scholar]
 Hadrava, P. 2005, Ap&SS, 296, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Harmanec, P. 1988, Bull. Astron. Inst. Czech., 39, 329 [NASA ADS] [Google Scholar]
 Harmanec, P., Mayer, P., & Šlechta, M. 2013, Massive Stars: From alpha to Omega, 70 [Google Scholar]
 Harvey, A. S., Stickland, D. J., Howarth, I. D., & Zuiderwijk, E. J. 1987, The Observatory, 107, 205 [NASA ADS] [Google Scholar]
 Harvin, J. A., Gies, D. R., Bagnuolo, W. G. Jr., Penny, L. R., & Thaller, M. L. 2002, ApJ, 565, 1216 [NASA ADS] [CrossRef] [Google Scholar]
 Heintz, W. D. 1980, ApJS, 44, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Hill, G. 1982, Publ. Dominion Astrophys. Obs. Victoria, 16, 67 [NASA ADS] [Google Scholar]
 Hummel, C. A., White, N. M., Elias, N. M. I., Hajian, A. R., & Nordgren, T. E. 2000, ApJ, 540, L91 [NASA ADS] [CrossRef] [Google Scholar]
 Ishihara, D., Onaka, T., Kataza, H., et al. 2010, A&A, 514, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kaufer, A., Stahl, O., Tubbesing, S., et al. 1999, The Messenger, 95, 8 [Google Scholar]
 Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanz, T., & Hubený, I. 2003, ApJS, 146, 417 [NASA ADS] [CrossRef] [Google Scholar]
 Lanz, T., & Hubený, I. 2007, ApJS, 169, 83 [CrossRef] [Google Scholar]
 Leahy, D. A., & Leahy, J. C. 2015, Comput. Astrophys. Cosmol., 2, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Lenz, P., & Breger, M. 2004, in The AStar Puzzle, eds. J. Zverko, J. Ziznovsky, S. J. Adelman, & W. W. Weiss, 224, 786 [Google Scholar]
 Leone, F., Bohlender, D. A., Bolton, C. T., et al. 2010, MNRAS, 401, 2739 [NASA ADS] [CrossRef] [Google Scholar]
 Maíz Apellániz, J., Trigueros Páez, E., Negueruela, I., et al. 2019, A&A, 626, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martins, F., Schaerer, D., & Hillier, D. J. 2005, A&A, 436, 1049 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mason, B. D., Martin, C., Hartkopf, W. I., et al. 1999, AJ, 117, 1890 [NASA ADS] [CrossRef] [Google Scholar]
 Mason, B. D., Wycoff, G. L., Hartkopf, W. I., Douglass, G. G., & Worley, C. E. 2001, AJ, 122, 3466 [Google Scholar]
 Maucó, K., Briceño, C., Calvet, N., et al. 2018, ApJ, 859, 1 [CrossRef] [Google Scholar]
 Mayer, P., Harmanec, P., Wolf, M., Božić, H., & Šlechta, M. 2010, A&A, 520, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moultaka, J., Ilovaisky, S. A., Prugniel, P., & Soubiran, C. 2004, PASP, 116, 693 [NASA ADS] [CrossRef] [Google Scholar]
 Nelder, J. A., & Mead, R. 1965, Comput. J., 7, 308 [Google Scholar]
 Nemravová, J. A., Harmanec, P., Brož, M., et al. 2016, A&A, 594, A55 [Google Scholar]
 Neugebauer, G., Habing, H. J., van Duinen, R., et al. 1984, ApJ, 278, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Nichols, J., Huenemoerder, D. P., Corcoran, M. F., et al. 2015, ApJ, 809, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Niesten, L. 1904, Annales de l’Observatoire Royal de Belgique Nouvelle serie, 8, 1 [Google Scholar]
 Oplištilová, A., Harmanec, P., Mayer, P., et al. 2020, Contrib. Astron. Obs. Skalnate Pleso, 50, 585 [Google Scholar]
 Pablo, H., Richardson, N. D., Moffat, A. F. J., et al. 2015, ApJ, 809, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Pablo, H., Whittaker, G. N., Popowicz, A., et al. 2016, PASP, 128 [Google Scholar]
 Parenago, P. P. 1954, Trudy Gosudarstvennogo Astronomicheskogo Instituta, 25, 3 [NASA ADS] [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Perryman, M. A. C., & ESA 1997, in The HIPPARCOS and TYCHO catalogues, (Noordwijk, Netherlands: ESA Publications Division), ESA SP Ser., 1200 [Google Scholar]
 Pigulski, A. 2018, in 3rd BRITE Science Conference, eds. G. A. Wade, D. Baade, J. A. Guzik, & R. Smolec, 8, 175 [NASA ADS] [Google Scholar]
 Prša, A., & Zwitter, T. 2005, ApJ, 628, 426 [Google Scholar]
 Puebla, R. E., Hillier, D. J., Zsargó, J., Cohen, D. H., & Leutenegger, M. A. 2016, MNRAS, 456, 2907 [NASA ADS] [CrossRef] [Google Scholar]
 Richardson, N. D., Moffat, A. F. J., Gull, T. R., et al. 2015, ApJ, 808, 88 [NASA ADS] [CrossRef] [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, eds. J. Oschmann, M. Jacobus, M. Clampin, G. G. Fazio, & H. A. MacEwen, SPIE Conf. Ser., 9143, 914320 [Google Scholar]
 Robert, C., & Casella, G. 2011, Stat. Sci., 26, 102 [Google Scholar]
 Rowan, N. 1990, PhD Thesis, Univ. Texas Austin, USA [Google Scholar]
 SchmidtKaler, T. 1982, Bulletin d’Information du Centre de Donnees Stellaires, 23, 2 [Google Scholar]
 Shabun, K., Richichi, A., Munari, U., Siviero, A., & Paczynski, B. 2008, in A Giant Step: from Milli to Microarcsecond Astrometry, eds. W. J. Jin, I. Platais, & M. A. C. Perryman, 248, 118 [NASA ADS] [Google Scholar]
 Shenar, T., Oskinova, L., Hamann, W. R., et al. 2015, ApJ, 809, 135 [Google Scholar]
 SimónDíaz, S., Caballero, J. A., Lorenzo, J., et al. 2015, ApJ, 799, 169 [CrossRef] [Google Scholar]
 Škoda, P., Šlechta, M., & Honsa, J. 2002, Publ. Astron. Inst. Czech. Acad. Sci., 90 [Google Scholar]
 Standish, E. M., & Williams, J. G. 2006, in Orbital Ephemerides of the Sun, Moon, and Planets, ed. P. Seidelmann (University Science Books) [Google Scholar]
 Tierney, L. 1994, Ann. Stat., 1701 [Google Scholar]
 van Hamme, W. 1993, AJ, 106, 2096 [Google Scholar]
 Van Hamme, W., & Wilson, R. E. 2007, ApJ, 661, 1129 [NASA ADS] [CrossRef] [Google Scholar]
 Walborn, N. R. 1972, AJ, 77, 312 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, D., Mizuno, D., Buffington, A., et al. 2006, J. Geophys. Res., 111 [CrossRef] [Google Scholar]
 Wilson, R. E. 1979, ApJ, 234, 1054 [Google Scholar]
 Wilson, R. E. 1990, ApJ, 356, 613 [Google Scholar]
 Wilson, R. E., & Devinney, E. J. 1971, ApJ, 166, 605 [Google Scholar]
 Wilson, R. E., Van Hamme, W., & Terrell, D. 2010, ApJ, 723, 1469 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Secular rates for δ Ori A
We used the standard Gauss equations to estimate the secular rates for the δ Ori A system. A perturbation due to the quadrupole moment J_{2} ≡ −C_{20} of the primary induces a precession of the argument of pericentre:
and the longitude of the node:
where , , and is the inclination with respect to the primary equator. Otherwise, the notation is the same as in Table 11. For , e_{1} → 0, J_{2} > 0, we would have , , .
The quadrupole moment is determined by the internal structure:
where the Legendre polynomial . For a homogeneous body, it would be related to the ellipticity (Fitzpatrick 2012) ^{16}:
For a rotating body, it is related to the Love number ^{17}:
where Ω_{0} is the angular rotation frequency; n_{0}, the mean motion at the surface. Assuming J_{2} = 1.8 ⋅ 10^{−4} results in , , where the node circulates with respect to the primary equator, but it only librates with respect to the observer plane.
For the binary Aa1+Aa2 (acting on a massless particle), the equations should be modified as follows:
where the inclination is with respect to the binary. The effective quadrupole moment is:
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 , denote the reduced masses; , , 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, , , which is negligible on the observational time span. All these effects were nevertheless included in our numerical Nbody model (Sect. 9).
Appendix B: Supplementary figures and tables
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 threebody 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.
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 largeamplitude oscillations (not included in our model); uncertainties 0.01 mag were thus assigned to all data points. 
Fig. B.2. Contributions to χ^{2} for a set 100 bestfit models of δ Ori. Individual contributions are shown in the panels (from top left): astrometry (SKY), RV, eclipse timings (ETV), eclipse duration (ECL), normalised spectra (SYN), SED, light curve (LC), and total. 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. 
Details on spectra from the coudé focus of the Ondřejov 2m reflector in the blue region. The RVs were determined during the spectral disentangling (threestep method) in KOREL.
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 (threestep method) in KOREL.
Results of eight phoebe2 models of δ Ori. LCs from eight individual seasons S. (Table 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 Table 9.
All Tables
Atmospheric parameters derived with PYTERPOL from the blue spectral region 4271–4513 Å, with the mean resolution of 0.0144 Å.
Results of fitting three phoebe2 models for δ Ori Aa1+Aa2. LC from season 2016 and all RVs were used to constrain the models.
Details on spectra from the coudé focus of the Ondřejov 2m reflector in the blue region. The RVs were determined during the spectral disentangling (threestep method) in KOREL.
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 (threestep method) in KOREL.
Results of eight phoebe2 models of δ Ori. LCs from eight individual seasons S. (Table 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 Table 9.
All Figures
Fig. 1. Scheme of the multiple system δ Ori (HD 36486, ADS 4134, Mintaka). Orbital periods were taken from ^{*}Leone et al. (2010), ^{**}Mayer et al. (2010), and ^{* * *}this paper. 

In the text 
Fig. 2. Coverage of the orbital period P_{1} by blue spectra. Phases are determined with respect to time T_{0} = HJD 2454002.8735 (time of periastron passage determined by KOREL) for the eclipsing binary. 

In the text 
Fig. 3. Photometric data from MOST and BRITE displayed with respect to time. The BRITE data covers six consecutive seasons between 2013 and 2021. 

In the text 
Fig. 4. Sorted distances of objects from the Orion OB1b (blue) and OB1a (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. 

In the text 
Fig. 5. Orbit of 2body model. Top: orbit of (Aa1+Aa2)+Ab components in the (u–v) plane (green), calculated using the twobody 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 Table 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}. 

In the text 
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 twobody 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. 

In the text 
Fig. 7. Pseudoχ^{2} in Fourier space vs. mass ratio q_{1} = m_{2}/m_{1}. This is related to the twostep disentangling method, to its second step, when the signal of the secondary component (Aa2) was sought for in the residuals. 

In the text 
Fig. 8. Comparison of the disentangled spectra (blue) of the Aa1, Aa2, and Ab components with the bestfit 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 Table 7). 

In the text 
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 (forestgreen greenblack) and the secondary Aa2 (pypurple purpleblack). The grey points correspond to our model, the red lines to the residuals, or contributions to χ^{2}. 

In the text 
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 Table 9). Each subdiagonal panel shows a 2D histogram, the isolines corresponding to the confidence intervals, and the correlations between parameters. 

In the text 
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 burnin. 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 γ. 

In the text 
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). 

In the text 
Fig. 13. Relative differences of parameters derived for δ Ori for eight seasons (Table 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. Explanations of the parameters can be found in Table 9. 

In the text 
Fig. 14. Model from SED data. Top: comparison of the observed (blue) and synthetic (orange) SED of δ Ori. The residuals are plotted in redredblack. 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. 

In the text 
Fig. 15. Example of a convergence of the threebody 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 nonunit weights (see values in the main text). 

In the text 
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. 

In the text 
Fig. 17. Same as Fig. 16, with the distinguished contributions of individual components Aa1, Aa2, Ab of δ Ori. The wavelength range is limited to 430–440 nm. Only the sixth spectrum, 2455836.582700 is shown. 

In the text 
Fig. 18. Residuals of the MOST LC (red), and the Fourier series up to the order of ℓ = 10 (black; Table 13). Individual orders as they are summed up are shown in grey. The times of primary (P) and secondary (S) eclipses are indicated on top. 

In the text 
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. 

In the text 
Fig. 20. Hertzsprung–Russel diagram with the positions of the Aa1, Aa2, and Ab components of δ Ori and evolutionary tracks from Paxton et al. (2015). Numbers next to the main sequence indicate the theoretical masses, small numbers the theoretical gravitational acceleration log g. According to the 3body 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. 

In the text 
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 largeamplitude oscillations (not included in our model); uncertainties 0.01 mag were thus assigned to all data points. 

In the text 
Fig. B.2. Contributions to χ^{2} for a set 100 bestfit models of δ Ori. Individual contributions are shown in the panels (from top left): astrometry (SKY), RV, eclipse timings (ETV), eclipse duration (ECL), normalised spectra (SYN), SED, light curve (LC), and total. 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. 

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