© The Authors 2023

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)¹ is a detached eclipsing binary with a negligible mass transfer, the orbital period P₁  =  5.732436 d (Mayer et al. 2010), a slightly eccentric orbit (0.08), and apsidal motion (1.45°  yr⁻¹) (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₂ must be of the order of tens of thousands of days.

• B (V_B ≃ 14 mag) is a very faint distant companion² 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.

• Ca+Cb (V_Ca + Cb = 6.85 mag) is another distant companion³ to (Aa1+Aa2)+Ab that is a spectroscopic, non-eclipsing 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 sub-system δ 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₁ = 11.2 M_⊙ and m₂ = 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₂/m₁.

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 X-ray emission line widths. They identified two periods in the total X-ray flux: 4.76 ± 0.30 and 2.04 ± 0.05 days. 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⁻⁷ M_⊙ yr⁻¹. 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.

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.

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 φ₁ is illustrated in Fig. 2. The short period P₁ 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 reSPEFO2⁴ (written by Adam Harmanec).

Fig. 2. Coverage of the orbital period P1 by blue spectra. Phases are determined with respect to time T0 = HJD 2454002.8735 (time of periastron passage determined by KOREL) for the eclipsing binary.

2.2. 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).

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 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.

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⁻³.

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)⁵. The OB1a subgroup (north-west) 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.

4. Visual orbit of (Aa1+Aa2)+Ab

To determine the parameters of the long-period (P₂) orbit of the (Aa1+Aa2)+Ab system, we constructed a simplified two-body model in the Xitau program (Brož 2017; Brož et al. 2021, 2022a,b)⁶. 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⁷, 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 γ₁ of δ Ori (Harvey et al. 1987; Harvin et al. 2002), which vary between approximately 12 and 23 km s⁻¹ (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 γ₁, that is, blending with Ab (1910, 1948), low amplitude of RV curve K₁ (1951, 1969, 1981), and different γ₁ 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 χ² = 95, with contributions ${\chi }_{\text{SKY}}^2=60$ for astrometry and ${\chi }_{\text{RV}}^2=35$ for RVs. Although χ² > 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 γ₁.

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 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 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 P2.

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.

We estimated the uncertainties of parameters using χ² mapping and verified the MCMC method. According to the χ² statistics, a 1σ level corresponds to χ² ≐ 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₂. 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 χ², especially the ${\chi }_{\text{rv}}^2$ contribution. Moreover, in the mirror solution, the RVs of Ab are larger than γ₁, 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 long-period (P₂) 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 χ² 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 χ² sum is almost comparable to the noise. The mass ratio q₁ 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 (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.⁸ 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 line-free regions and assigned each spectrum a weight w according to the formula:

$$\begin{array}{c}\hfill w=\frac{{(S/N)}^2}{{(S/N)}_{\mathrm{mean}}^2},\end{array}$$(1)

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:

$$\begin{array}{c}\hfill \mathcal{F}{I}_i(y,t)={\displaystyle \sum _j{s}_{\mathit{ij}}\mathcal{F}{I}_{\mathit{ij}}·exp(i{\mathit{yv}}_{\mathit{ij}})\mathrm{for}\forall i,}\end{array}$$(2)

where j denotes the component; i, the spectrum; I, the normalised intensity; ℱI, its Fourier transform; y, the Fourier-transformed quantity x ≡ lnλ/λ₀, related to the wavelength λ; s_ij, the intensity factors (constant or variable); v_ij, the radial velocity.

5.1. 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₁, e₁, ω₁, K₁, T₀, while P₁, $\dot{\omega }$ 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, s₃ = 1) constant and assigned zero intensity to Aa2 (s₂ = 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₁ = 0, s₃ = 0) and the intensity factor of Aa2 was constant s₂ = 1. We fitted the spectrum of the Aa2 component by converging the mass ratio q and fixing T₀, e₁, ω₁, K₁, P₁, and ${\dot{\omega }}_1$. We successfully detected the desired spectrum of the Aa2 component. The determined parameters are summarised in Table 6. This method gave higher e₁, lower q₁, precise K₁ and ω₁, 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-χ² map (see Fig. 7).⁹; uncertainties as in Table 6.

Fig. 7. Pseudo-χ2 in Fourier space vs. mass ratio q1 = m2/m1. 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.

5.2. 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₁, s₃. We fixed P₁ and $\dot{{\omega }_1}$ of the close orbit and converged T₀, e₁, ω₁, K₁. 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₁, all parameters were fixed. Finally, we found the solution for all three components with constant s-factors, free T₀, e₁, ω₁, K₁, q₁, and fixed P₁, $\dot{{\omega }_1}$.

The resulting s₁(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₁ = 0.4517 is higher, compared to the two-step disentangling, while e₁ = 0.0761 is close to that found from the LC, and K₂ = 239.7 km s⁻¹.

In the three-step 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)¹⁰. 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 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 Table 7).

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¹¹, 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 χ² 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:

$$\begin{array}{c}\hfill P({\mathbf{\Theta }}_{\mathrm{M}}|\mathit{D})=\frac{P(\mathit{D}|{\mathbf{\Theta }}_{\mathrm{M}})P({\mathbf{\Theta }}_{\mathrm{M}})}{P(\mathit{D})},\end{array}$$(3)

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 χ²; 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:

$$\begin{array}{c}\hfill F_{c,\mathrm{norm}}=S_{\mathrm{c}}·\frac{F_{c,i}}{\max (F_{\mathrm{c},i})},\end{array}$$(4)

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 $\dot{\omega }=1.45°{\mathrm{y}}^{-1}$ (Pablo et al. 2015), and in some models also the effective temperature of the primary T₁ = 31 000 K, and the third light, additional to the components Aa1 and Aa2, l₃ = 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.

7.1. 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 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 χ². 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 T1 and the third light l3 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.

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₁ 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₁ and m_i show negative correlations. The value of the systemic velocity γ is a little problematic as the value of 21.96 km s⁻¹ was assumed and subtracted, then our model drifted to about −2.5 km s⁻¹, so that the resulting value is 18.5 km s⁻¹.

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 sub-diagonal 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 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 γ.

The detached binary system Aa1+Aa2 is shown in Fig. 12. In addition, we derived several parameters from the nominal phoebe2 model (χ² = 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⁻²]):

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{V}}={\mathrm{\Delta }}_{\mathrm{eff}}\omega {\displaystyle \sum _k{I}_{\lambda k}{S}_k{\mu }_k{\eta }_k,}\end{array}$$(5)

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).

where Δ_eff [m^−1] stands for the effective wavelength range; ω = 1 m²/d² [sr⁻¹], for the solid angle; d, for the distance of the system; ∑_k, for the summation over the triangular elements (grid); I_λ [W m⁻³ sr⁻¹], the monochromatic intensity on the stellar surface; S [m²], 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):

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\lambda ,\mathrm{cal}}={\mathrm{\Phi }}_{\nu ,\mathrm{cal}}\frac{c}{{\lambda }_{\mathrm{eff}}^2}=0.03669877\mathrm{W}{\mathrm{m}}^{-3}\end{array}$$(6)

and the Johnson V passband flux:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{V},\mathrm{cal}}={\displaystyle {\int }_{\lambda }f{\mathrm{\Phi }}_{\lambda ,\mathrm{cal}}\mathrm{d}\lambda \doteq {\mathrm{\Delta }}_{\mathrm{eff}}{\mathrm{\Phi }}_{\lambda ,\mathrm{cal}}\doteq 3.119396\times {10}^{-9}\mathrm{W}{\mathrm{m}}^{-2},}\end{array}$$(7)

where f(λ) denotes the filter transmission; λ_eff = 545 nm, the effective wavelength; and Δ_eff = 85 nm, the effective range.

The apparent magnitude V₀ (without absorption) of the primary Aa1 and Aa2 is then:

$$\begin{array}{c}\hfill V_0=0[\mathrm{mag}]-2.5log\frac{{\mathrm{\Phi }}_{\mathrm{V}}}{{\mathrm{\Phi }}_{\mathrm{V},\mathrm{cal}}}.\end{array}$$(8)

For the tertiary, we used the value of the third light:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{V},2}=({\mathrm{\Phi }}_{\mathrm{V},1}+{\mathrm{\Phi }}_{\mathrm{V},2})l_{3\mathrm{V}}.\end{array}$$(9)

Comparing V₀ 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₀ 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₁ and the third light l₃ 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 T1 of the primary (Aa1) and the third light l3 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)¹², 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 line-of-sight samples in the maps of Green et al. (2019)¹³.

From the photometric catalogues in http://vizier.u-strasbg.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 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).

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 (far-infrared). 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. 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:

$$\begin{array}{c}\hfill {\mathit{f}}_i={\displaystyle \sum _{j\ne i}\frac{G{m}_j}{r_{\mathit{ij}}^3}{\mathit{r}}_{\mathit{ij}}+{\mathit{f}}_{\mathrm{oblat}}+{\mathit{f}}_{\mathrm{ppn}}\text{for}\forall i,}\end{array}$$(10)

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 Ω₁, ϖ₁, Ω₂, ϖ₂, 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)¹⁴, which results in the observed value of precession $\dot{\omega }\simeq 1.45°{\mathrm{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⁻⁸, 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 χ² metrics were multiplied by weights:

$$\begin{array}{cc}& {\chi }^2=w_{\mathrm{SKY}}{\chi }_{\mathrm{SKY}}^2+w_{\mathrm{RV}}{\chi }_{\mathrm{RV}}^2+w_{\mathrm{ETV}}{\chi }_{\mathrm{ETV}}^2+w_{\mathrm{ECL}}{\chi }_{\mathrm{ECL}}^2+w_{\mathrm{LC}}{\chi }_{\mathrm{LC}}^2+\hfill \\ \hfill & +w_{\mathrm{SYN}}{\chi }_{\mathrm{SYN}}^2+w_{\mathrm{SED}}{\chi }_{\mathrm{SED}}^2,\hfill \end{array}$$(11)

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₃, T₃, log g₃). 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 i) the synthetic spectra (in particular, log g₂ or R₂) 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 χ² 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 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).

The observed and synthetic SEDs are compared in Fig. 14. Even though the corresponding contribution ${\chi }_{\text{LC}}^2$ 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₂ = 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₂, 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₁ + m₂ + m₃).

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₁ + m₂ + m₃ = 52.0 M_⊙ and performed a similar testing as in Fig. B.2. The overall best-fit has χ² = 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 (${\chi }_{\text{SYN}}^2$ = 69 441 vs. 44 795). This is directly related to the log g₃ 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 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), 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 fmost and its combinations with f1, etc.

For reference, the minimum frequency (also spacing) is given by the time span of observations, f_Δ = 1/Δ = 0.045 c d⁻¹, 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⁻¹. The orbital frequency of the binary Aa1+Aa2 is f_orb = 1/P₁ = 0.174 c d⁻¹.

In the MOST data, the first two frequencies are f₁ = 0.218 c d⁻¹, f₂ = 0.478 c d⁻¹. These correspond to the Chandra data (Nichols et al. 2015), f_X1 = 0.210 c d⁻¹, f_X2 = 0.490 c d⁻¹, 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₃ sin i₂) = 0.474 c d⁻¹, where we assumed an alignment.

Interestingly, f_rot1 = v_rot1/(2πR₁ sin i₁) = 0.174 c d⁻¹ derived from rotational broadening is not equal to f₁, 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₁ = f₁/f_orb = 1.250. For comparison, the periastron-synchronised value is F₁ = [(1 + e₁)/(1 − e₁)³]^1/2 = 1.181. In both cases, it means a minor modification of the LC.

The remaining frequencies (f₃, f₄, f₁₀) 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.

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 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 cross-correlation 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₃ 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₃ value (3.2) or to its m₃ 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.7, or m₃ ≃ 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 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.

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 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 g3 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ón-Díaz et al. 2015). All of its components seem to be less evolved¹⁵. 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.

Acknowledgments

This research was supported by the grants P209/10/0715, GA15-02112S, and GA19-01995S of the Czech Science Foundation. M.B. was supported by the Czech Science Foundation grant GA21-11058S. 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 Jean-Marie 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

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₂ ≡ −C₂₀ of the primary induces a precession of the argument of pericentre:

$$\begin{array}{c}\hfill {\dot{\omega }}_1=+3n_1{J}_2{\left(\frac{R_1}{a_1}\right)}^2\frac{5{cos}^2{\stackrel{\sim }{i}}_1-1}{4{\eta }_1^4},\end{array}$$(A.1)

and the longitude of the node:

$$\begin{array}{c}\hfill {\dot{\mathrm{\Omega }}}_1=-\frac{3}{2}{n}_1{J}_2{\left(\frac{R_1}{a_1}\right)}^2\frac{cos{\stackrel{\sim }{i}}_1}{{\eta }_1^4},\end{array}$$(A.2)

where $n_1=\sqrt{G(m_1+m_2)/a_1^3}$, ${\eta }_1\equiv \sqrt{1-e_1^2}$, and ${\stackrel{\sim }{i}}_1$ is the inclination with respect to the primary equator. Otherwise, the notation is the same as in Table 11. For ${\stackrel{\sim }{i}}_1\to 0$, e₁ → 0, J₂ > 0, we would have ${\dot{\omega }}_1>0$, ${\dot{\mathrm{\Omega }}}_1<0$, ${\dot{\omega }}_1=-2{\dot{\mathrm{\Omega }}}_1$.

The quadrupole moment is determined by the internal structure:

$$\begin{array}{c}\hfill J_2=-\frac{1}{m_1{R}_1^2}{\displaystyle {\int }_V\rho {|\mathit{r}|}^2{P}_2(cos\theta )\mathrm{d}V,}\end{array}$$(A.3)

where the Legendre polynomial $P_2(x)=\frac{1}{2}(3x^2-1)$. For a homogeneous body, it would be related to the ellipticity (Fitzpatrick 2012) ¹⁶:

$$\begin{array}{c}\hfill J_2=-\frac{2}{5}{\epsilon }_0.\end{array}$$(A.4)

For a rotating body, it is related to the Love number ¹⁷:

$$\begin{array}{c}\hfill J_2\simeq k_2{\left(\frac{{\mathrm{\Omega }}_0}{n_0}\right)}^2,\end{array}$$(A.5)

where Ω₀ is the angular rotation frequency; n₀, the mean motion at the surface. Assuming J₂ = 1.8 ⋅ 10⁻⁴ results in $\dot{\omega }=1.47\mathrm{deg}{\mathrm{y}}^{-1}$, $\dot{\mathrm{\Omega }}=-0.73\mathrm{deg}{\mathrm{y}}^{-1}$, 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 mass-less particle), the equations should be modified as follows:

$$\begin{array}{c}\hfill {\dot{\omega }}_2=+3n_2{J}_2{\left(\frac{a_1}{a_2}\right)}^2\frac{5{cos}^2{\stackrel{\sim }{i}}_2-1}{4{\eta }_2^4},\end{array}$$(A.6)

$$\begin{array}{c}\hfill {\dot{\mathrm{\Omega }}}_2=-\frac{3}{2}{n}_2{J}_2{\left(\frac{a_1}{a_2}\right)}^2\frac{cos{\stackrel{\sim }{i}}_2}{{\eta }_2^4},\end{array}$$(A.7)

where the inclination ${\stackrel{\sim }{i}}_2$ is with respect to the binary. The effective quadrupole moment is:

$$\begin{array}{c}\hfill J_2=\frac{1}{2}\frac{m_1{m}_2}{{(m_1+m_2)}^2},\end{array}$$(A.8)

because the respective radius (|r| = a₁) is the same as the reference radius (R = a₁). Given that J₂ = 0.109, and the ratio of a₁/a₂ = 0.002, the precession rates should be of the order of 10⁻⁶ deg y⁻¹.

However, for the massive triple system (Aa1+Aa2)+Ab, we used the theory of Breiter & Vokrouhlický (2015); for the longitudes of pericentre and node:

$$\begin{array}{cc}\hfill {\dot{\varpi }}_2=& \frac{3}{8}\frac{n_2}{{\eta }_2^3}\frac{m_3}{m_1+m_2+m_3}\frac{n_2}{n_1}\gamma ·\hfill \\ \hfill & [3{cos}^{2}J-1-\frac{\gamma sinJsin2J}{1+\gamma cosJ+\sqrt{1+{\gamma }^2+2\gamma cosJ}}],\hfill \end{array}$$(A.9)

$$\begin{array}{c}\hfill {\dot{\mathrm{\Omega }}}_2=\frac{3}{4}\frac{n_2}{{\eta }_2^3}\frac{m_3}{m_1+m_2+m_3}\frac{n_2}{n_1}cosJ\sqrt{1+{\gamma }^2+2\gamma cosJ},\end{array}$$(A.10)

where $m_1^{\prime }\text{}=m_1{m}_2/(m_1+m_2)$, $m_2^{\prime }\text{}=(m_1+m_2)m_3/(m_1+m_2+m_3)$ denote the reduced masses; $L_1=m_1^{\prime }{n}_1{a}_1^2$, $L_2=m_2^{\prime }{n}_2{a}_2^2$, the angular momenta; cos J = cos i₁ cos i₂ + sin i₁ sin i₂ cos(Ω₁ − Ω₂), the mutual inclination; and γ = L₁/(L₂η₂). Assuming parameters from Table 11, ${\dot{\varpi }}_2=3.4·{10}^{-6}\mathrm{deg}{\mathrm{y}}^{-1}$, ${\dot{\mathrm{\Omega }}}_2=8.5·{10}^{-5}\mathrm{deg}{\mathrm{y}}^{-1}$, which is negligible on the observational time span. All these effects were nevertheless included in our numerical N-body 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 three-body model (Sect. 9) is shown in Fig. B.1. Individual contributions to χ² 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 large-amplitude 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 best-fit 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 g1, log g3, 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 g1, log g3, 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 g1, log g3 due to the SYN, SED, or large log g1 due to the ECL, LC. The weighted very best fit is denoted by the red circle.

All Tables

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 P1 by blue spectra. Phases are determined with respect to time T0 = 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 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 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 P2.

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

	Fig. 7. Pseudo-χ2 in Fourier space vs. mass ratio q1 = m2/m1. 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.
In the text

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 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 T1 and the third light l3 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 sub-diagonal 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 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 γ.
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 T1 of the primary (Aa1) and the third light l3 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 (far-infrared). 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 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).
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 fmost and its combinations with f1, 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 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 g3 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 large-amplitude 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 best-fit 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 g1, log g3, 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 g1, log g3, 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 g1, log g3 due to the SYN, SED, or large log g1 due to the ECL, LC. The weighted very best fit is denoted by the red circle.

In the text